# Optimal Transport Approach to Sobolev Regularity of Solutions to the Weighted Least Gradient Problem

Samer Dweik, Wojciech Górny

► To cite this version:

Samer Dweik, Wojciech Górny. Optimal Transport Approach to Sobolev Regularity of Solutions to the Weighted Least Gradient Problem. SIAM Journal on Mathematical Analysis, 2023, 55 (3), 10.1137/21M1468358. hal-04258091

HAL Id: hal-04258091

https://hal.science/hal-04258091v1

Submitted on 25 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# OPTIMAL TRANSPORT APPROACH TO SOBOLEV REGULARITY OF SOLUTIONS TO THE WEIGHTED LEAST GRADIENT PROBLEM

SAMER DWEIK, WOJCIECH GÓRNY

**ABSTRACT.** We study the equivalence between the weighted least gradient problem and the weighted Beckmann minimal flow problem or equivalently, the optimal transport problem with Riemannian cost. Thanks to this equivalence, we prove existence and uniqueness of a solution to the weighted least gradient problem. Then, we show $L^p$ regularity on the transport density between two singular measures in the corresponding equivalent Riemannian optimal transport formulation. This will imply $W^{1,p}$ regularity of the solution of the weighted least gradient problem.

## 1. INTRODUCTION

In this paper, we study a variant of the least gradient problem. The classical version of the least gradient problem consists of minimizing the total variation of the vector measure $Du$ among all BV functions $u$ defined on an open bounded domain $\Omega$ with given boundary datum $g$:

$$ (1.1) \qquad \inf \left\{ \int_{\Omega} |Du| : u \in BV(\Omega), u|_{\partial\Omega} = g \right\}, $$

where $u|_{\partial\Omega}$ denotes the trace of $u$ on the boundary $\partial\Omega$. It was first considered in this form in [38], where the authors prove existence and uniqueness of a solution to (1.1) in the case where $\Omega$ is strictly convex and the boundary datum $g$ is continuous. The convexity assumption on the domain is crucial for the analysis performed in [38], and indeed the authors provide examples of boundary data for which there exist no solutions if it is violated. In the case when $\Omega$ is not strictly convex, one needs to introduce admissibility conditions on the boundary data to obtain existence and uniqueness of solutions, see [14, 34]. On the other hand, on strictly convex domains the continuity assumption on the boundary data can be relaxed; for instance, in [17] the author proved that in two dimensions Problem (1.1) admits a solution as soon as $g \in BV(\partial\Omega)$ (for generalisations to higher dimensions, see [19, 32]). This result is sharp in the sense that a counter-example given in [37] shows that if $g \notin BV(\partial\Omega)$, solutions to Problem (1.1) might not exist. Moreover, we lose uniqueness of the solution to Problem (1.1) if $g \notin C(\partial\Omega)$.

In [15, 21], the authors proved that in 2D, the problem (1.1) is closely related to the following minimal flow formulation (which is called the Beckmann problem; see [4])

$$ (1.2) \qquad \inf \left\{ \int_{\Omega} |v| : v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2), \text{ div } v = 0 \text{ and } v \cdot n = f := \partial_{\tau} g \text{ on } \partial\Omega \right\}, $$

provided that the domain $\Omega$ is convex and $g \in BV(\partial\Omega)$. Here, the finite signed measure $\partial_{\tau}g$ is the tangential derivative of $g$. Moreover, $\mathcal{M}(\bar{\Omega}, \mathbb{R}^2)$ denotes the set of vector measures over $\bar{\Omega}$ while the divergence constraint is understood in the distributional sense. In fact, the infima in the two problems coincide, and if $\Omega$ is strictly convex there is a one-to-one correspondence between vector measures $Du$ which are derivatives of solutions to (1.1) and vector measures $v$ which solve (1.2). This equivalence is formally given by $v := R_{\frac{1}{2}} Du$ on $\bar{\Omega}$. If $\Omega$ is only convex, the one-to-one correspondence is still true provided that for every optimal flow $v$ of the Beckmann problem, its total variation $|v|$ gives zero mass to the boundary. As a consequence, regularity estimates for one problem can be transferred to the other one, in particular the $L^p$ summability of the optimal flow $v$ yields $W^{1,p}$ regularity of the solution $u$ to Problem (1.1). We will describe this equivalence in more detail in Section 2.3.

<footnote>Date: December 30, 2021.</footnote>

<footnote n="2020">2020 Mathematics Subject Classification. 35J25, 35J75, 49Q22.</footnote>

<footnote>Key words and phrases. Least Gradient Problem, Optimal transport, Anisotropy.</footnote>
On the other hand, the Beckmann problem (1.2) is equivalent to the optimal transport problem with Euclidean cost (which is called the Monge-Kantorovich problem; see [24,30]) with source and target measures located on $∂Ω$ (see [36]):

$$ (1.3) \quad \min \left\{ \int_{\Omega \times \bar{\Omega}} |x-y| d\Lambda : \Lambda \in \mathcal{M}^+(\bar{\Omega} \times \bar{\Omega}), (\Pi_x)_\# \Lambda = f^+ \text{ and } (\Pi_y)_\# \Lambda = f^- \right\}, $$

where $f^+$ and $f^-$ are the positive and negative parts of $f$ (we note that $f^+$ and $f^-$ have the same total mass since $f$ is the tangential derivative of $g$ on the closed set $∂Ω$, so $f(∂Ω) = 0$). Moreover, it is well known (see, for instance, [36]) that the Monge-Kantorovich problem (1.3) has a dual formulation, which is the following:

$$ (1.4) \quad \sup \left\{ \int_{\bar{\Omega}} \psi \, d(f^+ - f^-) : \psi \in \text{Lip}_1(\bar{\Omega}) \right\}. $$

From [36], every optimal flow $v$ of Problem (1.2) is of the form $v = -|v|∇ψ$, where $ψ$ is a Kantorovich potential (i.e. a maximizer for the dual problem (1.4)). In addition, one can show that this nonnegative measure $|v|$ is a transport density between $f^+$ and $f^-$, i.e. there is an optimal transport plan $Λ$ such that

$$ |v|(A) = \int_{\bar{\Omega} \times \bar{\Omega}} \mathcal{H}^1([x, y] \cap A) \, d\Lambda(x, y), \text{ for all Borel set } A \subset \bar{\Omega}. $$

We see that $|v|(∂Ω) = 0$ as soon as $Ω$ is strictly convex and so, the equivalence between Problems (1.1) & (1.2) gives existence of a solution $u$ for Problem (1.1). Moreover, classical results imply that the transport density $|v|$ is unique (i.e., it does not depend on the choice of the optimal transport plan $Λ$) provided $f^+$ or $f^-$ is absolutely continuous with respect to the Lebesgue measure $L^d$. The $L^p$ summability of the transport density has been studied in several papers [9–11,35] and the results can be summarized as follows: for every $p \in [1, \infty]$, the transport density is in $L^p(Ω)$ provided that both $f^+$ and $f^-$ are in $L^p(Ω)$. However, none of these assumptions is satisfied in our setting: the source and target measures $f^+$ and $f^-$ in Problem (1.3) are concentrated on the boundary, so it is not clear whether the transport density $|v|$ is unique or not. Nonetheless, the authors of [15] prove uniqueness of the optimal transport plan $Λ$ (and so, uniqueness of the optimal flow $v$) under the assumption that $f^+$ or $f^-$ is atomless and $Ω$ is strictly convex. Moreover, they proved that if the domain is uniformly convex and both $f^+$ and $f^-$ are in $L^p(∂Ω)$ with $p \le 2$, then the transport density between them is in $L^p(Ω)$. Recalling again the equivalence between Problems (1.1) & (1.2), this implies that on uniformly convex domains the solution $u$ of the problem (1.1) is unique as soon as $g \in C(∂Ω)$ and (for $p \le 2$) it belongs to $W^{1,p}(Ω)$ as soon as $g \in W^{1,p}(∂Ω)$.

In this paper, we are interested in studying the planar weighted least gradient problem (see [23,31]):

$$ (wLGP) \quad \inf \left\{ \int_{\Omega} k(x) |Du| : u \in BV(\Omega), u|_{\partial\Omega} = g \right\} $$

where $k: Ω \to \mathbb{R}^+$ is a given smooth function. It is closely connected to the conductivity imaging problem (see for instance [23]), which is an inverse problem appearing with relation to medical imaging, and its purpose is to recover the conductivity $σ$ of a given body from a measurement of the magnitude of the current density $|J|$ inside the body and voltage $g$ on its boundary. Denote by $u$ the electrical potential corresponding to the voltage $g$, then $u$ formally satisfies the equation

$$ (1.5) \quad \left\{ \begin{array}{ll} \operatorname{div}(\sigma \nabla u) = 0 & \text{in } \Omega, \\ u = g & \text{on } \partial\Omega. \end{array} \right. $$

By Ohm's law, the current density equals $J = -\sigma \nabla u$, so we may formally rewrite the above problem as the weighted 1-Laplace equation

$$ (1.6) \quad \left\{ \begin{array}{ll} -\operatorname{div}(|J| \frac{\nabla u}{|\nabla u|}) = 0 & \text{in } \Omega, \\ u = g & \text{on } \partial\Omega, \end{array} \right. $$

which is the Euler-Lagrange equation of the weighted least gradient problem with weight $k = |J|$ (for a precise justification of this passage, see [31]). The conductivity imaging problem is in turn related to a variant of the Calderón problem, in which one wants to recover the conductivity of the body $Ω$ solely from the measurement at the boundary and the knowledge of a Dirichlet-to-Neumann map associated to the weighted 1-Laplace equation (1.6) (or equivalently, to (1.5)). Let
us note that the Dirichlet-to-Neumann map for the isotropic 1-Laplace operator was studied in [22].

In [23], the authors showed existence of a solution $u$ for Problem (wLGP) as soon as $\partial\Omega$ satisfies a positivity condition on a sort of generalized mean curvature related to $k$, and $g \in C(\partial\Omega)$. Moreover, they proved that if $k \in C^{1,1}(\bar{\Omega})$ is positive and bounded away from zero and if $g$ is continuous on $\partial\Omega$, then the weighted least gradient problem (wLGP) has at most one minimizer. They also showed that the condition $k \in C^{1,1}(\bar{\Omega})$ is sharp in the sense that uniqueness may fail if $k \in C^{1,\alpha}(\bar{\Omega})$ with $\alpha < 1$. In the present paper, we show existence and uniqueness of a solution $u$ for Problem (wLGP) via an optimal transport approach. More precisely, we will refine and generalize first the result of [21] by proving that if $\Omega \subset \mathbb{R}^2$ is contractible and if $g \in BV(\partial\Omega)$, then Problem (wLGP) is completely equivalent to the weighted Beckmann problem:

$$ (\text{wBP}) \qquad \inf \left\{ \int_{\bar{\Omega}} k(x) |v| : v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2), \ \text{div } v = 0, v \cdot n = f \text{ on } \partial\Omega \right\}, $$

where $f := \partial_{\tau} g$. In particular, we have that $u$ is a solution for Problem (wLGP) if and only if the corresponding $v := R_{\frac{1}{2}} Du$ is an optimal flow for Problem (wBP) with $|v|(\partial\Omega) = 0$. Yet, one can also show that Problem (wBP) is equivalent to the following optimal transport problem with Riemannian cost (see also [12, 33]):

$$ (1.7) \qquad \min \left\{ \int_{\bar{\Omega} \times \bar{\Omega}} d_k(x, y) \, d\Lambda : \Lambda \in \mathcal{M}^+(\bar{\Omega} \times \bar{\Omega}), (\Pi_x)_{\#} \Lambda = f^+ \text{ and } (\Pi_y)_{\#} \Lambda = f^- \right\}, $$

where $d_k$ denotes the distance generated by the conformally Riemannian metric $k$. In fact, we will show that if $\Omega$ is geodesically convex, then every optimal flow $v$ of Problem (wBP) is of the form $v = k^{-1}v_{\Lambda}$, for some optimal transport plan $\Lambda$ of Problem (1.7), where

$$ <v_{\Lambda}, \xi> = \int_{\bar{\Omega} \times \bar{\Omega}} \int_{0}^{1} k(\gamma_{x,y}(t)) \xi(\gamma_{x,y}(t)) \cdot \gamma'_{x,y}(t) \, dt \, d\Lambda(x,y), \text{ for all } \xi \in C(\bar{\Omega}, \mathbb{R}^2). $$

It is not difficult to see that $|v_{\Lambda}|(\partial\Omega) = 0$ as soon as $\Omega$ is geodesically strictly convex. Thanks to this fact, we will show that Problem (wLGP) reaches a minimum provided that $g \in BV(\partial\Omega)$. On the other hand, we will also show that Problem (1.7) has a unique optimal transport plan $\Lambda$ as soon as $f^+$ or $f^-$ is atomless and $\Omega$ is geodesically strictly convex. This yields to uniqueness of the solution $u$ of Problem (wBP) provided that $g \in C(\partial\Omega)$.

Moreover, we will study the $L^p$ summability of the transport density $|v_{\Lambda}|$ (where $\Lambda$ is the unique optimal transport plan in Problem (1.7)) between two measures $f^+$ and $f^-$ which are located on the boundary $\partial\Omega$. More precisely, we will extend the $L^p$ estimates of [15] on the transport density to the case where the transport cost is given by a Riemannian distance $d_k(x, y)$, provided that $\Omega$ is geodesically uniformly convex. In fact, proving $L^p$ estimates on this transport density is significantly more difficult than in [15], since the transport rays are now Riemannian geodesics; they are not necessarily straight lines like in the Euclidean case, and the key point in the proof would be to show a geometric result concerning the directions of the transport rays. Together with the equivalence between problems (wLGP) and (wBP), this yields $W^{1,p}$ regularity on the solution $u$ of Problem (wLGP). In short, the main result in the present paper can be stated as follows: if $p \le 2$, the domain $\Omega$ is geodesically uniformly convex and $g \in W^{1,p}(\partial\Omega)$, the weighted least gradient problem (wLGP) has a unique solution $u \in W^{1,p}(\Omega)$.

The paper is organised as follows. In Section 2, we recall the basic definitions concerning anisotropic BV functions and Anzellotti pairings, as well as the state of the literature concerning the weighted least gradient problem and the equivalence between the least gradient problem and the boundary-to-boundary optimal transport problem in the Euclidean case. In Section 3, we prove the equivalence between the anisotropic least gradient problem and anisotropic Beckmann problem under minimal assumptions. Then, in Section 4, we consider Riemannian metrics and show equivalence between the weighted Beckmann problem and the optimal transport problem with Riemannian cost. Then, in Section 5, we study some consequences of the results in the previous two Sections for the structure of solutions in the least gradient problem. The main results in this paper, concerning $L^p$ estimates for the transport density, are in Section 6. Finally, in Section 7, we apply the $L^p$ estimates for the transport density to obtain regularity estimates for solutions of the weighted least gradient problem.
## 2. PRELIMINARIES

In this paper, our main goal is to study regularity of solutions to the weighted least gradient problem. To this end, we first recall in Section 2.1 the notion of anisotropic BV spaces and some of their properties. We follow the presentation in [1]. We also recall the notion of Anzellotti pairings, weak normal traces and a corresponding version of the Gauss-Green formula (see [3,7]). Then, in Section 2.2, we recall the precise statement and known results on the weighted least gradient problem. Finally, in Section 2.3, we present known results on the equivalence between the two-dimensional least gradient problem and a version of the optimal transport problem in the Euclidean case. Throughout the whole paper, $\Omega \subset \mathbb{R}^2$ denotes an open bounded set with Lipschitz boundary.

### 2.1. Anisotropic BV spaces.

**Definition 2.1.** A continuous function $\phi : \bar{\Omega} \times \mathbb{R}^2 \to [0, \infty)$ is called a metric integrand if it satisfies the following conditions:

(1) $v \mapsto \phi(x, v)$ is a norm for every $x$,

(2) $\phi$ is bounded and uniformly elliptic in $\bar{\Omega}$, i.e.

$$ \exists \lambda, \Lambda > 0 \quad \text{s.t.} \quad \forall x \in \bar{\Omega}, \quad \forall \xi \in \mathbb{R}^2, \quad \lambda|\xi| \le \phi(x, \xi) \le \Lambda|\xi|. $$

Standard examples of metric integrands which appear with relation to the least gradient problem are: the Euclidean case (i.e. $\phi(x, p) = |p|$), the $l_1$ norm $\phi(x, p) = |p_1| + |p_2|$, and the weighted least gradient problem $\phi(x, p) = k(x)|p|$ with $k$ bounded and bounded away from zero (see [17, 23, 40]). This last case will be the main focus of our attention in the second half of the paper.

The definition below is specific to the two-dimensional case.

**Definition 2.2.** Let $\phi : \bar{\Omega} \times \mathbb{R}^2 \to [0, \infty)$ be a metric integrand. We denote its rotation by $-\frac{\pi}{2}$ in the second variable by

$$ \phi^{\perp}(x, \xi) := \phi(x, R_{-\frac{\pi}{2}}\xi). $$

It is clear that $\phi^{\perp} : \bar{\Omega} \times \mathbb{R}^2 \to [0, \infty)$ is also a metric integrand.

**Definition 2.3.** Let $\phi : \bar{\Omega} \times \mathbb{R}^2 \to [0, \infty)$ be a metric integrand. Its polar function is $\phi^0 : \bar{\Omega} \times \mathbb{R}^2 \to [0, \infty)$ defined as

$$ \phi^0(x, \xi^*) = \sup \{\langle \xi^*, \xi \rangle : \xi \in \mathbb{R}^2, \phi(x, \xi) \le 1\}. $$

**Definition 2.4.** Let $\phi : \bar{\Omega} \times \mathbb{R}^2 \to [0, \infty)$ be a metric integrand. For a given function $u \in L^1(\Omega)$, we define its $\phi$-total variation in $\Omega$ by the formula:

$$ \int_{\Omega} |Du|_{\phi} = \sup \left\{ \int_{\Omega} u \operatorname{div}(\mathbf{z}) \, dx : \phi^0(x, \mathbf{z}(x)) \le 1 \quad \text{a.e.,} \quad \mathbf{z} \in C_c^1(\Omega) \right\}. $$

Another popular notation for the $\phi$-total variation is $\int_{\Omega} \phi(x, Du)$. In fact, one can show that $\int_{\Omega} \phi(x, Du) = \int_{\Omega} \phi(x, \frac{dDu}{d|Du|})d|Du|$, where $\frac{dDu}{d|Du|}$ denotes the Radon-Nikodym derivative. We will say that $u \in BV_{\phi}(\Omega)$ if its $\phi$-total variation is finite in $\bar{\Omega}$; furthermore, let us define the $\phi$-perimeter of a set $E$ as

$$ P_{\phi}(E, \Omega) = \int_{\Omega} |D\chi_{E}|_{\phi}. $$

If $P_{\phi}(E, \Omega) < \infty$, we say that $E$ is a set of bounded $\phi$-perimeter in $\Omega$.

**Remark 2.5.** Let $\phi : \bar{\Omega} \times \mathbb{R}^2 \to [0, \infty)$ be a metric integrand. Since $\phi$ is bounded and uniformly elliptic, we have

$$ \lambda \int_{\Omega} |Du| \le \int_{\Omega} |Du|_{\phi} \le \Lambda \int_{\Omega} |Du|. $$

In particular, $BV_{\phi}(\Omega) = BV(\Omega)$ as sets. They are equipped with different (but equivalent) norms.

Now, we recall the definition and basic properties of Anzellotti pairings introduced in [3]. The construction below is usually performed under more general assumptions, but for simplicity we restrict ourselves to the setting we will use in the paper (in Section 5). For $p \ge 1$, denote

$$ X_p(\Omega) = \left\{ \mathbf{z} \in L^{\infty}(\Omega, \mathbb{R}^2) : \operatorname{div}(\mathbf{z}) \in L^p(\Omega) \right\}. $$
Given $\mathbf{z} \in X_p(\Omega)$ and $w \in BV(\Omega) \cap L^{p'}(\Omega)$, we define the functional $(\mathbf{z}, Dw) : C_c^\infty(\Omega) \to \mathbb{R}$ by the formula

$$ \langle (\mathbf{z}, Dw), \varphi \rangle = - \int_{\Omega} w \varphi \operatorname{div}(\mathbf{z}) dx - \int_{\Omega} w \mathbf{z} \cdot \nabla \varphi dx. $$

The distribution $(\mathbf{z}, Dw)$ turns out to be a Radon measure on $\Omega$. It is a generalization of the pointwise product $\mathbf{z} \cdot \nabla w$; namely, for $w \in W^{1,1}(\Omega) \cap L^{p'}(\Omega)$, we have

$$ \int_{\Omega} (\mathbf{z}, Dw) \varphi = \int_{\Omega} \mathbf{z} \cdot \nabla w \varphi dx, \quad \forall \varphi \in C_c^\infty(\Omega). $$

The following Proposition summarizes the most important properties of the pairing $(\mathbf{z}, Du)$.

**Proposition 2.6.** Suppose that $\phi$ is a metric integrand. Let $\mathbf{z} \in X_p(\Omega)$ and $u \in BV(\Omega) \cap L^{p'}(\Omega)$. Then, for any Borel set $B \subset \Omega$, we have

$$ \left| \int_B (\mathbf{z}, Du) \right| \le \|\phi^0(x, \mathbf{z}(x))\|_{L^\infty(\Omega)} \int_B |Du|_\phi. $$

In particular, $(\mathbf{z}, Du) \ll |Du|$ as measures in $\Omega$. Moreover, there exists a function $[\mathbf{z}, n] \in L^\infty(\partial\Omega)$ such that $\| [\mathbf{z}, n] \|_{L^\infty(\partial\Omega)} \le \| \mathbf{z} \|_{L^\infty(\Omega, \mathbb{R}^2)}$ and the following Gauss-Green formula holds:

$$ \int_{\Omega} u \operatorname{div}(\mathbf{z}) dx + \int_{\Omega} (\mathbf{z}, Du) = \int_{\partial\Omega} [\mathbf{z}, n] u d\mathcal{H}^1. $$

The function $[\mathbf{z}, n]$ has the interpretation of the normal trace of the vector field $\mathbf{z}$ at the boundary and it coincides with the classical normal trace if $\mathbf{z}$ is sufficiently smooth.

## 2.2. Weighted least gradient problem.
In this paper, we are interested in the following problem

$$ (2.1) \qquad \inf \left\{ \int_{\Omega} k(x) |Du| : u \in BV(\Omega), u|_{\partial\Omega} = g \right\}, $$

where $k$ is sufficiently smooth. In the isotropic case, a standard assumption for existence of solutions is strict convexity of $\Omega$. In the anisotropic case, the standard assumption is called the barrier condition and is a local property at every point $x_0 \in \partial\Omega$. In the definition below, denote $\phi(x, p) = k(x)|p|$. Some of the results below are valid for more general $\phi$, but for the sake of presentation, we restrict ourselves to the two-dimensional weighted case.

**Definition 2.7 (Barrier condition).** We say that $\Omega$ satisfies the barrier condition if the following condition holds: for every $x_0 \in \partial\Omega$, there exists $r_0 > 0$ such that for all $r < r_0$, if $V$ is a minimizer of

$$ (2.2) \qquad \inf \left\{ P_{\phi}(W, \mathbb{R}^2) : W \subset \Omega, (\Omega \setminus W) \subset B(x_0, r) \right\}, $$

then

$$ \partial V \cap \partial \Omega \cap B(x_0, r) = \emptyset. $$

Under this assumption, existence of solutions was proved for continuous boundary data by Jerrard, Moradifam and Nachman in [23] (actually, the result holds also when the boundary datum is continuous a.e. with respect to the codimension one Hausdorff measure on $\partial\Omega$, see [32] or [19]). However, in order to use optimal transport techniques, we will later need a slightly stronger assumption of strict geodesic convexity of $\Omega$. In the Euclidean case, these two concepts are closely related: in two dimensions, the barrier condition is equivalent to strict convexity of $\Omega$; in higher dimensions, the barrier condition is something between convexity and strict convexity of $\Omega$, see [19]. We will comment on the relationship between these two concepts at the end of Section 4.

The literature regarding regularity of solutions to the weighted least gradient problem is very thin. In the paper [23], the authors have shown that for sufficiently smooth $k$, in low dimensions continuity of boundary data implies continuity of solutions. This result was extended to any dimension by Zuniga in [40] under the assumption $k \in C^2(\bar{\Omega})$. However, if $k$ is not regular enough, then continuity of solutions may break down: several examples for Lipschitz weights appear in [25] and an example for $C^{1,\alpha}$ weights with $\alpha < 1$ appears in [23].

Finally, let us recall another formulation of the least gradient problem, first introduced by Mazón, Rossi and Segura de León in the isotropic case (see [28]).
**Definition 2.8.** We say that $u \in BV(\Omega)$ is a function of $\phi$-least gradient if for all $v \in BV(\Omega)$ such that $v|_{\partial\Omega} = 0$, we have

$$ \int_{\Omega} |Du|_{\phi} \le \int_{\Omega} |D(u+v)|_{\phi}. $$

Using this definition, we may reformulate the weighted least gradient problem as follows: the goal is to find a $\phi$-least gradient function with prescribed trace. The following characterisation in terms of Anzellotti pairings was proved in [28] in the isotropic case; in the form below it was proved by Mazón in [27].

**Theorem 2.9.** Suppose that $u \in BV(\Omega)$ is a function of $\phi$-least gradient. Then, there exists a vector field $z \in L^{\infty}(\Omega, \mathbb{R}^2)$ such that $\|\phi^0(x, z(x))\|_{\infty} \le 1$, $\mathrm{div}(z) = 0$ in the distributional sense and

$$ (z, Du) = |Du|_{\phi} \quad \text{as measures in } \Omega. $$

We also recall the celebrated Bombieri-de Giorgi-Giusti theorem. It was first proved in the Euclidean case in [5] and in the form below it was proved by Mazón in [27].

**Theorem 2.10.** Suppose that $u \in BV(\Omega)$ is a $\phi$-least gradient function. Then, for all $t \in \mathbb{R}$, the function $\chi_{\{u \ge t\}}$ is also a function of $\phi$-least gradient.

**2.3. Equivalence in the Euclidean case.** In this Section, we present the relationship between the least gradient problem, the Beckmann problem and the classical Monge-Kantorovich problem in the Euclidean case, together with some properties of solutions to these problems. Recall that throughout the paper $\Omega \subset \mathbb{R}^2$ denotes an open bounded set with Lipschitz boundary. In this subsection, we assume additionally that $\Omega$ is convex. We start with the relationship between the least gradient problem

$$ (2.3) \qquad \inf \left\{ \int_{\Omega} |Du| : u \in BV(\Omega), u|_{\partial\Omega} = g \right\}, $$

and the Beckmann problem

$$ (2.4) \qquad \min \left\{ \int_{\bar{\Omega}} |v| : v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2), \mathrm{div} \, v = f \right\}. $$

Here, $g \in BV(\partial\Omega)$, and the boundary condition in (2.3) is understood in the sense of traces. The divergence condition in (2.4) is understood in the distributional sense; we have $\mathrm{div}(v) = 0$ in $\Omega$ and $[v, n] = f$ on $\partial\Omega$, where $[v, n]$ denotes the weak normal trace of a vector field whose divergence is integrable (see Section 2.1). Suppose that $f \in \mathcal{M}(\partial\Omega)$ satisfies a mass balance condition $\int_{\partial\Omega} df = 0$. Then, Problem (2.4) admits a solution (see [36]). Moreover, Problem (2.3) reaches a minimum as soon as $\Omega$ is strictly convex (see [17]). Notice that such $f$ and $g$ are in a one-to-one correspondence (up to an additive constant in $g$) via the relation $f = \partial_{\tau}g$. It was observed in [21] that Problems (2.3) and (2.4) are closely related. Namely, if we take an admissible function $u \in BV(\Omega)$ in (2.3), then $v = R_{\frac{\pi}{2}}Du$ is admissible in (2.4). Indeed, in dimension two, a rotation of a gradient by $\frac{\pi}{2}$ is a divergence-free field in $\Omega$ and it interchanges the normal and tangent components at the boundary. In the other direction, given a vector field $v \in L^1(\Omega, \mathbb{R}^2)$ admissible in (2.4), we can recover $u \in W^{1,1}(\Omega)$ admissible in (2.3). In [15], the authors noticed that this result may be improved: if $v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2)$ is such that $|v|(\partial\Omega) = 0$, then there exists $u \in BV(\Omega)$ such that $v = R_{\frac{\pi}{2}}Du$ and $[v, n] = \partial_{\tau}(Tu)$. In particular, notice that $|v| = |Du|$ as measures on $\bar{\Omega}$; hence, the infimal values in (2.3) and (2.4) coincide. We will see later that a solution $v$ for Problem (2.4) gives zero mass to the boundary as soon as $\Omega$ is strictly convex (this implies existence of a solution for Problem (2.3); see [15]). Let us sum up these considerations as follows.

**Theorem 2.11.** Suppose that $\Omega$ is convex. Then, the problems (2.3) and (2.4) are equivalent in the following sense:

(1) Their infimal values coincide, i.e. $\inf(2.3) = \min(2.4)$.

(2) Given a solution $u \in BV(\Omega)$ of (2.3), we can construct a solution $v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2)$ of (2.4); moreover, $v = R_{\frac{\pi}{2}}Du$.

(3) Given a solution $v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2)$ of (2.4) with $|v|(\partial\Omega) = 0$, we can construct a solution $u \in BV(\Omega)$ of (2.3); moreover, $v = R_{\frac{\pi}{2}}Du$.
Now, we turn to the equivalence between the Beckmann problem (2.4) and the following Monge-Kantorovich problem:

$$ (2.5) \quad \min \left\{ \int_{\overline{\Omega} \times \overline{\Omega}} |x-y| \, d\Lambda : \Lambda \in \mathcal{M}^+(\overline{\Omega} \times \overline{\Omega}), (\Pi_x)_\# \Lambda = f^+ \text{ and } (\Pi_y)_\# \Lambda = f^- \right\}. $$

This equivalence is standard in the optimal transport theory (see for instance [36, Chapter 4]), but we will describe it shortly for completeness. Here, $\Omega \subset \mathbb{R}^2$ is a bounded open convex set, $f \in \mathcal{M}(\overline{\Omega})$ and $f = f^+ - f^-$ is its decomposition into a positive and negative part. In order for the problem to be well defined, we assume that the mass balance condition $f^+(\overline{\Omega}) = f^-(\overline{\Omega})$ holds. In particular, this setting covers our case, when the measures are concentrated on $\partial\Omega$.

First, let us recall a few standard results in the optimal transport theory. The fact that

$$ (2.6) \quad \sup \left\{ \int_{\overline{\Omega}} \psi \, d(f^+ - f^-) : \psi \in \text{Lip}_1(\overline{\Omega}) \right\} $$

is the dual problem to (2.5) is well-known, see for instance [36,39]. As a corollary of its proof, we get that there exist solutions to both problems, and that any optimal transport plan $\Lambda$ and any Kantorovich potential $\psi$ satisfy the following equality:

$$ \int_{\overline{\Omega} \times \overline{\Omega}} (|x-y| - (\psi(x) - \psi(y))) \, d\Lambda(x, y) = 0, $$

which implies that

$$ \psi(x) - \psi(y) = |x - y| \quad \text{on spt}(\Lambda). $$

If $\psi$ is a Kantorovich potential, we call any maximal segment $[x, y]$ satisfying $\psi(x) - \psi(y) = |x - y|$ a transport ray. The Kantorovich potential is in general not unique, but frame of transport rays does not depend on the choice of $\psi$ (at least on the support of $\Lambda$). Moreover, an optimal transport plan $\Lambda$ has to move the mass along the transport rays.

Now, we turn our attention to the equivalence between the Kantorovich problem (2.5) and the Beckmann problem (2.4). First, let us see that for any $v \in \mathcal{M}(\overline{\Omega}, \mathbb{R}^d)$ admissible in the Beckmann problem (2.4), we have that for any $C^1$ function $\phi$ with $|\nabla\phi| \le 1$,

$$ |v|(\overline{\Omega}) = \int_{\overline{\Omega}} 1 \, d|v| \ge \int_{\overline{\Omega}} (-\nabla\phi) \cdot dv = \int_{\overline{\Omega}} \phi \, df. $$

Hence, $\min(2.4) \ge \max(2.6) = \min(2.5)$ (the supremum in (2.6) is taken for Lipschitz functions, but we may approximate them uniformly by $C^1$ functions). On the other hand, given an optimal transport plan $\Lambda$, we may construct a vector measure $v_\Lambda \in \mathcal{M}(\overline{\Omega}, \mathbb{R}^d)$ defined by the formula

$$ \langle v_\Lambda, \xi \rangle := \int_{\overline{\Omega} \times \overline{\Omega}} \int_0^1 \xi(\omega_{x,y}(t)) \cdot \omega'_{x,y}(t) \, dt \, d\Lambda(x, y), $$

for all $\xi \in C(\overline{\Omega}, \mathbb{R}^d)$. Here, $\omega_{x,y}(t) = (1-t)x + ty$ is the constant-speed parametrisation of $[x, y]$. By taking $\xi = \nabla\phi$, it is immediate that $v_\Lambda$ satisfies the divergence constraint. The total mass of $v_\Lambda$ will be estimated using the transport density $\sigma_\Lambda \in \mathcal{M}^+(\overline{\Omega})$, which is defined by the formula

$$ \langle \sigma_\Lambda, \phi \rangle := \int_{\overline{\Omega} \times \overline{\Omega}} \int_0^1 \phi(\omega_{x,y}(t)) |\omega'_{x,y}(t)| \, dt \, d\Lambda(x, y), $$

for all $\phi \in C(\overline{\Omega})$. The vector measure $v_\Lambda$ and the scalar measure $\sigma_\Lambda$ are related in the following way: if $\psi$ is a Kantorovich potential, we have

$$ \omega'_{x,y}(t) = y - x = -|x - y| \frac{x - y}{|x - y|} = -|x - y| \nabla\psi(\omega_{x,y}(t)), $$

for all $t \in (0, 1)$ and $(x, y) \in \text{spt}(\Lambda)$. Thus, $\langle v_\Lambda, \xi \rangle = \langle \sigma_\Lambda, -\xi \cdot \nabla\psi \rangle$, for all $\xi \in C(\overline{\Omega}, \mathbb{R}^d)$, so we have

$$ v_\Lambda = -\nabla\psi \cdot \sigma_\Lambda. $$

Hence, $v_\Lambda$ is absolutely continuous with respect to $\sigma_\Lambda$ and $|v_\Lambda| = \sigma_\Lambda$. Then, one has

$$ \min (2.5) = \int_{\overline{\Omega} \times \overline{\Omega}} |x - y| \, d\Lambda = \int_{\overline{\Omega} \times \overline{\Omega}} \int_0^1 |\omega'_{x,y}(t)| \, dt \, d\Lambda(x, y) = \sigma_\Lambda(\overline{\Omega}) = |w_\Lambda|(\overline{\Omega}) \ge \min (2.4). $$

This implies that $\min(2.5) = \min(2.4)$ and that from an optimal transport plan $\Lambda$ we can construct a solution to the Beckmann problem (2.4). On the other hand, it can be shown that every solution to (2.4) is of the form $v = v_\Lambda$ for some optimal transport plan $\Lambda$, see [36, Theorem 4.13]. We summarize the above discussion in the following Theorem.
**Theorem 2.12.** Suppose that $\Omega$ is convex. Then, the problems (2.5) and (2.4) admit solutions, and are equivalent in the following sense:

(1) Their minimal values coincide, i.e. $\min(2.5) = \min(2.4)$.

(2) Given an optimal transport plan $\Lambda \in \mathcal{M}^+(\bar{\Omega} \times \bar{\Omega})$ in (2.5), we can construct a solution $v_\Lambda \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^d)$ to (2.4).

(3) Given a solution $v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2)$ to (2.4), we can construct an optimal transport plan $\Lambda \in \mathcal{M}^+(\bar{\Omega} \times \bar{\Omega})$ in (2.5) such that $v = v_\Lambda$.

The equivalence between the two-dimensional least gradient problem (2.3) and the Monge-Kantorovich problem (2.5) comes from combining Theorems 2.11 and 2.12. Given a solution to (2.3), we may construct an optimal transport plan $\Lambda$ for (2.5) with $f^\pm = (\partial_\tau g)^\pm$; in the other direction, since $|v_\Lambda| = \sigma_\Lambda$, we may recover a solution to (2.3) from an optimal transport plan $\Lambda$ as soon as the transport density $\sigma_\Lambda$ gives no mass to the boundary, i.e. $\sigma_\Lambda(\partial\Omega) = 0$. An important special case is when $\Omega$ is strictly convex. Using the equivalent formula for transport density, i.e. for every Borel set $A \subset \bar{\Omega}$,

$$ \sigma_\Lambda(A) = \int_{\bar{\Omega} \times \bar{\Omega}} \mathcal{H}^1([x, y] \cap A) d\Lambda(x, y), $$

we have that if $\Omega$ is strictly convex, then for any optimal transport plan $\Lambda$ we have that $\sigma_\Lambda(\partial\Omega) = 0$, and the correspondence between problems is one-to-one. This link between the two problems was exploited for the first time in [15]. For a strictly convex domain $\Omega$, the authors studied the boundary-to-boundary optimal transport problem. If at least one of the measures $f^\pm$ is atomless, then the optimal transport plan is unique and induced by a map, and proved several variants of regularity estimates on the transport density. The main application of these results so far is the $W^{1,p}$ regularity of solutions to the least gradient problem: the estimates on the transport density imply that if $g \in W^{1,p}(\partial\Omega)$ for $p \in [1, 2]$, then the unique solution $u$ to (2.3) lies in $W^{1,p}(\Omega)$, provided that $\Omega$ is uniformly convex.

### 3. EQUIVALENCE BETWEEN THE ANISOTROPIC LEAST GRADIENT PROBLEM AND THE ANISOTROPIC BECKMANN PROBLEM

In this Section, we will prove the equivalence between the anisotropic least gradient problem

$$ \inf \left\{ \int_{\Omega} \phi(x, Du) : u \in BV(\Omega), u|_{\partial\Omega} = g \right\}, $$

and the anisotropic Beckmann problem

$$ \min \left\{ \int_{\bar{\Omega}} \phi^{\perp}(x, v) : v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2), \operatorname{div} v = f \right\}. $$

Recall that throughout the paper $\Omega \subset \mathbb{R}^2$ denotes an open bounded set with Lipschitz boundary. In this Section, we will additionally require that $\Omega$ is contractible (or equivalently, since we study the planar case, that $\partial\Omega$ is connected). This is unlike the results in [21] for the Euclidean case, which require convexity of the domain. We will generalize the results from [21] in two ways: by allowing anisotropic settings and a wider class of domains. Relaxing the convexity assumption on $\Omega$ is particularly relevant in the anisotropic case; throughout most of the paper, we will instead require that the domain $\Omega$ is geodesically convex. Hence, we need to prove the equivalence for a more general class of domains.

**Proposition 3.1.** Suppose that $\Omega$ is contractible. Let $v \in L^1(\Omega, \mathbb{R}^2)$ be such that $\operatorname{div}(v) = 0$ in the sense of distributions. Then, there exists $u \in W^{1,1}(\Omega)$ such that $v = R_{\frac{\pi}{2}}\nabla u$. In particular, for any metric integrand $\phi$, we have

$$ \int_{\Omega} \phi^{\perp}(x, v) dx = \int_{\Omega} \phi(x, \nabla u) dx. $$

Moreover, if $v \cdot n|_{\partial\Omega} = f$ and $Tu = g$, then we have $f = \partial_\tau g$.

**Proof.** Step 1. Definition of an approximating sequence $u^\epsilon$. Fix any $x_0 \in \Omega$. Denote $v = (v_1, v_2)$. Let $\varphi_\epsilon$ (with $\epsilon > 0$) be a sequence of mollifiers and set $v^\epsilon$ to be the mollification of $v$, i.e. $v^\epsilon = v * \varphi_\epsilon := (v_1^\epsilon, v_2^\epsilon)$. Since $\operatorname{div}(v) = 0$, we also have $\operatorname{div}(v^\epsilon) = 0$. We define a sequence of mollified differential forms $\omega^\epsilon = v_2^\epsilon dx_1 - v_1^\epsilon dx_2$; since $\operatorname{div}(v^\epsilon) = 0$, we also have $d\omega^\epsilon = 0$. The
domain is contractible, so the differential form $\omega^\epsilon$ uniquely (up to an additive constant depending on the choice of $x_0$) defines a function $u^\epsilon \in C^\infty(\bar{\Omega})$ via

$$u^{\epsilon}(x) = \int_{\gamma_x} \omega^{\epsilon},$$

with $\gamma_x$ is any smooth oriented curve from $x_0$ to $x$. In particular, we have $\nabla u^\epsilon = R_{-\frac{\pi}{2}} v^\epsilon$.

**Step 2. Convergence in $L^1(\Omega)$.** We claim that a modified version of the sequence $u^\epsilon$ converges in $L^1(\Omega)$. Denote by $(u^\epsilon)_\Omega$ the mean value of $u^\epsilon$ on $\Omega$. Set $\tilde{u}^\epsilon = u^\epsilon - (u^\epsilon)_\Omega$ and notice that $\nabla \tilde{u}^\epsilon = \nabla u^\epsilon = R_{-\frac{\pi}{2}} v^\epsilon$. Using the Poincaré inequality, we get the following estimate:

$$(3.1) \qquad \int_{\Omega} |\tilde{u}^{\epsilon}| dx = \int_{\Omega} |u^{\epsilon} - (u^{\epsilon})_{\Omega}| dx \le C(\Omega) \int_{\Omega} |\nabla u^{\epsilon}| dx = C(\Omega) \int_{\Omega} |v^{\epsilon}| dx.$$

Yet, the sequence $v^\epsilon$ is bounded in $L^1(\Omega, \mathbb{R}^2)$ with $\|v^\epsilon\|_{L^1} \le \|v\|_{L^1}$. Hence, we have

$$\|\tilde{u}^{\epsilon}\|_{W^{1,1}(\Omega)} \le (C(\Omega) + 1)\|v^{\epsilon}\|_{L^1(\Omega)} \le M.$$

Then, there exists a subsequence $\tilde{u}^{\epsilon_n} \to u$ in $L^1(\Omega)$.

**Step 3. Convergence in $W^{1,1}(\Omega)$**. We recall that $\nabla \tilde{u}^\epsilon = R_{-\frac{\pi}{2}} v^\epsilon$. Since $v^\epsilon$ is a Cauchy sequence in $L^1(\Omega, \mathbb{R}^2)$, $\nabla \tilde{u}^{\epsilon_n}$ is also a Cauchy sequence in $L^1(\Omega, \mathbb{R}^2)$, because

$$\|\nabla \tilde{u}^{\epsilon_n} - \nabla \tilde{u}^{\epsilon_m}\|_{L^1(\Omega, \mathbb{R}^2)} = \|R_{-\frac{\pi}{2}}(v^{\epsilon_n} - v^{\epsilon_m})\|_{L^1(\Omega, \mathbb{R}^2)} = \|v^{\epsilon_n} - v^{\epsilon_m}\|_{L^1(\Omega, \mathbb{R}^2)}.$$

Hence, $\tilde{u}^{\epsilon_n}$ is a Cauchy sequence in $W^{1,1}(\Omega)$, so it converges to $u$ in $W^{1,1}(\Omega)$. In particular, we have $\nabla u = R_{-\frac{\pi}{2}} v$ and,

$$\int_{\Omega} \phi(x, \nabla u) \, dx = \int_{\Omega} \phi^{\perp}(x, v) \, dx.$$

**Step 4. The boundary condition.** Finally, we study the relationship between the trace of $f$ and the normal trace of $g$; this part is similar to part of the proof in [21], but we write it here for completeness. If $\phi \in \text{Lip}(\mathbb{R}^2)$, and $\phi$ is its restriction on $\partial\Omega$, then

$$\begin{aligned} \langle f, \phi \rangle &= \langle v \cdot n|_{\partial\Omega}, \phi \rangle = \int_{\Omega} v \cdot \nabla \phi \, dx = \lim_{\epsilon_n \to 0} \int_{\Omega} v^\epsilon \cdot \nabla \phi \, dx = \lim_{\epsilon_n \to 0} \int_{\Omega} R_{\frac{\pi}{2}} \nabla \tilde{u}^{\epsilon_n} \cdot \nabla \phi \, dx \\ &= \lim_{\epsilon_n \to 0} \int_{\partial\Omega} R_{\frac{\pi}{2}} \nabla \tilde{u}^{\epsilon_n} \cdot n \, \phi \, d\mathcal{H}^1. \end{aligned}$$

Due to the smoothness of $\tilde{u}^{\epsilon_n}$, we have $R_{\frac{\pi}{2}} \nabla \tilde{u}^{\epsilon_n} \cdot n = \partial_{\tau} \tilde{u}^{\epsilon_n}$. Hence,

$$\langle f, \phi \rangle = \lim_{\epsilon_n \to 0} \int_{\partial \Omega} R_{\frac{\pi}{2}} \nabla \tilde{u}^{\epsilon_n} \cdot n \, \phi \, d\mathcal{H}^1 = \lim_{\epsilon_n \to 0} \int_{\partial \Omega} \partial_{\tau} \tilde{u}^{\epsilon_n} \, \phi \, d\mathcal{H}^1 = - \lim_{\epsilon_n \to 0} \int_{\partial \Omega} T \tilde{u}^{\epsilon_n} \partial_{\tau} \phi \, d\mathcal{H}^1.$$

Since $\tilde{u}^{\epsilon_n}$ converges to $u$ in $W^{1,1}(\Omega)$, then the traces converge in $L^1(\partial\Omega)$, i.e., $T \tilde{u}^{\epsilon_n} \to Tu$. Thus,

$$\langle f, \phi \rangle = - \int_{\partial \Omega} T u \partial_{\tau} \phi \, d\mathcal{H}^1 = - \int_{\partial \Omega} g \partial_{\tau} \phi \, d\mathcal{H}^1 = \langle \partial_{\tau} g, \phi \rangle. \quad \square$$

In a similar way, we may enlarge the class of admissible vector fields $v$. We will require that $v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2)$ is such that $|v|(\partial\Omega) = 0$. The outline of this proof will be very similar to the proof of Proposition 3.1.

**Proposition 3.2.** Suppose that $\Omega$ is contractible. Let $v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2)$ be such that $\text{div}(v) = 0$ in the sense of distributions and $|v|(\partial\Omega) = 0$. Then, there exists $u \in \text{BV}(\Omega)$ such that $v = R_{\frac{\pi}{2}} Du$. In particular, for any metric integrand $\phi$, we have

$$\int_{\Omega} \phi^{\perp}(x, v) \, dx = \int_{\Omega} \phi(x, Du) \, dx.$$

In addition, if $v \cdot n|_{\partial\Omega} = f$ and $Tu = g$, then we have $f = \partial_{\tau} g$.

**Proof. Step 1. Definition of an approximating sequence $u^\epsilon$**. We proceed exactly as in the proof of Proposition 3.1; we also keep the same notation.

**Step 2. Convergence in $L^1(\Omega)$**. The estimate (3.1) is proved on the level of the approximation $\tilde{u}^\epsilon$ and remains the same when $v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2)$. When $v$ is a measure, we also have $v^\epsilon$ is bounded in $L^1(\Omega, \mathbb{R}^2)$ with $\|v^\epsilon\|_{L^1(\Omega, \mathbb{R}^2)} \le \|v\|_{\mathcal{M}(\bar{\Omega}, \mathbb{R}^2)}$. Hence

$$\|\tilde{u}^{\epsilon}\|_{W^{1,1}(\Omega)} \le M.$$
Then, there exists a subsequence $\tilde{u}^{\varepsilon_n}$ such that $\tilde{u}^{\varepsilon_n} \to u$ weakly* in $BV(\Omega)$, i.e. there is a function $u \in BV(\Omega)$ such that $u^{\varepsilon_n} \to u$ in $L^1(\Omega)$ and $\nabla u^{\varepsilon_n} d\mathcal{L}^2 \to Du$ weakly* as measures. In particular, we have $v = R_{\frac{\pi}{2}}Du$.

**Step 3.** **Convergence of the anisotropic total variations.** Now, recall that $|v|(\partial\Omega) = 0$; then, also $|\phi^{\perp}(x, v)|(\partial\Omega) = 0$. Hence, one has

$$ \int_{\Omega} \phi(x, Du) = \int_{\Omega} \phi^{\perp}(x, v). $$

Moreover, we have

$$ \int_{\Omega} |v| = \lim_{\varepsilon_n \to 0} \int_{\Omega} |v^{\varepsilon_n}| dx = \lim_{\varepsilon_n \to 0} \int_{\Omega} |\nabla \tilde{u}^{\varepsilon_n}| dx \geq \int_{\Omega} |Du| = \int_{\Omega} |v|, $$

and so, the inequality is in fact an equality and the subsequence $\tilde{u}^{\varepsilon_n} \to u$ strictly in $BV(\Omega)$.

**Step 4. The boundary condition.** This step is handled analogously to the corresponding part in the proof of Proposition 3.1. The only difference is that the convergence $Tu^{\varepsilon_n} \to Tu$ follows from strict convergence of the sequence $\tilde{u}^{\varepsilon_n}$ to $u$ in place of convergence in norm. $\square$

Finally, as a consequence of Step 4 in the proof of Proposition 3.1, we get a converse result to Proposition 3.1. Actually, the next result is valid without the contractibility assumption on $\Omega$, and a proof can be found in [14, Proposition 3.1].

**Proposition 3.3.** Let $u \in BV(\Omega)$ with trace $Tu = g$. Then, $v = R_{\frac{\pi}{2}}\nabla u$ is a vector-valued measure such that $\text{div}(v) = f$, where $f = \partial_{\tau}g$. In particular, it is an admissible function in (2.4).

Now, we state the main Theorem.

**Theorem 3.4.** Suppose that $\Omega$ is contractible. Then, the problems (2.3) and (2.4) are equivalent in the following sense:

(1) Their infimal values coincide, i.e. $\inf(2.3) = \min(2.4)$.

(2) Given a solution $u \in BV(\Omega)$ of (2.3), we can construct a solution $v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2)$ of (2.4); moreover, $v = R_{\frac{\pi}{2}}Du$.

(3) Given a solution $v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2)$ of (2.4) with $|v|(\partial\Omega) = 0$, we can construct a solution $u \in BV(\Omega)$ of (2.3); moreover, $v = R_{\frac{\pi}{2}}Du$.

**Proof.** (1) Suppose that $(u_n)_n \subset BV(\Omega)$ is a minimizing sequence for the problem (2.3). Then, taking $v_n = R_{\frac{\pi}{2}}Du_n$, we have by Proposition 3.3 that the functions $v_n$ are admissible in Problem (2.4) and we have

$$ \inf(2.3) \leftarrow \int_{\Omega} \phi(x, Du_n) = \int_{\Omega} \phi^{\perp}(x, v_n) \geq \min(2.4). $$

Similarly, suppose that $v_n \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2)$ is a minimizing sequence for the problem (2.4). Mollifying this sequence if necessary, we may take another minimizing sequence $\tilde{v}_n \in L^1(\Omega, \mathbb{R}^2)$ for this problem. By Proposition 3.1, there exist functions $u_n \in W^{1,1}(\Omega)$ admissible in Problem (2.3) such that $\tilde{v}_n = R_{\frac{\pi}{2}}\nabla u_n$. Hence,

$$ \min(2.4) \leftarrow \int_{\Omega} \phi^{\perp}(x, \tilde{v}_n) = \int_{\Omega} \phi(x, \nabla u_n) \geq \inf(2.3). $$

We obtained that $\inf(2.3) = \min(2.4)$.

(2) Let $u \in BV(\Omega)$ be a minimizer of Problem (2.3). Let $v = R_{\frac{\pi}{2}}\nabla u$; by Proposition 3.3, it is an admissible vector field in Problem (2.4). We have

$$ \inf(2.3) = \int_{\Omega} \phi(x, Du) = \int_{\Omega} \phi^{\perp}(x, v) \geq \min(2.4). $$

Hence, $v$ is a minimizer for Problem (2.4).

(3) Let $v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2)$ be a minimizer of Problem (2.4) such that $|v|(\partial\Omega) = 0$. By Proposition 3.2, there exists $u \in BV(\Omega)$ admissible in Problem (2.3) and $v = R_{\frac{\pi}{2}}\nabla u$. Hence, one has

$$ \min(2.4) = \int_{\Omega} \phi^{\perp}(x, v) = \int_{\Omega} \phi(x, Du) \geq \inf(2.3). $$

Consequently, this implies that the function $u$ is in fact a solution for the problem (2.3). $\square$
This concludes the equivalence between the anisotropic least gradient problem and the anisotropic Beckmann problem for any contractible domain. Let us note that if the domain is not contractible, then the results in this section are no longer true even in the Euclidean case (see for instance [14], where the authors study the least gradient problem on annulus).

### 4. OPTIMAL TRANSPORT PROBLEM WITH RIEMANNIAN COST

Let $k$ be a smooth (say $C^{1,1}$) positive function on $\mathbb{R}^2$. We denote by $d_k$ the Riemannian metric associated with $k$:

$$d_k(x, y) := \min \left\{ \int_0^1 k(\gamma(t)) |\gamma'(t)| dt : \gamma \in \text{Lip}([0, 1], \mathbb{R}^2), \gamma(0) = x \text{ and } \gamma(1) = y \right\}.$$

Let $\Omega$ be an open and geodesically convex domain in $\mathbb{R}^2$ (i.e., for any two points $x$ and $y$ in $\bar{\Omega}$, there is a unique geodesic $\gamma$ contained within $\bar{\Omega}$ that joins $x$ and $y$). Let $f^+$ and $f^-$ be two nonnegative Borel measures on $\bar{\Omega}$ such that $f^+(\bar{\Omega}) = f^-(\bar{\Omega})$. Then, we consider the Kantorovich problem:

$$(4.1) \quad \min \left\{ \int_{\bar{\Omega} \times \bar{\Omega}} d_k(x, y) d\Lambda(x, y) : \Lambda \in \mathcal{M}^+(\bar{\Omega} \times \bar{\Omega}), (\Pi_x)_\# \Lambda = f^+ \text{ and } (\Pi_y)_\# \Lambda = f^- \right\}.$$

In fact, Problem (4.1) is the relaxed version of the Monge problem with Riemannian cost:

$$(4.2) \quad \min \left\{ \int_{\bar{\Omega}} d_k(x, T(x)) df^+(x) : T_\# f^+ = f^- \right\}.$$

The authors of [16] prove existence of an optimal transport map $T$ or equivalently, an optimal transport plan $\Lambda$ which is concentrated on a map, under the assumption that $f^+ \in L^1(\Omega)$. From [36, 39], the problem (4.1) admits a dual formulation

$$(4.3) \quad \sup \left\{ \int_{\bar{\Omega}} \psi d(f^+ - f^-) : |\psi(x) - \psi(y)| \le d_k(x, y), \forall x, y \in \bar{\Omega} \right\}.$$

Notice that a function $\psi$ is 1-Lipschitz with respect to the geodesic distance $d_k$ if and only if $|\nabla\psi(x)| \le k(x)$ for almost every $x$. So, the idea in [16] was the following: to obtain an optimal transport map $T$ in the problem (4.2), it is sufficient to start from a Kantorovich potential $\psi$ (i.e., a maximizer of the dual problem (4.3)) and then, construct a transport map $T$ in such a way that $\psi$ and $T$ satisfy together the following:

$$\psi(x) - \psi(T(x)) = d_k(x, T(x)), \text{ for } f^+ \text{ a.e. } x \in \bar{\Omega}.$$

Anyway, we see that if $\Lambda$ is an optimal transport plan in Problem (4.1) and if $\psi$ is a Kantorovich potential in Problem (4.3), then from the duality $\min(4.1) = \sup(4.3)$, we infer that

$$\psi(x) - \psi(y) = d_k(x, y), \text{ for all } (x, y) \in \text{spt}(\Lambda).$$

We call any maximal geodesic $\gamma_{x,y}$ between $x$ and $y$ that satisfies the equality $\psi(x)-\psi(y) = d_k(x, y)$ a transport ray. This means that any optimal transport plan $\Lambda$ moves the mass along the transport rays. An important fact is that two different transport rays cannot intersect at an interior point of one of them; to see that, assume $\gamma^+ := \gamma_{x^+,y^+}$ and $\gamma^- := \gamma_{x^-,y^-}$ are two different transport rays. Without loss of generality, assume that they intersect at $z$, an interior point of $\gamma^-$. The point $z$ either is an interior point of $\gamma^+$ or it is one of its endpoints (by symmetry, we assume it to be $x^+$). We have

$$\psi(z) - \psi(y^+) = d_k(z, y^+) \quad \text{and} \quad \psi(x^-) - \psi(z) = d_k(x^-, z).$$

Then,

$$\psi(x^-) - \psi(y^+) = d_k(x^-, z) + d_k(z, y^+).$$

Yet, $\psi$ is 1-Lipschitz with respect to $d_k$. Hence, $z$ belongs to a geodesic $\gamma$ from $x^-$ to $y^+$, which is a concatenation of fragments of $\gamma^+$ and $\gamma^-$. But, this is a contradiction as there is a unique minimizing geodesic starting at $x^-$ with initial velocity $\gamma'(0)$ (thanks to the $C^{1,1}$ regularity of the Riemannian metric $k$; see [6, Chapter 8] or [13, Chapter 7]).

In optimal transport theory it is classical to associate with any optimal transport plan $\Lambda$ a non-negative measure $\sigma_\Lambda$ on $\bar{\Omega}$, called the *transport density*, which represents the amount of transport taking place in each region of $\Omega$. This measure $\sigma_\Lambda$ is defined as follows

$$(4.4) \quad \langle \sigma_\Lambda, \phi \rangle := \int_{\bar{\Omega} \times \bar{\Omega}} \int_0^1 \phi(\gamma_{x,y}(t)) k(\gamma_{x,y}(t)) |\gamma'_{x,y}(t)| dt d\Lambda(x, y), \text{ for all } \phi \in C(\bar{\Omega}),$$
where $\gamma_{x,y}$ is the unique geodesic between $x$ and $y$. We note that the measure $\sigma_\Lambda$ is well defined thanks to the fact that $\Omega$ is geodesically convex. Denote by $\mathcal{H}_k^1$ the weighted Hausdorff measure, i.e. $\mathcal{H}_k^1(\gamma) := \int_0^1 k(\gamma(t)) |\gamma'(t)| dt$ for a Lipschitz curve $\gamma$. From (4.4), we can see easily that

$$ (4.5) \qquad \sigma_\Lambda(A) := \int_{\Omega \times \bar{\Omega}} \mathcal{H}_k^1(\gamma_{x,y} \cap A) \, d\Lambda(x,y), \quad \text{for all Borel set } A \subset \bar{\Omega}. $$

On the other hand, we define a vector measure $v_\Lambda$, which is the vector version of $\sigma_\Lambda$, as follows

$$ (4.6) \qquad \langle v_\Lambda, \xi \rangle := \int_{\bar{\Omega} \times \bar{\Omega}} \int_0^1 \xi(\gamma_{x,y}(t)) \cdot k(\gamma_{x,y}(t)) \gamma'_{x,y}(t) \, dt \, d\Lambda(x,y), \quad \text{for all } \xi \in C(\bar{\Omega}, \mathbb{R}^2). $$

Let $\psi$ be a Kantorovich potential in the dual problem (4.3). Thanks to [16, Lemma 10], $\psi$ is differentiable at any interior point of a transport ray $\gamma_{x,y}$ and, we have

$$ \gamma'_{x,y}(t) = -k^{-1}(\gamma_{x,y}(t)) d_k(x,y) \frac{\nabla \psi(\gamma_{x,y}(t))}{|\nabla \psi(\gamma_{x,y}(t))|}. $$

This implies that

$$ v_\Lambda = -\sigma_\Lambda \frac{\nabla \psi}{|\nabla \psi|}. $$

In particular, we have $|v_\Lambda| = \sigma_\Lambda$. From (4.4) and the optimality of the transport plan $\Lambda$, we infer that

$$ \int_{\bar{\Omega}} |v_\Lambda| = \sigma_\Lambda(\bar{\Omega}) = \int_{\bar{\Omega} \times \bar{\Omega}} d_k(x, y) \, d\Lambda(x, y) = \min (4.1). $$

Fix $\phi \in C^1(\bar{\Omega})$. Taking $\xi = k^{-1}\nabla\phi \in C(\bar{\Omega}, \mathbb{R}^2)$ as a test function in the definition (4.6) of $v_\Lambda$, then we have

$$ \begin{aligned} \langle v_\Lambda, k^{-1}\nabla\phi \rangle &= \int_{\bar{\Omega} \times \bar{\Omega}} \int_0^1 \nabla\phi(\gamma_{x,y}(t)) \cdot \gamma'_{x,y}(t) \, dt \, d\Lambda(x,y) \\ &= \int_{\bar{\Omega} \times \bar{\Omega}} \int_0^1 \left[ \frac{d}{dt}\phi(\gamma_{x,y}(t)) \right] \, dt \, d\Lambda(x,y) \\ &= \int_{\bar{\Omega} \times \bar{\Omega}} [\phi(y) - \phi(x)] \, d\Lambda(x,y) \\ &= \int_{\bar{\Omega}} \phi \, d(f^{-} - f^{+}), \end{aligned} $$

This implies that $\text{div}[k^{-1}v_\Lambda] = f^{+} - f^{-}$. Now, let us consider the weighted Beckmann problem

$$ (4.7) \qquad \min \left\{ \int_{\bar{\Omega}} k|v| : v \in \mathcal{M}(\bar{\Omega}, \mathbb{R}^2), \text{ div } v = f \right\}. $$

Similarly to the Euclidean case, on geodesically convex domains we have an equivalence between the Beckmann problem and the optimal transport problem. By the reasoning above, we infer that

$$ \min (4.7) \le \int_{\bar{\Omega}} k |k^{-1} v_{\Lambda}| = \int_{\bar{\Omega}} |v_{\Lambda}| = \min (4.1). $$

Yet, it is easy to see that the reverse inequality also holds, since for any vector measure $v$ such that $\text{div}(v) = f$ and any smooth function $\psi$ such that $|\nabla\psi| \le k$, one has

$$ \int_{\bar{\Omega}} \psi \, d(f^{+} - f^{-}) = - \int_{\bar{\Omega}} \nabla \psi \cdot dv \le \int_{\bar{\Omega}} k \, d|v|. $$

Hence, we get that

$$ \min (4.7) = \min (4.1). $$

Moreover, we infer that $k^{-1}v_\Lambda$ solves the Beckmann problem (4.7). In fact, one can prove (see [15, Proposition 2.3] in the Euclidean case) that every optimal flow $v$ of Problem (4.7) is of the form $v = k^{-1}v_\Lambda$, for some optimal transport plan $\Lambda$ of Problem (4.1). Let us note that in [33] the author proved the equivalence between the Beckmann problem and the optimal transport problem in a more general setting, but then the definition of the transport density has to be altered to allow for splitting of the transport along multiple transport rays. Since in our regularity results we rely on the uniqueness of the geodesic between any two given points, let us present a sketch of the proof of an analogue of [15, Proposition 2.3] in the Riemannian case.
Let $\mathcal{C}$ be the set of absolutely continuous curves $\gamma : [0, 1] \mapsto \bar{\Omega}$. We call a traffic plan any nonnegative measure $\eta$ on $\mathcal{C}$ such that $(e_0)_{\#}\eta = f^+$ and $(e_1)_{\#}\eta = f^-$, where for all $t \in [0, 1]$ we denote $e_t(\gamma) := \gamma(t)$. We define the traffic flow $v_\eta$ and the traffic intensity $i_\eta$ as follows:

$$< v_\eta, \xi > := \int_C \int_0^1 \xi(\gamma(t)) \cdot k(\gamma(t)) \gamma'(t) \, dt \, d\eta(\gamma), \quad \text{for all } \xi \in C(\bar{\Omega}, \mathbb{R}^2),$$

and

$$< i_\eta, \phi > := \int_C \int_0^1 \phi(\gamma(t)) k(\gamma(t)) |\gamma'(t)| \, dt \, d\eta(\gamma), \quad \text{for all } \phi \in C(\bar{\Omega}).$$

We note that $\text{div}[k^{-1}v_\eta] = f^+ - f^-;$ in particular, it is admissible in Problem (4.7). Let $v$ be an optimal flow for Problem (4.7). Similarly to [15, Lemma 2.2], one can show that there is a traffic plan $\eta$ such that

$$\int_{\bar{\Omega}} d|kv - v_{\eta}| + \int_{\bar{\Omega}} d i_{\eta} = \int_{\bar{\Omega}} k \, d|v|.$$

Note that $|v_\eta| \le i_\eta$. Since $k^{-1}v_\eta$ is admissible in Problem (4.7) while $v$ minimizes Problem (4.7), then we get

$$(4.8) \qquad \int_{\bar{\Omega}} d i_{\eta} \ge \int_{\bar{\Omega}} k \, d|v|.$$

Moreover, we have

$$\int_{\bar{\Omega}} k \, d|v| \le \int_{\bar{\Omega}} d|kv - v_{\eta}| + \int_{\bar{\Omega}} d|v_{\eta}| \le \int_{\bar{\Omega}} d|kv - v_{\eta}| + \int_{\bar{\Omega}} d i_{\eta} = \int_{\bar{\Omega}} k \, d|v|.$$

Hence, all the inequalities above are in fact equalities. But then, equation (4.8) implies that

$$v = k^{-1}v_{\eta} \quad \text{and} \quad |v_{\eta}| = i_{\eta}.$$

Set $\Lambda_\eta := (e_0, e_1)_{\#}\eta$. Then, it is clear that $\Lambda_\eta$ is a transport plan between $f^+$ and $f^-$. In addition, we have

$$\int_{\bar{\Omega}} d|v_{\eta}| = \int_C \int_0^1 k(\gamma(t)) |\gamma'(t)| \, dt \, d\eta(\gamma) \ge \int_C d_k(\gamma(0), \gamma(1)) \, d\eta(\gamma) = \int_{\bar{\Omega}} d_k(x, y) \, d\Lambda_{\eta}.$$

Yet, $\min(4.7) = \min(4.1)$ and $k^{-1}v_\eta$ minimizes Problem (4.7). Hence, the inequality above is an equality, which implies that for $\eta$-a.e. $\gamma \in \mathcal{C}$, $d_k(\gamma(0), \gamma(1)) = \int_0^1 k(\gamma(t))|\gamma'(t)| \, dt$ and so, $\gamma$ is the unique geodesic between $\gamma(0)$ and $\gamma(1)$. Consequently, $v = k^{-1}v_\eta = k^{-1}v_{\Lambda_\eta}$ and $\Lambda_\eta$ is an optimal transport plan for Problem (4.1).

The $L^p$ summability of the transport density $\sigma_\Lambda$ was already considered in many papers [9–11, 35]. However, in all these works, the authors just consider the Euclidean case (i.e., $k \equiv 1$). More precisely, they show that for all $p \in [1, \infty]$ the transport density $\sigma_\Lambda$ belongs to $L^p(\Omega)$ as soon as $f^+ \in L^p(\Omega)$. Moreover, the authors of [15] have already considered the case where the source and target measures are singular ($f^+$ are two measures concentrated on $\partial\Omega$): they show (again, under the assumption that $k \equiv 1$) that if $f^+ \in L^p(\partial\Omega)$, then the transport density $\sigma_\Lambda$ between them is in $L^p(\Omega)$ provided that $p \le 2$ and $\Omega$ is uniformly convex. In this paper, we will study the $L^p$ summability of the transport density $\sigma_\Lambda$ in the case where $f^+$ are concentrated on the boundary $\partial\Omega$ and $k$ is a positive $C^{1,1}$ Riemannian metric.

Let us come back to the question of existence and uniqueness of an optimal transport map $T$ for the Monge problem (4.2). We note that this map $T$ (when it exists) is not necessarily unique and, the Kantorovich problem (4.1) may have many different solutions as well. Under the assumption that $f^+ \in L^1(\Omega)$, one can pick a special optimal transport plan $\Lambda$ which will be induced by a transport map $T$, and this map turns out to be a solution for the Monge problem (4.2). Yet, in this paper, we consider the case where $f^+$ and $f^-$ are two nonnegative measures which are concentrated on $\partial\Omega$ and so, it is not clear if the Monge problem (4.2) reaches a minimum or not. Before studying existence and uniqueness of minimizers for Problems (4.1) & (4.2), we need to introduce the following:

**Definition 4.1.** We say that $\Omega$ is geodesically strictly convex if given any two points in $\bar{\Omega}$, the unique minimizing geodesic between them lies in the interior of $\Omega$, possibly except for its endpoints.
Then, we have the following existence and uniqueness result (the proof is similar to the one given in [15, Proposition 2.5], but for the sake of completeness we will introduce here its adaptation to the Riemannian case):

**Proposition 4.2.** Assume that $\Omega$ is geodesically strictly convex. Then, there is a unique optimal transport plan $\Lambda$, between $f^+$ and $f^-$, which will be induced by a transport map $T$, provided that $f^+$ is non-atomic. This map $T$ turns out to be the unique optimal transport map in the Monge problem (4.2).

*Proof.* Let $\Lambda$ be an optimal transport plan between $f^+$ and $f^-$. Let us denote by $D$ the set of points whose belong to different transport rays. Since two transport rays cannot intersect at an interior point of either of them, we have that $D \subset \partial\Omega$. Fix $x \in D$ and let $\gamma_x^\pm$ be two different transport rays starting at $x$. Let $E_x$ be the region delimited by $\gamma_x^+$, $\gamma_x^-$ and $\partial\Omega$. Then, we easily see that these sets $\{E_x\}_{x \in D}$ must be essentially disjoint with $\mathcal{L}^2(E_x) > 0$, for every $x \in D$. This implies that the set $D$ is at most countable. Yet, $f^+$ is non-atomic and so, $f^+(D) = 0$. On the other hand, we have that for $f^+-$almost every $x \notin D$, there is a unique transport ray $\gamma_x$ starting at $x$ and, thanks to the fact that $\Omega$ is geodesically strictly convex, this geodesic $\gamma_x$ intersects $\text{spt}(f^-)$ at exactly one point $T(x)$. This implies that $\Lambda = (Id, T)_{\#}f^+$. Yet, this is sufficient to infer that $\Lambda$ is the unique optimal transport plan in the Kantorovich problem (4.1) since, if $\Lambda' = (Id, T')_{\#}f^+$ is another optimal transport plan then $\Lambda'' = (\Lambda + \Lambda')/2$ is also optimal and so, $\Lambda''$ must be induced by a transport map $T''$. But, this yields to a contradiction as soon as $T \neq T'$, since it is not possible to have the following equality for all $\varphi \in C(\bar{\Omega} \times \bar{\Omega})$:

$$ \int_{\bar{\Omega}} \varphi(x, T''(x)) \, df^+(x) = \frac{1}{2} \left[ \int_{\bar{\Omega}} \varphi(x, T(x)) \, df^+(x) + \int_{\bar{\Omega}} \varphi(x, T'(x)) \, df^+(x) \right]. \quad \square $$

<figure><img src="image_2.png" /></figure>

Fix $\tau \in [0, 1]$. Then, we define the partial transport densities $\sigma_\Lambda^+$ and $\sigma_\Lambda^-$ as follows:

$$ (4.9) \qquad < \sigma_{\Lambda}^{+}, \phi > := \int_{\bar{\Omega} \times \bar{\Omega}} \int_{0}^{\tau} \phi(\gamma_{x,y}(t)) k(\gamma_{x,y}(t)) |\gamma'_{x,y}(t)| dt d\Lambda(x,y), \quad \text{for all } \phi \in C(\bar{\Omega}), $$

and

$$ (4.10) \qquad < \sigma_{\Lambda}^{-}, \phi > := \int_{\bar{\Omega} \times \bar{\Omega}} \int_{\tau}^{1} \phi(\gamma_{x,y}(t)) k(\gamma_{x,y}(t)) |\gamma'_{x,y}(t)| dt d\Lambda(x,y), \quad \text{for all } \phi \in C(\bar{\Omega}). $$

It is clear that $\sigma_\Lambda = \sigma_\Lambda^+ + \sigma_\Lambda^-$. To prove $L^p$ estimates on the transport density $\sigma_\Lambda$, the idea in the next sections will be to prove $L^p$ summability on $\sigma_\Lambda^\pm$.

In what follows, we will also need the following stability result.

**Proposition 4.3.** Assume that $f^\pm \in \mathcal{M}^+(\partial\Omega)$ and let $(f_n^-)_n \subset \mathcal{M}^+(\partial\Omega)$ be such that $f_n^- \to f^-$. Let $\Lambda_n$ be an optimal transport plan between $f^+$ and $f_n^-$. Then, $\Lambda_n \to \Lambda$, where $\Lambda$ is an optimal transport plan between $f^+$ and $f^-$. Moreover, $\sigma_{\Lambda_n}^\pm \to \sigma_\Lambda^\pm$. Finally, if $\Lambda_n := (I, T_n)_{\#}f^+$ and $\Lambda := (I, T)_{\#}f^+$, then (up to taking a subsequence) for $f^+-$a.e. $x \in \partial\Omega$, we have $T_n(x) \to T(x)$.

*Proof.* Let $\psi_n$ be a Kantorovich potential between $f^+$ and $f_n^-$. Then, from the duality min (4.1) = sup (4.3), we have

$$ \int_{\bar{\Omega} \times \bar{\Omega}} d_k(x, y) \, d\Lambda_n(x, y) = \int_{\bar{\Omega}} \psi_n \, d(f^+ - f_n^-). $$

Yet, it is clear that, up to a subsequence, $\psi_n \to \psi$ uniformly. So, passing to the limit when $n \to \infty$, we get

$$ \int_{\bar{\Omega} \times \bar{\Omega}} d_k(x, y) \, d\Lambda(x, y) = \int_{\bar{\Omega}} \psi \, d(f^+ - f^-). $$

This implies that $\Lambda$ is an optimal transport plan between $f^+$ and $f^-$, while $\psi$ is a Kantorovich potential between them. The second statement follows directly from the definitions (4.9) & (4.10) of $\sigma_\Lambda^\pm$. The last statement follows directly from the fact that

$$ \int_{\bar{\Omega}} \xi(x, T_n(x)) \, df^+(x) \to \int_{\bar{\Omega}} \xi(x, T(x)) \, df^+(x), \quad \text{for all } \xi \in C(\bar{\Omega} \times \bar{\Omega}). $$

Take $\xi(x, y) = \xi^+(x) \cdot y$ and $\xi(x, y) = |y|^2$ respectively, we infer that $T_n \to T$ in $L^2$ and $\|T_n\|_{L^2} \to \|T\|_{L^2}$ and so, $T_n \to T$ in $L^2$. In particular, it converges $f^+-$a.e. on a subsequence. $\square$

<figure><img src="image_3.png" /></figure>
We will conclude this section by commenting on the relationship between geodesic convexity of the domain and the barrier condition (see Definition 2.7), which is the assumption used for existence of solutions in the least gradient problem in [23]. It turns out that the assumption of geodesic strict convexity of $\Omega$, required to use optimal transport techniques, implies the barrier condition.

**Proposition 4.4.** Suppose that $\Omega$ is geodesically strictly convex. Then, it satisfies the barrier condition.

*Proof.* Let $x_0 \in \partial\Omega$. Because $\Omega$ has Lipschitz boundary, we may choose $r > 0$ small enough so that $\partial\Omega \cap \partial B(x_0, r)$ consists of two points; we denote them by $x_1$ and $x_2$. Denote by $\gamma_{x_1x_2}$ the geodesic from $x_1$ to $x_2$; because $\Omega$ is geodesically strictly convex, we have $\gamma_{x_1x_2} \subset \Omega$.

Now, we consider the problem (2.2) in the definition of the barrier condition. Using the direct method of the calculus of variations, we easily see that it admits a solution $V$. The points $x_1, x_2$ separate $\partial\Omega$ into two arcs: $\Gamma_1$ is the arc which contains $x_0$ and $\Gamma_2$ is the arc which does not contain $x_0$. On the other hand, notice that $\gamma_{x_1x_2}$ separates $\Omega$ into two regions; let us denote by $V'$ the region whose closure contains $\Gamma_2$. Hence, we may take $V'' = V' \cup (\Omega \setminus B(x_0, r))$ as a competitor in (2.2) and see that $V = V''$; this happens because $V'$ minimizes perimeter in a larger class, namely it is a solution of

$$ (4.11) \qquad \inf \left\{ P_{\phi}(W, \mathbb{R}^2) : W \subset \Omega, T\chi_W \ge \chi_{\Gamma_2} \right\}, $$

so as $V$ is admissible in (2.2), then it is also admissible in (4.11); in particular, $P_{\phi}(V', \mathbb{R}^2) \le P_{\phi}(V, \mathbb{R}^2)$. This inequality still holds after taking a union with $\Omega \setminus B(x_0, r)$, see for instance [26] (or simply notice that $\partial V''$ is a Lipschitz curve which arises from replacing a part of $\partial V'$ in case when it leaves $B(x_0, r)$ with a part of $\partial B(x_0, r)$). In particular, $\partial V$ consists of parts of $\gamma_{x_1,x_2}$ and parts of $\partial B(x_0, r)$, so $\partial V \cap \partial \Omega \cap B(x_0, r) = \emptyset$ and the barrier condition is satisfied. $\square$

On the other hand, the barrier condition is only a local property near $\partial\Omega$ and does not imply geodesic convexity of $\Omega$, because the geodesic between any two given points may fail to be unique. This can be seen by considering the examples in [25] (in this regard, they may be modified for a smooth weight). However, the barrier condition means that the boundary is not locally area-minimizing with respect to internal variations, so once we know that the geodesic between any two given points in $\bar{\Omega}$ is unique, it means that the domain is strictly geodesically convex.

## 5. STRUCTURE OF SOLUTIONS

In this section, we extend the relationship between the weighted least gradient problem (wLGP) and the Kantorovich problem with Riemannian cost (4.1) to their respective dual problems. Namely, we study the relationship between the maximization problem (see [32])

$$ (5.1) \qquad \sup \left\{ \int_{\partial\Omega} [\mathbf{z}, n] g \, d\mathcal{H}^1 : \mathbf{z} \in \mathcal{Z} \right\}, $$

where $g \in BV(\partial\Omega)$ and

$$ \mathcal{Z} = \left\{ \mathbf{z} \in L^{\infty}(\Omega, \mathbb{R}^2), \ \text{div} \mathbf{z} = 0, \ |\mathbf{z}(x)| \le k(x) \text{ for a.e. } x \in \Omega \right\}, $$

with the maximization problem (4.3)

$$ \sup \left\{ \int_{\bar{\Omega}} \psi \, d(f^+ - f^-) : |\psi(x) - \psi(y)| \le d_k(x, y), \forall x, y \in \bar{\Omega} \right\}. $$

Here, the normal trace $[\mathbf{z}, n]$ is understood in the weak sense (see Section 2.1). Using a standard reasoning in duality theory, one can see that both problems admit solutions; see [32] for problem (5.1) and [36, 39] for problem (4.3). Since the infimal values in the primal problems are equal, the supremal values in the dual problems also coincide. We now show that these problems are equivalent in the sense that from a solution of one problem we may recover a solution of the other problem. In the Euclidean case, such a result was shown in [20]. We use a similar technique in the weighted case; the main differences are that instead of uniform bounds on $\mathbf{z}$ we now have pointwise bounds and the domain is no longer convex, so we need to rely on Proposition 3.1 (valid for any contractible domain) to recover the Kantorovich potential.
**Theorem 5.1.** Suppose that $\Omega$ is geodesically convex. Then, the problems (4.3) and (5.1) are equivalent in the following sense:

(1) Their supremal values coincide, i.e. $\sup(4.3) = \sup(5.1)$.

(2) Given a maximizer $\psi \in \text{Lip}(\bar{\Omega})$ of (4.3), we can construct a maximizer $\mathbf{z} \in L^{\infty}(\Omega, \mathbb{R}^2)$ of (5.1); moveover, $\mathbf{z} = R_{-\frac{\pi}{2}}\nabla\psi$ in $\Omega$.

(3) Given a maximizer $\mathbf{z} \in L^{\infty}(\Omega, \mathbb{R}^2)$ of (5.1), we may construct a maximizer $\psi \in \text{Lip}(\bar{\Omega})$ of (4.3) and, we have $\mathbf{z} = R_{-\frac{\pi}{2}}\nabla\psi$ in $\Omega$.

Notice that the direction of the rotation is opposite to the direction of rotation in Theorem 2.11. We need geodesic convexity of $\Omega$, because in the computation of supremal values we need to pass through the primal problems, and the equivalence for primal problems holds on geodesically convex domains.

*Proof.* (1) This follows immediately from the equivalence between the weighted least gradient problem, weighted Beckmann problem, and the Monge-Kantorovich problem. We have

$$ \sup(5.1) = \inf(\text{wLGP}) = \inf(4.1) = \sup(4.3). $$

(2) Suppose that $\psi \in \text{Lip}(\bar{\Omega})$ is a maximizer in Problem (4.3). Take $\mathbf{z} = R_{-\frac{\pi}{2}}\nabla\psi \in L^{\infty}(\Omega, \mathbb{R}^2)$. Then, $\mathbf{z}$ is admissible in (5.1), since $\text{div}(\mathbf{z}) = 0$ as distributions and $|\mathbf{z}(x)| = |\nabla\psi(x)| \le k(x)$ a.e. in $\Omega$. Since $\psi$ is Lipschitz and $g \in BV(\partial\Omega)$, we have $\psi g \in BV(\partial\Omega)$ and a mass balance condition holds:

$$ (5.2) \qquad 0 = \int_{\partial\Omega} d[\partial_{\tau}(\psi g)] = \int_{\partial\Omega} \psi d(\partial_{\tau}g) + \int_{\partial\Omega} (\partial_{\tau}\psi) g d\mathcal{H}^1. $$

Now, $\mathbf{z} = R_{-\frac{\pi}{2}}\nabla\psi = -R_{\frac{\pi}{2}}\nabla\psi$, so by Proposition 3.1 we have $[\mathbf{z}, n] = -\partial_{\tau}\psi$. By (5.2),

$$ \sup(5.1) = \sup(4.3) = \int_{\partial\Omega} \psi df = \int_{\partial\Omega} \psi d(\partial_{\tau}g) = \int_{\partial\Omega} (-\partial_{\tau}\psi) g d\mathcal{H}^1 = \int_{\partial\Omega} [\mathbf{z}, n] g d\mathcal{H}^1. $$

Hence, $\mathbf{z}$ is a maximizer for Problem (5.1).

(3) Suppose that $\mathbf{z} \in L^{\infty}(\Omega, \mathbb{R}^2)$ is a maximizer in (5.1). By Proposition 3.1 there exists $\overline{\psi} \in W^{1,1}(\Omega)$ such that $\mathbf{z} = R_{\frac{\pi}{2}}\nabla\overline{\psi}$. Since $\mathbf{z} \in L^{\infty}(\Omega, \mathbb{R}^2)$ with $|\mathbf{z}(x)| \le k(x)$ a.e. in $\Omega$, we also have $\nabla\overline{\psi} \in L^{\infty}(\Omega, \mathbb{R}^2)$ with $|\nabla\overline{\psi}(x)| \le k(x)$ a.e. in $\Omega$. Hence, $\overline{\psi}$ is 1-Lipschitz with respect to the geodesic distance $d_k$, so it is admissible in Problem (4.3). Moreover, Proposition 3.1 implies that $[\mathbf{z}, n] = \partial_{\tau}\overline{\psi}$. Again using equation (5.2), we get

$$ \sup(4.3) = \sup(5.1) = \int_{\partial\Omega} [\mathbf{z}, n] g d\mathcal{H}^1 = \int_{\partial\Omega} (\partial_{\tau}\overline{\psi}) g d\mathcal{H}^1 = - \int_{\partial\Omega} \overline{\psi} d(\partial_{\tau}g) = \int_{\partial\Omega} (-\overline{\psi}) df. $$

Hence, $\psi = -\overline{\psi}$ is a maximizer for Problem (4.3). In particular, $\mathbf{z} = R_{\frac{\pi}{2}}\nabla\overline{\psi} = R_{-\frac{\pi}{2}}\nabla\psi$. $\square$

<figure><img src="image_4.png" /></figure>

Hence, we may express the solution to the dual of the weighted least gradient problem (5.1) via a Kantorovich potential of the corresponding optimal transport problem and vice versa. In general, we cannot expect solutions to any of these problems to be unique. However, since the Kantorovich potentials are differentiable in the interiors of the transport rays and their gradient is uniquely defined, their structure is somewhat prescribed by the boundary data. We will now study some consequences of this result for the structure of solutions to the weighted least gradient problem. For the remainder of this section, for $u \in BV(\Omega)$, we denote $E_t := \{u \ge t\}$.

**Corollary 5.2.** Suppose that $\Omega$ is geodesically convex. Let $u$ be a solution to the weighted least gradient problem with boundary data $g \in BV(\partial\Omega)$. Then, for all $t \in \mathbb{R}$, every connected component of $\partial E_t$ is a transport ray.

This result is a precise statement of the informal observation that boundaries of superlevel sets of a solution to the least gradient problem are transport rays, which was made in the Euclidean case in [15].

*Proof.* By Theorem 2.10, characteristic functions of superlevel sets $\chi_{E_t}$ are functions of least gradient, for all $t \in \mathbb{R}$. By Theorem 2.9, there exists a vector field $\mathbf{z} \in L^{\infty}(\Omega, \mathbb{R}^2)$ with $|\mathbf{z}| \le k$ and zero divergence such that $(\mathbf{z}, Du) = k|Du|$; in particular, $\mathbf{z}$ maximizes (5.1) since $u$ is a solution for Problem (wLGP). Notice that by the co-area formula (see [3, Proposition 2.7]) the vector field $\mathbf{z}$ satisfies $(\mathbf{z}, D\chi_{E_t}) = k|D\chi_{E_t}|$; technically, after the application of the co-area formula this result holds for almost all $t$, but since the vector field $\mathbf{z}$ is fixed we may approximate $\chi_{E_t}$ with $\chi_{E_{t_n}}$ in $L^1$
norm and obtain the result for all $t$. Now, since the distributional derivative of $\chi_{E_t}$ is concentrated on its jump set ($|D\chi_{E_t}| = \mathcal{H}_{|\partial E_t|}^1$), then by [8, Theorem 3.3], we have $[z, n^{E_t}] = k \mathcal{H}^1$-a.e. on $\partial E_t$; if we understand by $z$ its precise representative, it means that $z = k \cdot n^{E_t} \mathcal{H}^1$-a.e. on $\partial E_t$. Hence, by Theorem 5.1, there is a Kantorovich potential $\psi$ such that $\nabla\psi = k(n^{E_t})^\perp \mathcal{H}^1$-a.e. on $\partial E_t$, so every connected component of $\partial E_t$ is in fact a transport ray. $\square$

For the remainder of this section, let us introduce the following notation. We will work under the assumption that $\Omega$ is geodesically convex. For two points $x, y \in \bar{\Omega}$, denote by $\gamma_{x,y}$ the unique geodesic between these two points. Given $g \in BV(\partial\Omega)$, we denote by $D$ the set of discontinuity points of $g$, which is at most a countable set. We also denote by $\mathcal{D}$ the union of all geodesics between points in $D$. Since there is exactly one geodesic between two given points in $D$, the set $\mathcal{D}$ is a union of at most countably many geodesics and it has Hausdorff dimension one (possibly with infinite measure).

The following several lemmas are steps in the proof of Theorem 5.6, but are also of independent interest as results on the boundary behaviour and the jump set of a solution. Let us note that under the stronger assumption that $\Omega$ is geodesically strictly convex, it satisfies the barrier condition, and the first two results follow from the analysis in [19, 32].

**Lemma 5.3.** Suppose that $\Omega$ is geodesically convex. Let $u$ be a solution to the weighted least gradient problem with boundary data $g \in BV(\partial\Omega)$. Then, for all $t \in \mathbb{R}$ and every $x \in \partial\Omega$, there is at most one connected component $R_x$ of $\partial E_t$ such that $x \in R_x$.

The proof is similar to the one given in [17, Proposition 3.5] in the Euclidean case, but we will present here its adaptation to the Riemannian case in order to highlight the main differences and difficulties. Before we proceed, let us note that since $k \in C^{1,1}(\bar{\Omega})$, geodesic convexity of $\Omega$ implies that its boundary cannot have corners with angle larger than $\pi$: the geodesics between any two points in a small neighbourhood on $\partial\Omega$ of a corner point $x_0$ are uniformly $C^{1,1}$ and contained in $\bar{\Omega}$.

*Proof.* Suppose that there are multiple connected components of $\partial E_t$ ending at $x$. Since all of them are transport rays, they cannot intersect in $\Omega$. Therefore, any connected component $R$ of $\partial E_t$ splits $\Omega$ into two open sets, and any other connected component $R'$ of $\partial E_t$ lies entirely in one of these sets. Let $R_1$ and $R_2$ be two adjacent connected components (if it is not possible to find adjacent connected components of $\partial E_t$, because between any given two there are infinitely many, then the perimeter of $E_t$ is infinite, which contradicts Theorem 2.10). Take any $x_1 \in R_1 \cap \Omega$ and $x_2 \in R_2 \cap \Omega$. Consider the curvilinear triangle $\Delta_{x_1,x_2}^x$ with sides $\gamma_{x,x_1}$, $\gamma_{x,x_2}$ and $\gamma_{x_1,x_2}$. It does not reduce to a one-dimensional structure: by geodesic convexity of $\Omega$ and uniqueness of geodesic from $x$ with given initial velocity, $R_1$ and $R_2$ make a positive angle with each other and with the boundary $\partial\Omega$ at $x$. Therefore, the angle between $R_1$ and $R_2$ is smaller than $\pi$, so the concatenation of $\gamma_{x,x_1}$ and $\gamma_{x,x_2}$ cannot be a minimizing curve from $x_1$ to $x_2$; otherwise, this contradicts the regularity of Riemannian geodesics. But then, we can modify the set $E_t$ by adding to it (or removing from it, depending on the position of $R_1$ and $R_2$ with respect to $E_t$) the curvilinear triangle $\Delta_{x_1,x_2}^x$, without changing the trace of $\chi_{E_t}$ and reducing the perimeter of $E_t$. This again contradicts Theorem 2.10, so there can be only one connected component $R_x$ of $\partial E_t$ such that $x \in R_x$. $\square$

**Lemma 5.4.** Suppose that $\Omega$ is geodesically convex. Let $u$ be a solution to the weighted least gradient problem with boundary data $g \in BV(\partial\Omega)$. Then, for all $t \in \mathbb{R}$, we have $\partial E_t \cap \partial\Omega \subset g^{-1}(t) \cup D$.

*Proof.* Suppose that $x \in \partial E_t \cap (\partial\Omega \setminus D)$. Denote by $R_x$ the unique (due to Lemma 5.3) connected component of $\partial E_t$ passing through $x$. Then, by geodesic convexity of $\Omega$ and uniqueness of geodesic from $x$ with given initial velocity, $R_x$ makes a positive angle with $\partial\Omega$ at $x$. Yet, $g$ is continuous at $x$ and so, we may proceed in the same way as in the proof of Lemma [21, Lemma 3.8] (see also [31, Lemma 3.4]) to obtain that $g(x) = t$. $\square$

**Lemma 5.5.** Suppose that $\Omega$ is geodesically convex. Let $u$ be a precise representative of a solution to the weighted least gradient problem with boundary data $g \in BV(\partial\Omega)$. Then, $u$ is continuous at every point in $\Omega \setminus \mathcal{D}$.

*Proof.* Suppose that $u$ is not continuous at $x \in \Omega$. Then, there exist $t, s \in \mathbb{R}$ such that $t < s$ and $x \in \partial E_t \cap \partial E_s$. However, there is at most one transport ray $R_x$ passing through $x$. Since the connected components of $\partial E_t$ and $\partial E_s$ are transport rays, these sets coincide on the whole
connected component passing through $x$ and are equal to $R_x$. By Lemma 5.4, this is only possible if this transport ray has both ends in the discontinuity set $D$, so $x \in R_x \subset D$. $\square$

The following result is the main result in this section. It is a generalisation of [18, Theorem 1.1] to the weighted least gradient problem. Notice that due to the fact that we use optimal transport methods, our proof is significantly shorter and simpler.

**Theorem 5.6.** Suppose that $\Omega$ is geodesically convex. Let $u, v$ be two solutions to the weighted least gradient problem with boundary data $g$. Then, $u = v$ on $\Omega \setminus (C \cup D)$, where both $u$ and $v$ are locally constant on $C$.

*Proof.* Consider any point $x \in \Omega \setminus D$. By Lemma 5.5, both functions $u$ and $v$ are continuous at $x$. Recall that $E_t = \{u \ge t\}$ and set $F_t = \{v \ge t\}$. There are four possibilities:

1. We have $x \in \partial E_t \cap \partial F_t$. Then, since $u$ and $v$ are continuous at $x$, we have $u(x) = v(x) = t$.

2. We have $x \in \partial E_t \cap \partial F_s$, where $t \neq s$. We argue similarly as in the proof of Lemma 5.5. There is at most one transport ray $R_x$ passing through $x$, and because the connected components of $\partial E_t$ and $\partial F_s$ are transport rays, they need to coincide on the whole connected component passing through $x$ and so, they are equal to $R_x$. But, by Lemma 5.4, this is only possible if the transport ray has both ends in the discontinuity set $D$, in which case by uniqueness of the geodesic between any two given points we have $R_x \subset D$, so $x \in \mathcal{N}$, a contradiction.

3. We have $x \in \partial E_t$, but $x \notin \partial F_s$ for any $s \in \mathbb{R}$. Then, $v$ is constant in a ball $B(x, r)$ with value $s \in \mathbb{R}$. If $t = s$, the claim is proved; suppose otherwise. Since $u$ is continuous at $x$ with $x \in \partial E_t$, there exists an interval $I_\varepsilon = [t - \varepsilon, t + \varepsilon]$ such that for all $\tau \in I_\varepsilon$ a connected component $\Gamma_\tau$ of $\partial E_\tau$ intersects the ball $B(x, r)$. Denote by $x_\tau$ and $y_\tau$ the endpoints of $\Gamma_\tau$ and denote by $C_s$ the connected component of $\{v \equiv s\}$ containing $B(x, r)$. Since transport rays cannot intersect inside $\Omega$, we have $\Gamma_\tau \cap \partial C_s = \emptyset$ for all $\tau \in I_\varepsilon$; then, by checking the trace of $v$, we see that on the arc $(x_{t-\varepsilon}, x_{t+\varepsilon})$ on $\partial\Omega$ we have $g \equiv s$. On the other hand, since $x \notin D$, we may assume that $g$ is continuous at $x_\tau$ and by Lemma 5.4 we have that $g(x_\tau) = t$, a contradiction.

4. We have $x \notin \partial E_t$ for any $t \in \mathbb{R}$ and $x \notin \partial F_s$ for any $s \in \mathbb{R}$. Then, both $u$ and $v$ are constant in a ball $B(x, r)$ (possibly with different constants), so $x \in C$. $\square$

## 6. $L^p$ ESTIMATES ON THE TRANSPORT DENSITY

In this section, we study the $L^p$ summability of the transport density $\sigma$ between $f^+$ and $f^-$. Assume that $\Omega$ is geodesically strictly convex. From Proposition 4.2, we recall that there is a unique optimal transport plan $\Lambda$ and so, the transport density $\sigma$ is unique, provided that $f^+$ or $f^-$ is non-atomic. To prove $L^p$ estimates on $\sigma$, the strategy will be the following: we decompose $\sigma$ into two parts $\sigma^+$ and $\sigma^-$, then we approach the target (resp. source) measure $f^-$ (resp. $f^+$) with atomic measures, and under some geometric assumptions on $\Omega$, we prove that $\sigma^+$ (resp. $\sigma^-$) is in $L^p(\Omega)$ as soon as $f^+$ (resp. $f^-$) belongs to $L^p(\partial\Omega)$ with $p \le 2$ (in order to obtain $L^p$ estimates on $\sigma$ with $p > 2$, we need to assume more regularity on $f^\pm$). We note that this $L^p$ summability of the transport density between two measures $f^+$ and $f^-$ which are in $L^p(\partial\Omega)$ was already proved in the Euclidean case (i.e., $k \equiv 1$) in [15]; more precisely, the authors of [15] proved that if $f^\pm \in L^p(\partial\Omega)$ then $\sigma \in L^p(\Omega)$ provided that $p \le 2$ and $\Omega$ is uniformly convex.

First, we are going to introduce a geometric lemma that will be the key point in the proof of the $L^p$ estimates on the transport density $\sigma$. So, let the manifold $\mathbb{R}^2$ be equipped with the conformal metric $d_k$. Fix $x_0 \in \partial\Omega$ and let $\alpha(s)$, $s \in [-\varepsilon, \varepsilon]$, be an arc in $\partial\Omega$. For each $s \in [-\varepsilon, \varepsilon]$, we denote by $\gamma_s$ the geodesic between $\alpha(s)$ and $x_0$. Let $\nu(s)$ be the initial unit tangent vector to $\gamma_s$ and $\tau(s)$ be the (weighted) length of $\gamma_s$.

**Lemma 6.1.** Assume that $\Omega$ is a bounded and geodesically convex domain, and $k \in C^{1,1}(\bar{\Omega})$ with $0 < k_{\min} \le k \le k_{\max} < \infty$. Then, there exists a constant $C < \infty$ depending only on $\text{diam}(\Omega)$, $k_{\min}$, $k_{\max}$, $||\nabla k||_\infty$ and $||D^2k||_\infty$ such that, for a.e. $s \in \partial\Omega$, we have

$$ (6.1) \qquad \det(D_{(s,t)}\gamma_s(t)) \ge (1-t)^C k^{-1}(s) \tau(s)[\nu(s) \cdot n(s)], \quad \text{for all } t \in [0,1]. $$

*Proof.* For $s \in [-\varepsilon, \varepsilon]$, let $(e_1, e_2)$ be a basis such that $e_1 = \alpha'(s)$ and $e_2 = \nu(s)$. Let us parallel-transport along the geodesic $\gamma_s$ to define a new family of basis $(e_1(t), e_2(t))$. Let $[\cdot, \cdot]$ be the Lie
bracket and $\nabla$ be the Levi-Civita connection. Let us denote by $\exp$ the Riemannian exponential map. For all $s \in [-\varepsilon, \varepsilon]$ and $t \in [0, 1]$, set

$$ \Psi(s, t) := \gamma_s(t) = \exp_{\alpha(s)}[t \tau(s) \nu(s)]. $$

Now, we define the vector fields $J_1$ and $J_2$ as follows:

$$ J_1(s, t) = \frac{d}{d\delta}|_{\delta=0} \Psi(s + \delta, t) \quad \text{and} \quad J_2(s, t) = \frac{d}{d\delta}|_{\delta=0} \Psi(s, t + \delta). $$

Set

$$ J(s, t) = (J_1(s, t), J_2(s, t)) = D\Psi(s, t) \quad \text{and} \quad \mathcal{J}(s, t) = \det[J(s, t)]. $$

Thanks to [39, Theorem 11.3], one can show that this Jacobian $\mathcal{J}$ cannot vanish, except possibly at the endpoints of the geodesic $\gamma_s$. So, we have

$$ \mathcal{J}'(s, t) = \operatorname{tr}[J'(s, t) J(s, t)^{-1}] \mathcal{J}(s, t). $$

Let us denote by $d\Psi$ the differential map of $\Psi$. The fact that $[\partial_{e_1}, \partial_{e_2}] = 0$ implies that $J_1$ and $J_2$ commute, since

$$ [J_1, J_2] = [d\Psi(\partial_{e_1}), d\Psi(\partial_{e_2})] = d\Psi[\partial_{e_1}, \partial_{e_2}] = 0. $$

Then, we have

$$ \nabla_{J_1} J_2 = \nabla_{J_2} J_1. $$

Yet,

$$ J_1(s, t) = J_{11}(s, t)e_1(t) + J_{21}(s, t)e_2(t). $$

Hence,

$$ \nabla_{J_2} J_1 = \nabla_{\gamma_s'} J_1 = J'_{11}(s, t)e_1(t) + J_{11}(s, t)e_1'(t) + J'_{21}(s, t)e_2(t) + J_{21}(s, t)e_2'(t). $$

But, $e_2'(t) = \nabla_{\gamma_s'} \gamma_s' = 0$ since $\gamma_s'$ is a geodesic. In addition, we have $\nabla_{\gamma_s'} e_1(t) = \Gamma_{21}^1 e_1(t) + \Gamma_{21}^2 e_2(t)$, where $\Gamma_{21}^1$ and $\Gamma_{21}^2$ denote the Christoffel symbols. Then, we get

$$ (6.2) \qquad \nabla_{J_2} J_1 = (J'_{11}(s,t) + J_{11}(s,t)\Gamma^1_{21})e_1(t) + (J'_{21}(s,t) + J_{21}(s,t)\Gamma^2_{21})e_2(t). $$

On the other hand, we have

$$ \nabla_{J_1} J_2 = J_{11}(s,t) \nabla_{e_1(t)} J_2 + J_{21}(s,t) \nabla_{e_2(t)} J_2. $$

Let $A(t)$ be the matrix, in the basis $\{e_1(t), e_2(t)\}$, associated with the endomorphism $X \mapsto \nabla_X J_2$. Then, one has

$$ \begin{aligned} & J_{11}(s,t)\nabla_{e_1(t)} J_2 + J_{21}(s,t)\nabla_{e_2(t)} J_2 \\ &= J_{11}(s,t)[A_{11}(t)e_1(t) + A_{21}(t)e_2(t)] + J_{21}(s,t)[A_{12}(t)e_1(t) + A_{22}(t)e_2(t)] \\ &= [J_{11}(s,t)A_{11}(t) + J_{21}(s,t)A_{12}(t)]e_1(t) + [J_{11}(s,t)A_{21}(t) + J_{21}(s,t)A_{22}(t)]e_2(t). \end{aligned} $$

From (6.2), we get

$$ \begin{aligned} & (J'_{11}(s,t) + J_{11}(s,t)\Gamma^1_{21})e_1(t) + (J'_{21}(s,t) + J_{11}(s,t)\Gamma^2_{21})e_2(t) \\ &= [J_{11}(s,t)A_{11}(t) + J_{21}(s,t)A_{12}(t)]e_1(t) + [J_{11}(s,t)A_{21}(t) + J_{21}(s,t)A_{22}(t)]e_2(t). \end{aligned} $$

Hence,

$$ J' = [A - B]J, $$

where

$$ B = \begin{pmatrix} \Gamma_{21}^1 & 0 \\ \Gamma_{21}^2 & 0 \end{pmatrix}. $$

This implies that

$$ \mathcal{J}'(s, t) = \operatorname{tr}[A - B] \mathcal{J}(s, t). $$

A direct computation shows that

$$ \Gamma_{21}^1 = \partial_{e_2} \log[\sqrt{1 - (k(\gamma_s(t))^2 e_1(t) \cdot e_2(t))^2}] = 0, $$

so $\operatorname{tr}[B] = 0$. Moreover, we have $J_2(s,t) = \gamma_s'(t) = -k^{-1}(\gamma_s(t)) \tau(s) \frac{\nabla_{d_k(\gamma_s(t),x_0)}}{\nabla_{d_k(\gamma_s(t),x_0)}}$. Now, it is well known that the distance function $d_k(\cdot, x_0)$ to $x_0$ is locally semiconcave in $\mathbb{R}^2\setminus\{x_0\}$ with
$D^2[d_k(x, x_0)] \le \frac{C}{d_k(x, x_0)}I$ for some constant $C$ depending on diam($\Omega$), $k_{\min}$, $k_{\max}$, $||\nabla k||_{\infty}$ and $||D^2k||_{\infty}$ (see, for instance, [6]). This yields that

$$ \mathrm{tr}[A] \ge \frac{-C}{(1-t)}. $$

Hence,

$$ \mathcal{J}'(s, t) \ge \frac{-C}{(1-t)} \mathcal{J}(s, t). $$

We infer that

$$ \log[\mathcal{J}(s, t)] - \log[\mathcal{J}(s, 0)] \ge C \log(1-t). $$

Finally, we get

$$ \mathcal{J}(s, t) \ge \mathcal{J}(s, 0) (1-t)^C. $$

Yet, $J_1(s, 0) = \alpha'(s)$ and $J_2(s, 0) = k^{-1}(s) \tau(s) \nu(s)$. Hence, $\mathcal{J}(s, 0) = k^{-1}(s) \tau(s)[\nu(s) \cdot n(s)]$. Consequently, we get the following estimate:

$$ \mathcal{J}(s, t) \ge (1-t)^C k^{-1}(s) \tau(s)[\nu(s) \cdot n(s)]. \quad \square $$

In order to prove $L^p$ summability on the transport density $\sigma$, we need also to introduce a definition that generalizes the notion of uniform convexity in the Riemannian case. Let us assume that for all $x, y \in \bar{\Omega}$, there is a unique geodesic $\gamma_{x,y}$ from $x$ to $y$. For $x \in \partial\Omega$, let us denote by $\gamma_x$ a geodesic starting at $x$ with initial unit tangent vector $\nu(x)$ such that the endpoint $\gamma_x(1) \in \partial\Omega$ and by $\tau(x)$ the weighted length of this geodesic $\gamma_x$. It is not difficult to see that a domain $\Omega$ is geodesically convex (resp. geodesically strictly convex) if and only if $\nu(x) \cdot n(x) \ge 0$ (resp. $\nu(x) \cdot n(x) > 0$), for all $x \in \partial\Omega$ and all geodesics $\gamma_x$, where $-n(x)$ is any unit vector in the exterior normal cone to $\Omega$ at $x$. Now, we define the uniform geodesic convexity as follows:

**Definition 6.2.** We say that $\Omega$ is geodesically uniformly convex if there is a constant $c > 0$ such that, for any $x \in \partial\Omega$ and all geodesics $\gamma_x$, $\nu(x) \cdot n(x) \ge c \tau(x)$.

Our main result in this section is the following theorem.

**Theorem 6.3.** Suppose that $\Omega$ is strictly geodesically convex. Then, if $f^+ \in L^1(\partial\Omega)$, the transport density $\sigma$ belongs to $L^1(\Omega)$. Moreover, for all $p \in [1, 2]$, we have $\sigma \in L^p(\Omega)$ as soon as $f^+$ and $f^-$ are both in $L^p(\partial\Omega)$ and $\Omega$ is uniformly geodesically convex.

*Proof.* We will use a similar strategy as in the proof of [15, Proposition 3.1]. First, assume that the target measure $f^-$ is atomic with $n$ atoms $\{x_i : 1 \le i \le n\}$. Let $\Gamma_i$ be the set of points $x$ on spth($f^+$) which are transported to the atom $x_i$ and, let us denote by $\Omega_i$ the set of points on the geodesics between $\Gamma_i$ and $x_i$. As $\Omega$ is geodesically convex, then all these sets $\Omega_i$ are subsets of $\bar{\Omega}$. Moreover, thanks to the fact that the transport rays cannot intersect at an interior point of either of them, the sets $\Omega_i$ are essentially disjoint. Let us decompose the transport density $\sigma$ into two parts $\sigma = \sigma^+ + \sigma^-$, where $\sigma^+$ and $\sigma^-$ are defined as follows: for a fixed $\tau \in [0, 1]$, we set

$$ < \sigma^+, \phi > := \int_{\bar{\Omega} \times \bar{\Omega}} \int_0^\tau \phi(\gamma_{x,y}(t)) d_k(x, y) \, dt \, d\Lambda(x, y), \quad \text{for all } \phi \in C(\bar{\Omega}), $$

and

$$ < \sigma^-, \phi > := \int_{\bar{\Omega} \times \bar{\Omega}} \int_{\tau}^{1} \phi(\gamma_{x,y}(t)) d_k(x, y) \, dt \, d\Lambda(x, y), \quad \text{for all } \phi \in C(\bar{\Omega}). $$

First, we will show $L^p$ estimates on $\sigma^+$. Recalling Proposition 4.2, there is a unique optimal transport map $T$ from $f^+$ to $f^-$. Then, we have

$$ < \sigma^+, \phi > = \int_{\partial\Omega} \int_0^\tau \phi(\gamma_{x,T(x)}(t)) d_k(x, T(x)) \, dt \, df^+(x), \quad \text{for all } \phi \in C(\bar{\Omega}). $$

For every $x \in \Gamma_i$, we have $T(x) = x_i$. Hence, we find that

$$ < \sigma^+, \phi > = \sum_{i=1}^{n} \int_{\Gamma_i} \int_{0}^{\tau} \phi(\gamma_{x,x_i}(t)) d_k(x, x_i) \, dt \, df^+(x), \quad \text{for all } \phi \in C(\bar{\Omega}). $$
Fix $i \in \{1, ..., n\}$ and consider $\sigma_i^+$ the restriction of $\sigma^+$ to $\Omega_i$. Let $\alpha(s)$, where $s \in [-\varepsilon, \varepsilon]$, be a parameterization of $\Gamma_i$. Then, we have

$$< \sigma_i^+, \phi > = \int_{-\varepsilon}^{\varepsilon} \int_{0}^{\tau} \phi(\gamma_s(t)) d_k(\alpha(s), x_i) f^+(\alpha(s)) |\alpha'(s)| dt ds, \text{ for all } \phi \in C(\bar{\Omega}),$$

where $\gamma_s$ denotes the geodesic from the point $\alpha(s)$ to the atom $x_i$. Set $y = \gamma_s(t)$, for every $s \in [-\varepsilon, \varepsilon]$ and $t \in [0, \tau]$. Then, one has

$$< \sigma_i^+, \phi > = \int_{\Omega_i} \phi(y) d_k(\alpha(s), x_i) f^+(\alpha(s)) |\alpha'(s)| \mathcal{J}(s,t)^{-1} dy, \text{ for all } \phi \in C(\bar{\Omega}),$$

where

$$\mathcal{J}(s,t) := \det[D_{(s,t)}\gamma_s(t)].$$

Hence,

$$\sigma_i^+(y) = \frac{|\alpha'(s)| d_k(\alpha(s), x_i) f^+(\alpha(s))}{\mathcal{J}(s,t)}, \quad \text{for a.e. } y \in \Omega_i.$$

Set $\tau(s) = d_k(\alpha(s), x_i)$. Then,

$$||\sigma_i^+||_{L^p(\Omega_i)}^p$$

$$(6.3) \qquad = \int_{\Omega_i} \frac{|\alpha'(s)|^p \tau(s)^p f^+(\alpha(s))^p}{\mathcal{J}(s,t)^p} dy = \int_{-\varepsilon}^{\varepsilon} \int_{0}^{\tau} \frac{|\alpha'(s)|^p \tau(s)^p f^+(\alpha(s))^p}{\mathcal{J}(s,t)^{p-1}} dt ds.$$

For $p = 1$, this gives us the following estimate:

$$||\sigma_i^+||_{L^1(\Omega_i)} \leq C ||f^+||_{L^1(\Gamma_i)}.$$

Hence, after taking a sum with respect to $i$ and $\tau = 1$, we infer that $\sigma = \sigma^+$ is in $L^1(\Omega)$ and we have

$$||\sigma||_{L^1(\Omega)} \leq C ||f^+||_{L^1(\partial\Omega)}.$$

Now, let $\nu(s)$ be the initial tangent vector to $\gamma_s$. From Lemma 6.1, there is a constant $C < \infty$ such that the following estimate holds

$$(6.4) \qquad \mathcal{J}(s,t) \geq C^{-1} (1-t)^C \tau(s) [\nu(s) \cdot n(s)].$$

If $\Omega$ is additionally uniformly geodesically convex (see Definition 6.2), then there is a constant $c > 0$ such that $\nu(s) \cdot n(s) \geq c\tau(s)$. So, we get that

$$\mathcal{J}(s,t) \geq C^{-1} (1-t)^C \tau(s)^2.$$

Thus, we continue the computation in (6.3) and we choose $\tau < 1$. Then, we get that

$$||\sigma_i^+||_{L^p(\Omega_i)}^p \leq C^{p-1} \int_0^\tau \frac{1}{(1-t)^{C(p-1)}} dt \int_{-\varepsilon}^\varepsilon \tau(s)^{2-p} |\alpha'(s)| f^+(\alpha(s))^p ds$$

$$(6.5) \qquad \leq C \int_{\Gamma_i} d_k(x, x_i)^{2-p} f^+(x)^p d\mathcal{H}^1(x).$$

For all $p \in [1, 2]$, this implies that

$$||\sigma_i^+||_{L^p(\Omega_i)} \leq C ||f^+||_{L^p(\Gamma_i)}.$$

Hence,

$$||\sigma^+||_{L^p(\Omega)} \leq C ||f^+||_{L^p(\partial\Omega)}.$$

Now, assume that $f^- \in \mathcal{M}^+(\partial\Omega)$ and let $(f_n^-)_n$ be a sequence of atomic measures such that $f_n^- \to f^-$. Let $\sigma_n$ be the transport density between $f^+$ and $f_n^-$. We have

$$(6.6) \qquad ||\sigma_n^+||_{L^p(\Omega)} \leq C ||f^+||_{L^p(\partial\Omega)}.$$

Thanks to Proposition 4.3, passing to the limit in (6.6) when the number of atoms $n \to \infty$, we infer that the positive measure $\sigma^+$ between $f^+$ and $f^-$ is in $L^p(\Omega)$ as soon as $f^+ \in L^p(\partial\Omega)$. Moreover, $\sigma^+$ satisfies the following estimate:

$$||\sigma^+||_{L^p(\Omega)} \leq C ||f^+||_{L^p(\partial\Omega)}.$$
But now, thanks to the uniqueness of the optimal transport plan $\Lambda$ between $f^+$ and $f^-$, it is not difficult to see that if $f^- \in L^p(\partial\Omega)$, then by approximating the source measure $f^+$ with atomic measures, one can obtain $L^p$ estimates on $\sigma^-$. And, we get that

$$||\sigma^-||_{L^p(\Omega)} \leq C||f^-||_{L^p(\partial\Omega)}.$$

Finally, we infer that

$$||\sigma||_{L^p(\Omega)} \leq C||f||_{L^p(\partial\Omega)}, \quad \text{for every } p \in [1, 2]. \quad \square$$

In addition, under the assumption that the source and target measures $f^\pm$ are smooth enough, one we can prove $L^p$ estimates on the transport density $\sigma$ for large $p$. More precisely, we have the following

**Proposition 6.4.** Suppose that $\Omega$ is uniformly geodesically convex. Assume that $f^\pm \in C^{0,\alpha}(\partial\Omega)$ with $0 < \alpha \leq 1$. Then, the transport density $\sigma$ between $f^+$ and $f^-$ belongs to $L^p(\Omega)$ with $p = \frac{2}{1-\alpha}$.

*Proof.* Let $T^+$ be the optimal transport map from $f^+$ to $f^-$ and set $T^- := [T^+]^{-1}$. From (6.5) and thanks to Proposition 4.2, we have the following estimate:

$$||\sigma||_{L^p(\Omega)}^p \leq C \left[ \int_{\Gamma} d_k(x, T^+(x))^{2-p} f^+(x)^p d\mathcal{H}^1(x) + \int_{\Gamma} d_k(T^-(y), y)^{2-p} f^-(y)^p d\mathcal{H}^1(y) \right].$$

As $f^\pm \in C^{0,\alpha}(\Omega)$, then we get that

$$f^{\pm}(x) = f^{\pm}(x) - f^{\pm}(T^{\pm}(x)) \leq C|x - T^{\pm}(x)|^{\alpha} \leq C d_k(x, T^{\pm}(x))^{\alpha}.$$

Hence,

$$||\sigma||_{L^p(\Omega)}^p \leq C^p \int_{\Gamma} d_k(x, T^+(x))^{2-p+p\alpha} d\mathcal{H}^1(x) < \infty$$

as soon as $p \leq 2/(1-\alpha)$. $\square$

On the other hand, it is possible to prove $L^p$ estimates on the transport density $\sigma$ for large $p$ (i.e., $p > 2$) if $f^\pm \in L^p(\partial\Omega)$ but under the assumption that $\text{spt}(f^+)$ and $\text{spt}(f^-)$ do not intersect.

**Proposition 6.5.** Suppose that $\Omega$ is strictly geodesically convex. Let $f^\pm \in L^p(\partial\Omega)$ be such that $\text{spt}(f^+) \cap \text{spt}(f^-) = \emptyset$. Then, the transport density $\sigma$ belongs to $L^p(\Omega)$, for all $p \in [1, \infty]$.

*Proof.* The proof follows directly from the previous estimates in Theorem 6.3. From (6.3), we have

$$||\sigma^+||_{L^p(\Omega)}^p \leq \int_{-\varepsilon}^{\varepsilon} \int_{0}^{\tau} \frac{|\alpha'(s)|^p \tau(s)^p f^+(\alpha(s))^p}{\mathcal{J}(s, t)^{p-1}} dt \, ds.$$

From (6.4) and thanks to the strict geodesic convexity of $\Omega$, we infer that there is a uniform constant $C < \infty$ depending only on the geometry of the domain $\Omega$ such that the following estimate holds

$$\mathcal{J}(s, t) \geq C^{-1}(1 - t)^C.$$

For $\tau < 1$, this implies that

$$||\sigma^+||_{L^p(\Omega)} \leq C||f^+||_{L^p(\partial\Omega)}.$$

In the same way, we prove that

$$||\sigma^-||_{L^p(\Omega)} \leq C||f^-||_{L^p(\partial\Omega)}.$$

Consequently, we get that the transport density $\sigma$ between $f^+$ and $f^-$ is in $L^p(\Omega)$. Moreover, we have the following estimate:

$$||\sigma||_{L^p(\Omega)} \leq C||f||_{L^p(\partial\Omega)}.$$
## 7. APPLICATIONS TO THE WEIGHTED LEAST GRADIENT PROBLEM

Thanks to the results in the previous sections, we prove existence and uniqueness of a solution $u$ to the weighted least gradient problem. Moreover, we show $W^{1,p}$ regularity on this solution $u$ from the $L^p$ estimates on the transport density $\sigma$. First, we have the following

**Theorem 7.1.** Assume that $\Omega$ is strictly geodesically convex. Then, for any boundary datum $g \in BV(\partial\Omega)$, the weighted least gradient problem (wLGP) admits a solution.

*Proof.* Set $f = \partial_{\tau}g$. Let $\Lambda$ be an optimal transport plan between $f^+$ and $f^-$. Let $\sigma_{\Lambda}$ be the transport density associated with $\Lambda$ and let $v_{\Lambda}$ be the vector version of $\sigma_{\Lambda}$ (see (4.6)). Then,

$$ \sigma_{\Lambda}(\partial\Omega) = \int_{\partial\Omega \times \partial\Omega} \mathcal{H}_{k}^{1}(\gamma_{x,y} \cap \partial\Omega) d\Lambda(x, y), $$

where $\mathcal{H}_{k}^{1}(\gamma) := \int_{0}^{1} k(\gamma(t))|\gamma'(t)|dt$. Since $\Omega$ is strictly geodesically convex, for any two points $x, y \in \partial\Omega$, the geodesic $\gamma_{x,y}$ lies in the interior of $\Omega$ except for its endpoints, so $\mathcal{H}_{k}^{1}(\gamma_{x,y} \cap \partial\Omega) = 0$. This implies that $\sigma_{\Lambda}(\partial\Omega) = 0$. From Theorem 3.4, we infer that there is a function $u \in BV(\Omega)$ such that $v_{\Lambda} = R_{\frac{\pi}{2}}Du$ and it is a solution for Problem (wLGP). $\square$

**Theorem 7.2.** Assume that $\Omega$ is strictly geodesically convex and $g \in BV(\partial\Omega) \cap C(\partial\Omega)$. Then, the weighted least gradient problem (wLGP) has a unique solution.

*Proof.* Since $g \in BV(\partial\Omega) \cap C(\partial\Omega)$, its tangential derivative is a non-atomic measure and so, by Proposition 4.2, the problem (4.1) has a unique optimal transport plan $\Lambda$. From the equivalence between Problems (4.1) and (4.7), we infer that $v_{\Lambda}$ is the unique minimizer of Problem (4.7). Thanks to Theorem 3.4, if $u$ is any solution for Problem (wLGP) then the measure $v := R_{\frac{\pi}{2}}Du$ solves Problem (wBP). This implies that also the solution to Problem (wLGP) is unique. $\square$

The main results in this paper are the following $W^{1,p}$ regularity estimates on the solution $u$ of Problem (wLGP):

**Theorem 7.3.** Suppose that $\Omega$ is strictly geodesically convex and $g \in W^{1,1}(\partial\Omega)$. Then, the solution $u$ of Problem (wLGP) belongs to $W^{1,1}(\Omega)$.

*Proof.* If $g \in W^{1,1}(\partial\Omega)$, then $f = \partial_{\tau}g \in L^1(\partial\Omega)$. From Theorem 6.3, we infer that the transport density $\sigma$ between $f^+$ and $f^-$ is in $L^1(\Omega)$. But, thanks to Theorem 3.4, we have $|Du| = \sigma$, so $u$ belongs to $W^{1,1}(\Omega)$. $\square$

**Theorem 7.4.** Suppose that $\Omega$ is uniformly geodesically convex and $g \in W^{1,p}(\partial\Omega)$ with $p \in [1, 2]$. Then, the solution $u$ of Problem (wLGP) is in $W^{1,p}(\Omega)$.

*Proof.* If $g \in W^{1,p}(\partial\Omega)$, then $f = \partial_{\tau}g \in L^p(\partial\Omega)$. From Theorem 6.3, we infer that the transport density $\sigma$ between $f^+$ and $f^-$ is in $L^p(\Omega)$. So, using Theorem 3.4, we have $|Du| = \sigma$ and we get $u \in W^{1,p}(\Omega)$. $\square$

**Theorem 7.5.** Assume that $\Omega$ is uniformly geodesically convex and $g \in C^{1,\alpha}(\partial\Omega)$ with $\alpha \in (0, 1]$. Then, the solution $u$ of Problem (wLGP) is in $W^{1,p}(\Omega)$ with $p = 2/(1 - \alpha)$.

*Proof.* If $g \in C^{1,\alpha}(\partial\Omega)$, then $f = \partial_{\tau}g \in C^{0,\alpha}(\partial\Omega)$. Then, the result follows immediately from Proposition 6.4. $\square$

**Theorem 7.6.** Suppose that $\Omega$ is strictly geodesically convex. Assume that $g \in W^{1,p}(\partial\Omega)$ and $g$ has flat parts separating those with positive and negative derivatives. Then, the solution $u$ of Problem (wLGP) belongs to $W^{1,p}(\Omega)$, for all $p \in [1, \infty]$.

*Proof.* Notice that $\text{spt}(f^+) \cap \text{spt}(f^-) = \emptyset$. If $g \in W^{1,p}(\partial\Omega)$, then $f = \partial_{\tau}g \in L^p(\partial\Omega)$, so the proof follows immediately from Proposition 6.5 and the fact that $|Du| = \sigma$. $\square$

We conclude with a generalisation of [20, Theorem 4.1] from the Euclidean case. Let us note that in this case the solution is not necessarily unique, but the result is valid for all the solutions.

**Theorem 7.7.** Suppose that $\Omega \subset \mathbb{R}^2$ is strictly geodesically convex. Let $g \in SBV(\partial\Omega)$. If $u \in BV(\Omega)$ is a solution to Problem (2.3) with boundary data $g$, then $u \in SBV(\Omega)$.
*Proof.* The proof of this result in the Euclidean case in [20] can shortly be described as follows. Since $g \in SBV(\partial\Omega)$, the measure $f = \partial_\tau g$ has no Cantor part, so $f^\pm$ is a union of an absolutely continuous part $f_{ac}^\pm$ and an atomic part $f_{at}^\pm$. We split the transport density into four parts, corresponding to the transport between each of the $f_{ac}^\pm$ and $f_{at}^\pm$. The transport density between the atomic parts is concentrated on a set of Hausdorff dimension one. The other parts are absolutely continuous, because we need absolute continuity of only one of the source and target measures in order to get $L^1$ estimates for the transport density. Therefore, one can also apply the same reasoning in the case of the weighted least gradient problem; we use the same decomposition of $f^\pm$ into atomic and absolutely continuous parts and the fact that in Theorem 6.3 one only needs absolute continuity of the source (or the target) measure. □

**Acknowledgements.** WG acknowledges support of the FWF grant I4354, the OeAD-WTZ project CZ 01/2021, and the grant 2017/27/N/ST1/02418 funded by the National Science Centre, Poland.

## REFERENCES

[1] M. Amar and G. Bellettini, *A notion of total variation depending on a metric with discontinuous coefficients*, Ann. Inst. H. Poincaré Anal. Non Linéaire **11** (1994), 91–133.

[2] L. Ambrosio, N. Fusco, and D. Pallara, *Functions of bounded variation and free-discontinuity problems*, Oxford Math. Monogr., Oxford, 2000.

[3] G. Anzellotti, *Pairings between measures and bounded functions and compensated compactness*, Ann. di Matematica Pura ed Appl. IV **135** (1983), 293–318.

[4] M. Beckmann, *A continuous model of transportation*, Econometrica **20** (1952), 643–660.

[5] E. Bombieri, E. de Giorgi, and E. Giusti, *Minimal cones and the Bernstein problem*, Invent. Math. **7** (1969), 243–268.

[6] P. Cannarsa and C. Sinestrari, *Semiconcave Functions, Hamilton—Jacobi Equations, and Optimal Control*, Nonlinear Differential Equations and Their Applications, Birkhäuser, Basel, 2004.

[7] G-Q. Chen and H. Frid, *Divergence-measure fields and hyperbolic conservation laws*, Arch. Rational Mech. Anal. **147** (1999), 89–118.

[8] G. Crasta and V. De Cicco, *Anzellotti’s pairing theory and the Gauss-Green theorem*, Adv. Math. **343** (2019), 935–970.

[9] L. De Pascale, L.C. Evans, and A. Pratelli, *Integral estimates for transport densities*, Bull. of the London Math. Soc. **36** (2004), no. 3, 383–395.

[10] L. De Pascale and A. Pratelli, *Regularity properties for monge transport density and for solutions of some shape optimization problem*, Calc. Var. Par. Diff. Eq. **14** (2002), no. 3, 249–274.

[11] __________, *Sharp summability for monge transport density via interpolation*, ESAIM Control Optim. Calc. Var. **10** (2004), no. 4, 549–552.

[12] S. Dweik, *Weighted Beckmann problem with boundary costs*, Quarterly of applied mathematics **76** (2018), 601–609.

[13] __________, *Transport and control problems with boundary costs: regularity and summability of optimal and equilibrium densities*, thesis, 2018.

[14] S. Dweik and W. Górny, *Least gradient problem on annuli*, Analysis & PDE, to appear.

[15] S. Dweik and F. Santambrogio, *$L^p$ bounds for boundary-to-boundary transport densities, and $W^{1,p}$ bounds for the BV least gradient problem in $2D$, Calc. Var. Partial Differential Equations **58** (2019), no. 1, 31.

[16] M. Feldman and R. McCann, *Monge’s transport problem on a Riemannian manifold*, Trans. Amer. Math. Soc. **354** (2002), 1667–1997.

[17] W. Górny, *Planar least gradient problem: existence, regularity and anisotropic case*, Calc. Var. Partial Differential Equations **57** (2018), no. 4, 98.

[18] __________, (*Non*)uniqueness of minimizers in the least gradient problem, J. Math. Anal. Appl. **468** (2018), 913–938.

[19] __________, *Existence of minimisers in the least gradient problem for general boundary data*, Indiana Univ. Math. J. **70** (2021), no. 3, 1003–1037.

[20] __________, *Applications of optimal transport methods in the least gradient problem*, preprint, available at arXiv:2102.05887 (2021).

[21] W. Górny, P. Rybka, and A. Sabra, *Special cases of the planar least gradient problem*, Nonlinear Anal. **151** (2017), 66–95.

[22] D. Hauer and J.M. Mazón, *The Dirichlet-to-Neumann operator associated with the 1-Laplacian and evolution problems*, Calc. Var. Partial Differential Equations, to appear.

[23] R.L. Jerrard, A. Moradifam, and A.I. Nachman, *Existence and uniqueness of minimizers of general least gradient problems*, J. Reine Angew. Math. **734** (2018), 71–97.

[24] L. Kantorovich, *On the transfer of masses*, Dokl. Acad. Nauk. USSR **37** (1942), 7–8.

[25] P. Lahti, L. Malý, N. Shanmugalingam, and G. Speight, *Domains in metric measure spaces with boundary of positive mean curvature, and the Dirichlet problem for functions of least gradient*, J. Geom. Anal. **29** (2019), no. 4, 3176–3220.

[26] F. Maggi, *Sets of finite perimeter and geometric variational problems: an introduction to Geometric Measure Theory*, Cambridge Studies in Advanced Mathematics, Cambridge University Press, Cambridge, 2012.
[27] J.M. Mazón, *The Euler-Lagrange equation for the anisotropic least gradient problem*, Nonlinear Anal. Real World Appl. **31** (2016), 452-472.

[28] J.M. Mazón, J.D. Rossi, and S. Segura de León, *Functions of least gradient and 1-harmonic functions*, Indiana Univ. Math. J. **63** (2014), 1067-1084.

[29] V. Maz'ya and S. Poborchi, *Differentiable functions on bad domains*, World Scientific Publishing Co., Inc., River Edge, NJ, 1997.

[30] G. Monge, *Mémoire sur la théorie des déblais et des remblais*, Histoire de l'Académie Royale des Sciences de Paris (1781), 666-704.

[31] A. Moradifam, A.I. Nachman, and A. Tamasan, *Uniqueness of weighted least gradient problems arising in conductivity imaging*, Calc. Var. Partial Differential Equations **57** (2018), no. 1, 6.

[32] A. Moradifam, *Existence and structure of minimizers of least gradient problems*, Indiana Univ. Math. J. **67** (2018), no. 3, 1025-1037.

[33] A. Pratelli, *Equivalence between some definitions for the optimal mass transport problem and for the transport density on manifolds*, Ann. Mat. Pura Appl. **184** (2005), no. 2, 215-238.

[34] P. Rybka and A. Sabra, *The planar Least Gradient problem in convex domains, the case of continuous datum*, Nonlinear Anal. **214** (2022), 112595.

[35] F. Santambrogio, *Absolute continuity and summability of transport densities: simpler proofs and new estimates*, Calc. Var. Partial Differential Equations **36** (2009), 343-354.

[36] __________, *Optimal transport for applied mathematicians*, Progress in Nonlinear Differential Equations and Their Applications 87, Birkhäuser, Basel, 2015.

[37] G. Spradlin and A. Tamasan, *Not all traces on the circle come from functions of least gradient in the disk*, Indiana Univ. Math. J. **63** (2014), 1819-1837.

[38] P. Sternberg, G. Williams, and W.P. Ziemer, *Existence, uniqueness, and regularity for functions of least gradient*, J. Reine Angew. Math. **430** (1992), 35-60.

[39] C. Villani, *Topics in optimal transportation*, American Mathematical Society, Graduate Studies in Mathematics Vol. 58, 2003.

[40] A. Zuniga, *Continuity of minimizers to the weighted least gradient problems*, Nonlinear Analysis **178** (2019), 86-109.

S. DWEIK: LABORATOIRE DE MATHÉMATIQUES DE VERSAILLES, UVSQ, UNIVERSITÉ PARIS-SACLAY, 45 AVENUE DES ETATS-UNIS, F-78035 VERSAILLES CEDEX.

W. GÓRNY: FACULTY OF MATHEMATICS, UNIVERSITY OF VIENNA, OSKAR-MORGERNSTERN-PLATZ 1, 1090 WIEN; FACULTY OF MATHEMATICS, INFORMATICS AND MECHANICS, UNIVERSITY OF WARSAW, BANACHA 2, 02-097 WARSAW, POLAND.

Email address: samer.dweik@uvsq.fr, wojciech.gorny@univie.ac.at