A Martin

M E Negretti

M Ungarish

T Zemach

Propagation of a continuously supplied gravity current head down bottom slopes

We present a combined theoretical-experimental investigation of the downslope propagation of a gravity current sustained by a source. The current propagates first on a horizontal bottom, then on a downslope. We focus on the case when the current at the ridge (point where donwslope begins) has a stable interface (Ri > 0.25) and is critical with F = 1, where Ri and F are the bulk Richardson and flow Froude numbers. We derive the equations that govern the nose propagation and speed using a shallow-water (SW) model, in which the nose is a jump matched to characteristics emitted at the ridge. This provides a self-contained prediction for the speed of propagation u N and position ξ N of the nose. The predicted u N increases with time and distance ξ from the ridge. Since Ri decreases with ξ in the tail behind the nose, appearance of instabilities at a certain traveled distance determines the domain of validity of the SW solution. A good agreement is reported with various experiments with different initial conditions at the ridge and slope angles (both fixed and changing with distance from the ridge). It is shown that the nose velocity is always less than the maximum velocity within the current head, which corresponds to the speed of the characteristics released at the ridge that catch on the current head.

I. INTRODUCTION

There is a long-standing interest in gravity driven flows since they are ubiquitous in natural and man-made environments. Examples include snow avalanches [START_REF] Hopfinger | Snow avalanche motion and related phenomena[END_REF], dust or pyrocastic flows [START_REF] Huppert | The slumping of gravity currents[END_REF], and bottom water formation in the oceans [START_REF] Sherwin | Mixing and advection of a cold water cascade over the Wyville Thomson ridge[END_REF].

The typical system of gravity current, which has received much attention in experimental, theoretical, and numerical investigations, is of lock-release of a dense fluid over a horizontal bottom in a tank filled with a slightly less dense fluid. Once the barrier is removed, the heavier fluid forms a well-defined frontal region after an initial accelerating transient regime. The Reynolds number is large, and the motion is governed by a buoyancy-inertia balance. After a short accelerating transient regime, the current propagates with a constant velocity u N for a period referred to as the slumping stage; then it enters into a self-similar stage with decaying u N and eventually the viscous forces begin to dominate inertia. These features are well documented by experimental investigations 1 FIG. 1. Sketch of the experimental setup and definition of notations. The horizontal component of the velocity is U and the along-slope component of the velocity is u. All of the notations are listed in Table I.

[e.g., [START_REF] Huppert | The slumping of gravity currents[END_REF][START_REF] Rottman | Gravity currents produced by instantaneous releases of a heavy fluid in a rectangular channel[END_REF] and are amenable to interpretation and prediction by theoretical models of the shallowwater (SW) type (see [START_REF] Ungarish | An Introduction to Gravity Currents and Intrusions[END_REF] and references herein).

An important property of the horizontal gravity current is the stability of the interface: the bulk Richardson number Ri = (g ′ h)/u 2 is fairly large for a long time and distance of propagation (see [START_REF] Johnson | Entraining gravity currents[END_REF]). Herein, g ′ = g ρ c -ρ a ρ a is the reduced gravity, where ρ c and ρ a are respectively the density of the gravity current and the ambient, h is the depth of the current, and u the depth-averaged speed of the current, as shown in Fig. 1. In such a configuration, the shallow water model combined with the characteristics method and a jump condition can be used to predict the nose behavior [START_REF] Klemp | On the dynamics of gravity currents in a channel[END_REF].

Several previous studies considered a continuous supply of the gravity current with the constant flux B 0 = g ′ q, where q is the initial flow rate per unit of width imposed by a source at the top of the slope (Britter and Linden [START_REF] Britter | The motion of the front of a gravity current travelling down an incline[END_REF] for the head; Ellison and Turner [START_REF] Ellison | Turbulent entrainment in stratified flows[END_REF] and Pawlak and Armi [START_REF] Pawlak | Mixing and entrainment in developing stratified currents[END_REF] for the tail). In these studies, the current emitted at the source is in the Ri < Ri c ≈ 0.25 regime at the ridge and Kelvin-Helmholtz vortices; entrainment and drag are present from the beginning. We can expect that the speed in the steady-state tail is, from the beginning, of the order of magnitude (B 0 /Ri c ) 1/3 ≈ 1.6B 1/3 0 (see Fig. 9 in Ellison and Turner [START_REF] Ellison | Turbulent entrainment in stratified flows[END_REF]). Britter and Linden [START_REF] Britter | The motion of the front of a gravity current travelling down an incline[END_REF] derived the semiempirical approximation for the (constant) nose propagation speed u N = 1.5B 1/3 0 . This result is in fair agreement with experiments (within about ±15%) for a wide range of slope angles. The flow is affected by entrainment and drag and, due to the interfacial instabilities, the initial conditions on the characteristics of a SW model are uncertain and need to be verified.

In Negretti et al. [START_REF] Negretti | Development of gravity currents on rapidly changing slopes[END_REF] and Martin et al. [START_REF] Martin | Development of gravity currents on slopes under different interfacial instability conditions[END_REF] (N17 and M19) the source was positioned at some upstream position (see Fig. 1), producing different initial conditions than those usually used in the literature [e.g., [START_REF] Pawlak | Mixing and entrainment in developing stratified currents[END_REF][START_REF] Turner | Buoyancy Effects in Fluids[END_REF]. In N17 and M19 the tail of the gravity current starts on the slope without any entrainment due to an initial stable interface of the current whereas. in the previous literature, the tail of the gravity current starts on the slope with a large entrainment due to an initial unstable interface of the current. Negretti et al. [START_REF] Negretti | Development of gravity currents on rapidly changing slopes[END_REF] demonstrate that beyond the straightforward separation into "formation" and "steady-state" stages, additional subdivisions for the second "steady-state" stage are required. They show that for a gravity flow released far upstream of the edge, the rear part of the gravity flow on the slope has a stable interface, while more downslope the flow displays an increasing speed and decreasing Ri, until a point where the interface becomes unstable. Negretti et al. [START_REF] Negretti | Development of gravity currents on rapidly changing slopes[END_REF] have shown for such a configuration that the motion of the nose can be well approximated by a "free-fall" behavior: assuming that the head is a rigid body, the potential energy in the reduced-gravity field is converted into kinetic energy as the nose advances from the initial height

Z i to Z N (t ) = Z i -tan θ [X N (t ) -X i ],
where (X, Z ) are the horizontal and vertical coordinates, θ the slope angle, as depicted in Fig. 1, and t denotes time. This consideration yields the formula

u N = 2S 2 g ′ (Z i -Z N ) + u 2 Ni 1/2 , (1)
where S 2 is a shape factor related to the shape of the velocity and density vertical profiles [START_REF] Turner | Turbulent entrainment: The development of the entrainment assumption, and its application to geophysical flows[END_REF] and can be interpreted as an adjustable constant. The advantage of this approximation is that it reproduces surprisingly well the maximal free stream velocity in the head [START_REF] Negretti | Development of gravity currents on rapidly changing slopes[END_REF]. The deficiency is that this is actually a curve fit which lacks theoretical justification, needs empirical inputs, and does not provide insight about the mechanism that governs the propagation of the nose. Here, we consider the complex system of Negretti et al. [START_REF] Negretti | Development of gravity currents on rapidly changing slopes[END_REF] and Martin et al. [START_REF] Martin | Development of gravity currents on slopes under different interfacial instability conditions[END_REF] in which the current is produced and sustained by a continuous supply of dense fluid (constant flux B 0 ) and the propagation is over an inclined bottom. In this case, the flow has two main stages: the time-dependent formation of a stream (tail) following the nose moving with u N , then a possible steady state along the entire slope or part of it. Intuitively, the high-Reynolds-number stream is expected to accelerate and contract during propagation. This suggests that time dependency and interface stability effects produce a different flow pattern from that of the standard system.

When the current emitted by the source has a stable interface, the standard shallow-water (or hydraulic) equations for the thickness h(ξ,t ) and depth-averaged u(ξ,t ) can be used, where ξ denotes the along-slope coordinate (see Fig. 1). The stable interface is a necessary condition for the propagation of the SW characteristics to the nose. The analysis of the characteristics indicates that at the ridge the flow is "critical" (the Froude number F = Ri -2 of the flow is equal to 1) with u i = (g ′ h i ) 1/2 , where the subscript i denotes the ridge position. Since u i h i = q, we obtain two useful results: (1) from the definition of B 0 , it follows that u i = B 1/3 0 = U where U is the characteristic velocity scale; and (2) from the definition of Ri it follows that Ri i = 1, which justifies the assumption of stable interface at the ridge (in agreement with experimental observations). In this case, for a while, the nose propagates in stable conditions and leaves behind a smooth steady-state tail in which u(ξ ) increases while h(ξ ) and Ri decreases. After the position where Ri is sufficiently small for instabilities to develop, the gravity current is affected by entrainment and the simple SW equations become invalid. Our analysis is focused only on the stable stage, and hence the SW equations are a good approximation.

The propagation of the head with stable initial conditions lacks theoretical understanding and means of prediction.

This paper attempts to close this gap of knowledge concerning the propagation of the nose of the current, ξ = ξ N (t ), and corresponding speed u N (t ) for a continuously supplied gravity current with a stable interface at the ridge during the "formation" phase by use of the SW framework for gravity currents.

The SW equations are a hyperbolic system, and there is good evidence from solutions of horizontal currents that the nose is a discontinuity where a jump condition can be applied [START_REF] Ungarish | An Introduction to Gravity Currents and Intrusions[END_REF]. We implement this idea in the present problem, and show that self-contained results for u N (t ) can be obtained. To corroborate the theoretical prediction and sharpen the understanding we compare with laboratory experimental velocity data collected using a particle image velocimetry (PIV) measurement technique.

The paper is organized as follows. In Sec. II we present the SW formulation and derive the theoretical prediction of u N . The experimental setup, and methodology for the detection of the propagation of the nose are described in Sec. III. The comparison between theory and measurements is performed in Sec. IV and discussed in Sec. V. Finally, some concluding remarks are given in Sec. VI. All of the symbols that are introduced for this study are listed in Table I to improve the readability of the paper.

II. SHALLOW WATER MODEL

We consider a layer of fluid of density ρ c which propagates over a bottom into an ambient fluid of density ρ a driven by the reduced gravity g ′ . The flow is two-dimensional. The coordinates (ξ,ζ) Nose vertical coordinate are along and normal to the sloping bottom, as depicted in Fig. 1. In this system, u is the velocity of the denser layer and h its depth (thickness). The SW model involves the following assumptions: the flow is in the inertial regime (large Reynolds number), and the horizontal length scales are much larger compared to the vertical scales. Vertical accelerations are small so that the pressure is well approximated by the hydrostatic balance. Also, we assume a Boussinesq current and a deep ambient so that the return flow above the current is negligible. Assuming a stable interface implies negligible entrainment, which justifies the SW approach. The governing equations (details of derivation can be found in [START_REF] Ungarish | An Introduction to Gravity Currents and Intrusions[END_REF]) give the volume, continuity, and momentum balance:

∂h ∂t + ∂hu ∂ξ = 0, (2)
∂u ∂t + u ∂u ∂ξ = g ′ - ∂h ∂ξ 1 cos θ + sin θ . ( 3
)
The system is hyperbolic. The characteristics are

du ± g ′ h 1/2 dh = g ′ tan θdt on dξ dt = c ± = u ± (g ′ h) 1/2 . ( 4
)
The initial condition for the gravity flow is given by a step function of the flux q at t = 0. Inspection of the characteristics indicates that at the ridge, denoted by the subscript i,t h e conditions are critical, so that c -= 0 and u i = √ g ′ h i . As for currents over horizontal boundaries, also the gravity currents over sloping boundaries require a jump condition at the nose that can be derived using force balances in a control volume (see [START_REF] Ungarish | An Introduction to Gravity Currents and Intrusions[END_REF][START_REF] Ungarish | Models of internal jumps and the fronts of gravity currents: Unifying twolayer theories and deriving new results[END_REF]) whose boundaries are separated by a small δx = δξ cos θ. By reviewing the analysis, we conclude that for moderate angles (θ<45 • say) the inclination of the bottom has negligible contribution to the control-volume balances about a thin jump. Therefore, we apply the standard jump conditions on the horizontal and on the slope:

U = F g ′ H and u = F (cos θ ) 3/2 g ′ h [ξ = ξ N (t )] (5)
respectively, where U and H are respectively the velocity and height of the current on the horizontal (X, Z). F = F (a) is the Froude number given by an "off-the-shelf" formula [START_REF] Huppert | The slumping of gravity currents[END_REF][START_REF] Benjamin | Gravity currents and related phenomena[END_REF] expressed in terms of the vertical thickness ratio of the current to the ambient a = h/(H cos θ ), with H the total water column thickness.

Benjamin [START_REF] Benjamin | Gravity currents and related phenomena[END_REF] proposed a relation for F (a) based on analytical considerations:

F (a) = (2 -a)(1 -a) 1 + a 1/2 ; ( 6 )
and Huppert and Simpson [START_REF] Huppert | The slumping of gravity currents[END_REF] proposed a relation for F (a) that relies on experimental data:

F (a) = 1 2 a -1/3
for a 0.075 and F (a) = 1.19 for a < 0.075.

The motion of the nose of the SW current is governed by two effects: (1) the information at the ridge ξ = 0: u = u i , h = h i ; and (2) the jump condition as given in Eq. ( 5), namely u = u N = φ(g ′ h) 1/2 at ξ N (t ), where φ = F/(cos θ ) 3/2 . The connection between these conditions is provided by a characteristic whose speed is c + = u + √ g ′ h and on which du = g ′ tan θ dt -√ g ′ dh/ √ h.The trajectories are sketched in Fig. 2. For concave slopes, the exact integration along the characteristic is cumbersome, and we take an averaged angle θ assuming that the variations along the characteristic are modest. Supposing that the characteristic is released at the time t k > 0(k ∈ N) and intersects the nose at time t, the integral along the characteristic k [cf. Eq. ( 4)] gives

u N = u i + g ′ tan θ (t -t k ) -2 g ′ ( h N -h i ). (8)
We recall that the subscript i denotes the ridge. The of the of the nose with the characteristic yields

1 2 (u N + u Ni )t = 1 2 [u i + u N + g ′ ( h i + h N )](t -t k ), (9)
where u Ni is the value of u N at t = 0 obtained from q = u Ni h Ni . Recalling that u Ni = φ √ g ′ h Ni and that u i = (g ′ h i ) 1/2 we rewrite the Eqs. ( 8) and (9)as

u N = 1 + 2 φ -1 [3u i + g ′ tan θ (t -t k )], (10)
t -t k = u N + u Ni 2u i + 1 + 1 φ u N . (11)
Eliminating (t -t k ) from Eqs.. ( 10) and [START_REF] Negretti | Development of gravity currents on rapidly changing slopes[END_REF], we obtain

u N = 1 + 2 φ -1 ⎡ ⎣ 3u i + g ′ u N + u Ni 2u i + 1 + 1 φ u N tan θ t ⎤ ⎦ , ( 12
)
where θ is the mean angle from the ridge position to the nose position approximated by θ = (θ 0 + θ )/2. The variable φ = F/ cos 3/2 θ is evaluated using both relations ( 6) and [START_REF] Klemp | On the dynamics of gravity currents in a channel[END_REF]. Combining one of these relations for F (a) with the condition of the Froude number, we obtain

u N = F (a) (cos θ ) 3/2 g ′ h N = φ(a) g ′ h N . (13)
Equation ( 12) can be numerically solved starting from the initial values and using a time advancement of δt, here chosen equal to 0.1 after testing the accuracy with smaller values.

III. EXPERIMENTAL PROCEDURES

The experiment, sketched in Fig. 1, and described in detail elsewhere [START_REF] Negretti | Development of gravity currents on rapidly changing slopes[END_REF], consists of a uniformdensity saline solution injected continuously from an external reservoir on a 25 cm wide and 20 cm TABLE II. Summary of the experimental parameters of the PIV velocity measurements. The total initial water depth in the channel is H 0 = 19 cm, the current release distances are D1 = 50 cm and D2 = 160 or 200 cm. X is the horizontal spatial coordinate defined in Fig. 1 and R is the curvature radius also defined in Fig. 1 II. A PIV measurement technique was used to obtain two-dimensional velocity fields in a field of view of M × N, where M and N are respectively 70 and 50 cm for experiments 2 and 3 and 86 and 65 cm for the other experiments (cf. Fig. 1). Further details are given in Negretti et al. [START_REF] Negretti | Development of gravity currents on rapidly changing slopes[END_REF] and Martin et al. [START_REF] Martin | Development of gravity currents on slopes under different interfacial instability conditions[END_REF]. The use of the PIV technique permits to have detailed information about the velocity distribution within the head and give some more options for defining the nose limits and the front speed [START_REF] Thomas | The structure of the head of an inertial gravity current determined by particle-tracking velocimetry[END_REF][START_REF] Nogueira | Dynamics of the head of gravity currents[END_REF]. Figure 3(a) gives an instantaneous image of the velocity field for experiment 6; the vertical gray line represents the position ξ (t ) at which the maximum velocity at the time t = 9 s is below the value of the reference velocity U . The velocity U is hence used as a threshold on the maximum velocity within the current at the time t to detect the nose position [vertical gray line in Fig. 3(b)], which is very close to the position of the nose calculated using the maximum spatial velocity gradient at each time t [red curve in Fig. 3(b)] but exhibiting larger fluctuations. Furthermore, Fig. 3(b) highlights that the maximum of the velocity within the head is upstream of the nose, as also reported by Thomas et al. [START_REF] Thomas | The structure of the head of an inertial gravity current determined by particle-tracking velocimetry[END_REF]. Figure 3

∂ u / ∂ x (s -1 ) (c)
FIG. 3. (a) Instantaneous (t = 9 s) along-slope velocity for a linear slope of 10 • (experiment 6), u = U (black contour). (b) Instantaneous free stream velocity with ξ (black) and spatial derivative of the free stream velocity (red). (c) Hovmöller diagram of the free stream velocity with ξ and t, time position of the snapshot (horizontal black line). The nose position is a grey line. The green cross represents the maximum velocity in the head as reported in [START_REF] Negretti | Development of gravity currents on rapidly changing slopes[END_REF]. . Model solution (continuous lines) and experiments (symbols) obtained from a threshold on the free stream velocity (blue triangles) and from the maximum spatial derivative of the free stream velocity (cyan triangles). Theoretical ũN and ξN are computed with F (a), which is estimated with the jump condition (6) (black line) and with Eq. ( 7) (red line). In (c) the theoretical ũN is estimated by using a constant averaged angle θ. function of space and time. The gray line gives the so obtained front advancement ξ N in time t.The nose propagation speed is defined as u N = dξ N /dt.

The values of ξ N are binned using a running window of dimensionless time width equal to 1. The bins are differentiated over a t = 1 such that

u N ( t ) = t+1 t ξN ( t′ )d t′ - t t-1 ξN ( t′ )d t′ , (14)
where variables are made dimensionless using L = (q 2 /g ′ ) 1/3 , T = (q/g ′2 ) 1/3 , and U = B 1/3 0 = (qg ′ ) 1/3 and are denoted by ˜. We highly smoothed the nose velocity u N with a running window of t = 7.

The rms of the difference between raw velocities and the smoothed velocities is smaller than 0.2cm/s.

Negretti et al. [START_REF] Negretti | Development of gravity currents on rapidly changing slopes[END_REF] defined the head velocity as the maximum PIV velocity within the head at each time t (cf. the cross in Fig. 3) and used it to compare with their free-fall velocity solution. This is highlighted by a cross in Figs. 3(a)-3(c). We point out that this is an alternative definition that is different from the Lagrangian velocity u N = dξ N /dt (see also Thomas et al. [17]). the D1 configuration (intermediate initial stable condition) is faster than that in the D2 configuration (initial stable condition). This difference may be ascribed to a higher dilution of the head for the D2 experiment due to the higher release distance.

IV. COMPARISON WITH EXPERIMENTS

The dimensionless nose velocity ũN = d ξN /d t is displayed in Fig. 4 (bottom panels). We see that ũN increases with t, with a rate of less than 20% for the experiments with 7 • and of 30% for the experiments on a larger slope angle.

In Fig. 4, for fixed-slope configurations, the theoretical nose propagation (top panels) and speed (bottom panels) are plotted as continuous lines using the two different relations for F (a): black lines represent Eq. ( 6) and red lines Eq. [START_REF] Klemp | On the dynamics of gravity currents in a channel[END_REF]. The nose displacement and the nose velocity reveal a strong dependency on the Froude condition with an impact of up to 10% for the displacement and up to 25% for the velocity.

The theoretical predictions reveal a continuously accelerating current in full accord with the experimental data.

In Fig. 5, for a changing-slope configuration, the theoretical nose propagation (a) and speed (b) are plotted as in Fig. 4. For experiment 4 the travel time of the nose is shorter than for the experiments of Fig. 4 because the field of view of the camera downstream of the ridge is smaller, the characteristic velocity of the nose is higher (due to a higher pump flow rate), and the slope is steeper. The theoretical solution shows a good agreement with the observations. The comparison of Figs. 5(b) and 5(c) shows that both the theoretical prediction for experiment 4 with a changing slope and the prediction with an average constant slope give a correct estimate of the nose displacement and velocity. Nevertheless, the theoretical prediction using a variable angle is much better than the prediction with an average constant slope, particularly for the derivative of the nose velocity. This demonstrates that the theoretical model is able to capture well the effects of the slope.

V. DISCUSSION

The presented theoretical SW model of the governing equations of motion combined with a nose-jump (of Benjamin type) condition can predict reasonably well the front displacement and velocity for most of the experiments with both constant and changing slope, paricularly in the first stages, as long as the effects of entrainment and associated drag are negligible. This is a significant improvement of the understanding and of the prediction tools for the motion of a sustained gravity current on a downslope of constant or changing-with-distance slope.

As pointed out in Sec. II, a necessary condition for validity of the SW equations and the corresponding characteristics of the nose is that the tail is not affected by entrainment and drag. number at the ridge and along the slope. We showed that for certain initial conditions the stability of the interface can be maintained for a significant distance of propagation.

We demonstrate that, in the case of a stable interface at the ridge, u N is a match between a jump condition at the moving nose and information carried by the characteristic from the ridge through the stable tail. This u N increases with time (and downslope distance), Ri decreases with time, and hence at some position the stable u N switches to the entrainment-drag dominated value 1.5B 1/3 0 . It is difficult to assess the accuracy of the theoretical prediction. First, there is some uncertainty concerning the Froude jump condition F = F (a). The semiempirical Huppert and Simpson [START_REF] Huppert | The slumping of gravity currents[END_REF] formula gives better agreement than the more rigorous Benjamin formula; this, however, is consistent with the behavior of gravity currents on a horizontal bottom. Second, and more important, are the limitations of the experimental data. The position ξ N of the nose is not a sharp experimental variable, and the calculation of the time derivative u N yields noisy data that need strong smoothing. We employed several methods of measurement and smoothing, but some discrepancies remain, which could be improved for example using high-resolution Navier-Stokes simulations.

The Lagrangian velocity of propagation of the nose, u N = dξ/dt, is significantly smaller than the interior maximal velocity in the head. This difference may be important in geophysical and environmental applications, where both position of the current and internal convection are of interest. The details of the transition from the SW regime to the entrainment-drag regime require a complex investigation, which is left for future work.

FIG. 2 .

 2 FIG.2. Sketch of the nose propagation (black line) and of the characteristics released at t 1 , t 2 , t 3 (blue dashed lines).

  (c) shows a Hovmöller diagram of the free stream velocity which is a

FIG. 4 .

 4 FIG. 4. Nose position (top panels) and nose velocity (bottom panels), for the configurations D1 [(a)-(c) and (e)-(g)] and D2 [(d) and (h)]. Model solution (continuous lines) and experiments (symbols) obtained from a threshold on the free stream velocity (blue triangles) and from the maximum spatial derivative of the free stream velocity (cyan triangles). Theoretical ũN and ξN are computed with F (a), which is estimated with the jump condition (6) (black line) and with Eq. (7) (red line). In (c) the theoretical ũN is estimated by using a constant averaged angle θ.

Figure 4 (FIG. 5 .

 45 Figure 4 (top panels) shows the dimensionless nose position ξN from the ridge against the dimensionless time t for experiments 8, 9, and 10. We see that the nose position ξ increases practically linearly with time. No significant change is reported when the slope angle changes during the first time period of 5T . The slightly different patterns of the nose trajectories observed between the experiments are induced by different initial stability of the tail at the ridge [see D1 configuration in Fig. 4(a) and D2 configuration in Fig. 4(d) for a linear slope of 15 • ]. The nose displacement for

TABLE I .

 I Notations.

	Symbol	Definition
	a	Vertical thickness ratio of the current to the ambient
	α	Coefficient of proportionality between nose velocity and buoyancy flux
	B 0	Buoyancy flux
	c ± D	Characteristics velocity Current release distance between the gate and the ridge
	F	Froude number
	g	Acceleration of gravity on the Earth's surface
	g ′	Reduced gravity of the gravity current
	g ′ e	Effective reduced gravity of the gravity current
	ζ	Normal coordinate to the sloping bottom
	H	Height of the gravity current on the vertical coordinate
	H	Water column thickness
	H 0	Water column thickness upstream of the ridge
	h	Height of the gravity current on the normal coordinate to the sloping bottom
	h 0	Current thickness at the gate
	h i	Height of the gravity current at the ridge
	h N	Height of the gravity current at the nose position
	h Ni	h N at the ridge
	θ	Slope angle
	θ	Mean slope angle from the ridge to the nose position
	θ end	Slope angle at the end of the slope
	L	Characteristic length scale
	M	Widthoffieldofview
	N	Height of field of view
	ξ	Along-slope coordinate
	ξ max	Length of the slope
	ξ N	Nose position along position
	q	Flow rate per unit of width
	q i	Flow rate per unit of width at the ridge
	q 0	Flow rate per unit of width at the gate
	R	Curvature radius of the concave slope
	Ri	Bulk Richardson number
	Ri c	Critical bulk Richardson number
	ρ a	Density of the ambient fluid
	ρ c	Density of the gravity current
	S 2	Shape factor in Eq. (1)
	T	Characteristic timescale
	t	Time
	U	Depth-averaged velocity of the gravity current on the horizontal coordinate
	U	Characteristic velocity scale
	u	Depth-averaged velocity of the gravity current on the along-slope coordinate
	u i	Depth-averaged velocity of the gravity current at the ridge
	u N	Nose velocity
	u Ni	Nose velocity at the ridge
	φ	Froude number in the slope coordinate
	X	Horizontal coordinate
	X N	Nose horizontal coordinate
	X i	X at the ridge
	Z	Vertical coordinate
	Z i	Z at the ridge
	Z N

  . θ • end is the slope angle at the end of the slope. a pump through a gate opening of h 0 = 6 cm. The total length of the horizontal channel is 230 cm and at its end different inclined boundaries (linear, concave) were employed to enable the flow to accelerate. Further variation parameters were the injected buoyancy flux, the gate position D on the horizontal portion [160 or 200 cm (D2) and 50 cm (D1) upstream of the ridge] that induce different initial stability conditions (interfacial Richardson numbers) at the ridge. All the parameters of the experiments are reported in Table

	Expt.	1	2	3	4	5	6	7	8	9	10
	θ • θ • end	15 15	22 22	17.2 -X/R 0	31.5 -X/R 0	15 5	10 10	7 7	15 15	10 10	7 7
	D (cm)	160		200			160			50
	q (cm 2 s -1 )3		2 4			1 6			1 6
	g ′ (cm 2 s -1	5.6		4.3			4.3			4.3
	B 0 (cm 3 s -3 )	179		103			69			69
	deep channel via

ACKNOWLEDGMENTS

This work has been partially supported by the LabEx Tec 21 (Investissements d'Avenir, Grant Agreement No. ANR-11-LABX-0030). We are grateful to P. Jarrin, S. Mercier, and S. Viboud for their technical support for the experiments. Special thanks go to E. J. Hopfinger for useful discussions.

Consequently, the relevant experiments are those for gravity currents with stable initial conditions at the ridge.

The SW equations ( 2) and (3) admit a steady-state solution, which corresponds to the tail left behind the nose. Letting u = u(ξ ), h = h(ξ ) = q/u(ξ ), and using the boundary conditions at ξ i ,we obtain the SW solution

where we set a constant θ and cos θ ≈ 1 for simplicity, without any loss of generality. The equation for u is implicit, but the behavior is clear, as follows. Defining ξ = 0 at the ridge we conclude that in the tail u increases monotonically with ξ ,fromthev alueB 1/3 0 . Formally, the value of u is bounded only by the length of the slope, ξ max . However, we must consider the stability of the interface of this tail. The local Ri(ξ ) = g ′ q/u 3 = B 0 /u 3 decreases with ξ , and hence the interface is expected to become unstable at some ξ c where Ri = Ri c < 0.25. Using [START_REF] Ungarish | Models of internal jumps and the fronts of gravity currents: Unifying twolayer theories and deriving new results[END_REF] we obtain the estimate for the length of the stable tail:

We conclude that the stability criterion limits the length of the SW smooth tail, and also the maximum value of u,t o( B 0 /Ri c ) 1/3 . With a safety margin, based on the results of Negretti et al. [START_REF] Negretti | Development of gravity currents on rapidly changing slopes[END_REF], we infer that the speed of the current in the stable tail is in the range [1, 2]B 1/3 0 and, since the instability produces entrainment and drag, we argue that the speed of the current in the unstable Ri < Ri c domain will not exceed 2B 1/3 0 .FromEq.( 16) the typical scaled ξ and t of the propagation in the stable domain are expected to be of the order of 5/ sin θ. For a slope of 15 • ,5/ sin θ = 19.

The experiments of Negretti et al. [START_REF] Negretti | Development of gravity currents on rapidly changing slopes[END_REF] confirm these predictions: observed currents starting with Ri ≈ 1a tξ i have a stable interface over a significant distance, but develop suddenly Kelvin-Helmholtz (KH) vortices for ξ>ξ c .Atξ c the experimental values of entrainment and bottom drag increase with ξ . The density of the current is diluted, and the effective reduced gravity g ′ e is smaller than the reduced gravity at the gate, g ′ 0 . They are related by the conservation of the buoyancy flux, g ′ e q i = g ′ 0 q 0 = B 0 , where q i is the flow rate per unit of width at the ridge, q 0 at the gate. If no mixing occurred between the gate and the ridge, then q i = q and g ′ = g ′ e , which is what it has been assumed within this paper. The behavior of the flow is very different from the simple SW results (which assume E = c D = 0); in particular, u decreases and h increases with ξξ c . Some distance after ξ c , the observed Kelvin-Helmholtz billows are so large that a depth-averaged approach becomes inappropriate.

A further case, as mentioned in the Introduction, is that of a current with unstable conditions at the ridge, as considered for example in Britter and Linden [START_REF] Britter | The motion of the front of a gravity current travelling down an incline[END_REF], Ellison and Turner [START_REF] Ellison | Turbulent entrainment in stratified flows[END_REF], and Pawlak and Armi [START_REF] Pawlak | Mixing and entrainment in developing stratified currents[END_REF]. Britter and Linden [START_REF] Britter | The motion of the front of a gravity current travelling down an incline[END_REF] reported that, given a buoyancy flux B 0 , after an initial short acceleration phase the forces are quickly balanced (within 10h i ) and the front enters an established regime where the velocity u N is constant, scaling as u N = α(B 0 ) 1/3 , with α = 1.5 and eventually with α a function of the Froude number at the nose and of the slope angle. Tokyay and García [START_REF] Tokyay | Effect of initial excess density and discharge on constant flux gravity currents propagating on a slope[END_REF] obtain α = 1.11-1.17 in their numerical simulations with slope angles θ<3 • . In our experiments we report α = 1.2-1.3f o r ξ>10. This difference may be attributed to the depth ratio, which is larger in the experiments of Britter and Linden, resulting in a larger Froude number at the nose as compared to our experiments.

VI. CONCLUDING REMARKS

We presented a theoretical model which explains and predicts the motion of the nose of a downslope gravity current sustained by a source, and compared the model with experiments. One important conclusion is that the propagation of the nose is strongly affected by the Richardson


Hamza Ennaji
email: hamza.ennaji@unilim.fr

NOUREDDINE IGBIDA † Van Thanh Nguyen
email: nguyenvanthanh@qnu.edu.vn.

BECKMANN-TYPE PROBLEM FOR DEGENERATE HAMILTON-JACOBI EQUATIONS

The aim of this note is to give a Beckmann-type problem as well as corresponding optimal mass transportation associated with a degenerate Hamilton-Jacobi equation H(x, ∇u) = 0, coupled with non-zero Dirichlet condition u = g on ∂Ω. In this article, Hamiltionian H is continuous in both arguments, coercive and convex in the second, but not enjoying any property of existence of a smooth strict sub-solution. We also provide numerical examples.

Introduction

Hamilton-Jacobi (HJ for short) equation, coupled with a Dirichlet condition, is a first order equation of the type

H(x, ∇u(x)) = 0 in Ω u = g on D ⊂ Ω. (1.1)
It has attracted the attention of mathematicians and physicists since it appears in many fields, in particular, classical mechanics, geometry, seismic analysis, optics, image processing, etc. Our main interest here lies in the study of the connection between these types of equation and the Beckmann problem as well as the Monge-Kantorovich problem in the case where H is continuous in both arguments, coercive and convex in the second, but not enjoying any property of existence of a smooth strict sub-solution. Recall that Beckmann's problem (called by Beckmann a continuous model of transportation [START_REF] Beckmann | A continuous model of transportation[END_REF]) is a divergence PDE-constrained optimization problem.

It is important in the study of transportation activities. As to the Monge-Kantorovich problem, it consists in finding the best way to push forward between two given measures (usually called goods and consumers) related to some cost function. It is a linear optimization problem which appears in the study of optimal transportation and allocation of resources.

1.1. Beckmann's problem and optimal transport. Beckmann's problem appears typically in the study of optimal transference between two distributions of masses. For instance in urban area represented by a bounded domain Ω ⊂ R N , N = 2, it can describe some kind of optimal traffic between two distributions of residents and services. These distributions may be represented by two given nonnegative Radon measures µ 1 and µ 2 , respectively. So, the signed measure µ := µ 2 -µ 1 represents the local measure of excess demand. The consumers traffic is given by a traffic flow field, i.e. a vector field Φ : Ω → R N whose direction indicates the consumers' travel direction and whose modulus |Φ| is the intensity of the traffic. The relationship between the excess demand and the traffic flow is obtained from the equilibrium condition:

-div(Φ) = µ in D (Ω).

This condition describes some kind of equilibrium: the outflow of consumers equals the excess demand in any subregion ω ⊂ Ω, i.e.

∂ω Φ • n ds = µ(ω). In standard Beckmann's problem, the measures µ 1 and µ 2 have equal masses (a balanced condition between residents and services), and the urban area is assumed to be isolated, i.e. no traffic flow should cross the boundary of the city Ω. In other words, Φ is subject to the boundary condition

Φ • n = 0 on ∂Ω,
where n denotes the outward normal vector to the boundary. Assuming the transportation cost per consumer is given by the quantity F (x, Φ(x)), for a given function

F : Ω × R N → R + ,
Beckmann therefore argued that one may define the transportation cost between µ 1 and µ 2 as the infimum of the total cost of the traffic

Ω F (x, Φ(x)) dx.
Recall that this problem is tightly related to the Monge-Kantorovich problem (also called optimal transport) which consists in finding the best way to push forward the measure µ 1 into µ 2 related to a specified cost function. For i = 1, 2, we use the convention of denoting by π i : R N ×R N → R N , i = 1, 2 the projections: π 1 (x, y) := x, π 2 (x, y) := y. Given a nonnegative Radon measure µ on R N × R N , its marginals are measures defined by proj x (µ

) := π 1 #µ, proj y (µ) := π 2 #µ, i.e. proj x (µ)(A) := µ(π -1 1 (A)) = µ(A × R N ) and proj y (µ)(B) = µ(R N × B) for any Borel subsets A, B ⊂ R N . Then, the Monge-Kantorovich (MK) problem consists in finding a nonnegative Radon measure µ on R N × R N such that π 1 #µ = µ 1 , π 2 #µ = µ 2 and minimizing the total work c(x, y) dµ(x, y),
where c is a given ground cost function describing the charge one needs to pay for transporting one unit of µ 1 into a unit of µ 2 . A minimizer µ of the above problem is called an optimal plan. One of basic concepts in the optimal transport theory is the Kantorovich duality that can be restated as follows (see e.g. [START_REF] Villani | Optimal transport: old and new[END_REF]Chapter 5]):

Theorem 1.1. Let c be a lower semicontinuous (l.s.c.) cost function and

µ 1 , µ 2 ∈ M + b (R N ) be nonnegative Radon measures such that µ 1 (R N ) = µ 2 (R N ).
Then (i) The (MK) problem has an optimal plan and the Kantorovich duality holds true, i.e. min c(x, y) dµ :

π 1 #µ = µ 1 , π 2 #µ = µ 2 = sup u dµ 1 + v dµ 2 : (u, v) ∈ S c (µ 1 , µ 2 ) ,
where

S c (µ 1 , µ 2 ) := (u, v) ∈ L 1 µ 1 (R N ) × L 1 µ 2 (R N ) : u(x) + v(y) ≤ c(x, y) ∀x, y ∈ R N . (ii) If c(x, y) ≤ C µ 1 (x) + C µ 2 (y) for some (C µ 1 , C µ 2 ) ∈ L 1 µ 1 × L 1 µ 2
, then the dual problem on the right-hand side (called Kantorovich dual problem) has an optimal solution. (iii) If the cost function c satisfies the triangle inequality and c(x, x) = 0 for any x ∈ R N , then the Kantorovich dual problem can be rewritten by using one variable u only:

sup R N ud(µ 2 -µ 1 ) : u ∈ L 1 µ 1 ∩ L 1 µ 2 and u(y) -u(x) ≤ c(x, y) ∀x, y ∈ R N .
The connection between Beckmann's problem and the Monge-Kantorovich problem is well known, for instance, in the case where c is given by the distance d k associated with the Riemannian metric k, i.e.

d k (x, y) = inf 1 0 k(ϕ(t)) |ϕ (t)| dt : ϕ ∈ Lip([0, 1], Ω), ϕ(0) = x, ϕ(1) = 1 . (1.2)
In this case, since c = d k is a distance, the Kantorovich dual problem reads as (see Theorem 1.1)

max Ω u d(µ 2 -µ 1 ) : u(y) -u(x) ≤ d k (x, y) for all x, y ∈ Ω , (1.3)
or equivalently,

max Ω u d(µ 2 -µ 1 ) : |∇u(x)| ≤ k(x) for all x ∈ Ω . (1.4)
Then, using Fenchel-Rockafellar duality, one leads to another dual problem to (1.4):

inf Ω k(x)|Φ|(x) dx : Φ ∈ L 1 (Ω) N , -div(Φ) = µ 2 -µ 1 in D (R N ) .
This is Beckmann's problem with F (x, p) = k(x)|p|. For the general case of non-degenerate

Finsler metric F , one has (see e.g. [START_REF] Benamou | A numerical solution to Monge's problem with a Finsler distance as cost[END_REF][START_REF] Igbida | Optimal partial mass transportation and obstacle Monge-Kantorovich equation[END_REF][START_REF] Santambrogio | Optimal transport for applied mathematicians[END_REF][START_REF] Igbida | Optimal mass transportation for costs given by Finsler distances via p-Laplacian approximations[END_REF]), with d F (x, y) := inf where d k is given by (1.2). Thanks to (1.3), it is meaningful to consider this problem as some kind of push-forward of µ 1 := L N |Ω the restriction of Lebesgue measure on Ω. The offset is clearly connected to the lack of a measure µ 2 which can fit out the problem with the balanced property. Actually, the linked Beckmann and Monge-Kantorovich problems aim to find moreover the optimal µ 2 concentrated on ∂Ω which will consume µ 1 := χ Ω . More precisely, we will prove that the problem is connected to

ϕ∈Lip([0,1],Ω) ϕ(0)=x, ϕ(1)=1 1 0 F (ϕ(t), ϕ (t))dt, min Ω×Ω d F (x, y) dµ : π 1 #µ = µ 1 , π 2 #µ = µ 2 = max Ω u d(µ 2 -µ 1 ) : u(y) -u(x) ≤ d F (x, y) for all x ∈ Ω = inf Ω F (x, Φ(x)) dx : Φ ∈ L 1 (Ω) N , -div(Φ) = µ 2 -µ 1 in D (R N ) .
min ν∈M b (∂Ω) max u Ω u (dx -dν) - ∂Ω g dν : u is 1 -Lipschitz w.r.t d k
which provides the following modified Beckmann problem

min ν∈M b (∂Ω) inf Φ∈L 1 (Ω) N Ω k(x) |Φ|(x) dx + ∂Ω g dν : -div(Φ) = 1 -ν in D (R N ) .
as well as to the following modified Monge-Kantorovich problem

min γ∈M + (Ω×Ω),ν∈M b (∂Ω) Ω×Ω d k (x, y)dγ(x, y) + ∂Ω g dν : π 1 γ = ν + , π 2 γ = 1 + ν -.
For the general case of degenerate Hamiltonian H, we will also obtain the corresponding results by means of the intrinsic distance d σ associated to the Hamiltonian H (see (2.7)).

The main contribution of our paper is the rigorous treatment for the case of degenerate Hamiltionian H and its degenerate intrinsic metric d σ , as well as non-zero Dirichlet condition. Moreover, we also illustrate numerical examples. This paper is organized as follows. In Section 2 we recall briefly some notions around HJ equation as well as the main results establishing the connection between HJ equation and Beckmanntype problem. Section 3 is dedicated to the proofs of the main results. Numerical examples are presented in Section 4. Finally, Section 5 contains some comments and extensions.

Preliminaries and main results

2.1. Preliminaries. All the properties and notions introduced in this subsection are more or less known. Proofs and more details can be found in [START_REF] Michael | User's guide to viscosity solutions of second order partial differential equations[END_REF][START_REF] Lions | Generalized solutions of Hamilton-Jacobi equations[END_REF][START_REF] Siconolfi | Hamilton-Jacobi equations and Weak KAM Theory[END_REF].

Given a regular connected open domain Ω ⊂ R N , and a continuous Hamiltonian H : Ω×R N → R satisfying, for any x ∈ Ω, Z(x) := {p ∈ R N , H(x, p) ≤ 0}:

(H1) coercivity: Z(x) is compact; (H2) convexity: Z(x) is convex; (H3) H(x, 0) ≤ 0, i.e. 0 ∈ Z(x), we consider the following HJ equation

H(x, ∇u) = 0, x ∈ Ω.
(2.6)

We recall briefly some definitions related to viscosity theory. A continuous function u : Ω → R is said to be a viscosity subsolution (respectively supersolution) of (2.6) if H(x, ∇φ(x)) ≤ 0 (respectively H(x, ∇φ(x)) ≥ 0) for any x ∈ Ω and any C 1 function φ such that u -φ has a strict local maximizer (respectively minimizer) at x. Finally, u is a viscosity solution of (2.6) if it is both a subsolution and a supersolution. We denote by S - H (Ω) (respectively S + H (Ω)) the family of viscosity subsolutions (respectively supersolutions) of (2.6).

For x ∈ Ω, we define the support function of the 0-sublevel set Z(x) by

σ(x, q) := sup q • Z(x) = sup{q • p | p ∈ Z(x)} for q ∈ R N .
The assumption (H1)-(H3) ensures that σ is a possibly degenerate Finsler metric, i.e., σ is a continuous nonnegative function in Ω × R N , convex and positively homogeneous with the second variable q. Due to the assumption (H3), σ(x, q) is possible to equal to 0 for q = 0, which leads to the degeneracy and its dual σ * , as defined below, may take the value +∞. Here, the dual σ * (also called polar) is defined by

σ * (x, p) := sup q {p • q | σ(x, q) ≤ 1}.
Under the assumption (H1)-(H3), it is not difficult to see that there exists

K > 0 0 ≤ σ(x, q) ≤ K|q| for x ∈ Ω, q ∈ R N .
We denote by Γ(x, y) the set of Lipschitz curves defined on [0, 1] joining x, y in Ω. We then define the intrinsic distance by

d σ (x, y) := inf ζ∈Γ(x,y) 1 0 σ(ζ(t), ζ(t))dt, (2.7)
which is a quasi-distance, i.e. satisfying d σ (x, x) = 0 and the triangular inequality, but not necessarily symmetric. We summarize some basic characterizations of subsolution in terms of intrinsic distance d σ .

Proposition 2.2. ([12])

1) d σ is a quasi-distance in the sense that for any x, y ∈ Ω d σ (x, y) ≥ 0 and d σ (x, x) = 0. Moreover, for all x, y, z ∈ Ω one has

d σ (x, y) ≤ d σ (x, z) + d σ (z, y). 2) For any x ∈ Ω one has d σ (x, .) ∈ S - H (Ω) ∪ S + H (Ω \ {x}). 3) Compatibility condition: v ∈ S - H (Ω) if and only if v(x)-v(y) ≤ d σ (y, x) for any x, y ∈ Ω.
Given a closed subset D ⊂ Ω (typically D = ∂Ω or D = ω for some ω ⊂⊂ Ω), we consider the following HJ equation

H(x, ∇u) = 0 in Ω u = g on D (2.8)
where g : D → R is a continuous function satisfying the compatibility condition

g(x) -g(y) ≤ d σ (y, x) for any x, y ∈ D.
Thanks to Proposition 2.2, the unique maximal viscosity subsolution of the equation (2.8) can be recovered via the following maximization problem

max Ω u dx : u(x) -u(y) ≤ d σ (y, x), ∀x, y ∈ Ω and u = g on D , (2.9)
where d σ (., .) is the intrinsic distance associated to the Hamiltonian defined by (2.7).

Note that in the case of non-degenerate Finsler metric σ (i.e., H(x, 0) < 0), the 1-Lipschitz condition u(x) -u(y) ≤ d σ (y, x) for all x, y ∈ Ω can be characterized by u ∈ W 1,∞ (Ω), σ * (x, ∇u(x)) ≤ 1 a.e. x ∈ Ω. This property still holds true for degenerate Hamiltionians H (i.e., H(x, 0) ≤ 0) (cf. [START_REF] Hamza Ennaji | Augmented Lagrangian method for Hamilton-Jacobi equtions[END_REF]).

Main results.

In order to prove the connection between (2.9) and a Beckmann-type problem, we will consider a slightly more general variant of (2.9) by considering for ρ ∈ M b (Ω), the following maximization problem

(M D ) : max Ω u dρ : u(y) -u(x) ≤ d σ (x, y), ∀x, y ∈ Ω and u = g on D .
Then, it is not difficult to see that the solution of (2.9) can be recovered by taking ρ ≡ 1.

The connection with the Beckmann problem is given in the following theorem. min

φ∈M b (Ω) N ,ν∈M b (D) Ω σ(x, φ |φ| (x))d|φ| + D gdν : -div(φ) = ρ -ν in D (R N ) .
Moreover, u and (φ, ν) are optimal solutions to (M D ) and (BK), respectively, if and only if

               -div(φ) = ρ -ν in D (R N ) φ(x) • ∇ |φ| u(x) = σ x, φ |φ| (x) |φ| -a.e. x u = g on D,
where ∇ |φ| u denotes the tangential gradient with respect to |φ|, the total variation of φ (cf. [START_REF] Bouchitté | Characterization of optimal shapes and masses through Monge-Kantorovich equation[END_REF]).

In particular, we have the following.

Corollary 2.4. Let u be the maximal viscosity subsolution to (2.8) and (φ, ν) an optimal solution to (BK) with ρ = 1, then

               -div(φ) = 1 -ν in D (R N ) φ(x) • ∇ |φ| u(x) = σ x, φ |φ| (x) |φ| -a.e. x u = g on D.
For the case of non-degenerate Finsler metric σ (i.e., H(x, 0) < 0) we know that, by (1.5), the minimal value of (BK) is the same with

(M N ) : min ν∈M b (D) max u Ω ud(ρ -ν) + D gdν : u is 1 -Lipschitz w.r.t d σ ,
as well as Monge-Kantorovich problem

(MK) : min γ∈M + (Ω×Ω),ν∈M b (D) Ω×Ω d σ (x, y)dγ(x, y)+ D gdν : π 1 γ = ρ -+ν + , π 2 γ = ρ + +ν -.
The following theorem ensures that the above-mentioned relations still hold true for the case of degenerate HJ equation, i.e., H(x, 0) ≤ 0.

Theorem 2.5. Under the assumptions (H1-H3), we have

max(M D ) = min(BK) = min(M N ) = min(MK).
As a typical example we will consider an HJ equation of Eikonal type, coupled with a zero Dirichlet boundary condtion,

|∇u(x)| = k(x) in Ω u = 0 on ∂Ω.
In other words,

H(x, p) = |p| -k(x)
where k is a continuous, nonnegative function on Ω. In this case, the problem (M D ) can be rewritten as

sup Ω udx : |∇u| ≤ k and u = 0 ∂Ω , or (P) : inf u∈V {F(u) + G(Λu)} where V = C 1 (Ω) ∩ H 1 0 (Ω), F(u) = - Ω udx, Λu = ∇u and G(q) = 0 if |q| ≤ k +∞ otherwise.
For the case of non-degeneracy, i.e., k(x) > 0 on Ω, the Fenchel-Rockafellar duality, since the qualification conditions are satisfied (see e.g. [START_REF] Ekeland | Convex Analysis and Variational Problems[END_REF]Theorem III 4.1]), gives

sup Ω udx : |∇u| ≤ k and u = 0 on ∂Ω = min φ∈M b (Ω) N Ω k d|φ| : -div(φ) = ρ in D (Ω) .
However, dealing with general degenerate Hamiltionians at least two difficulties arise. Firstly, the qualification conditions are not satisfied to apply directly the Fenchel-Rockafellar duality. Secondly, in this setting the problem (BK) is not coercive, it follows that the existence of an optimal solution to (BK) is not trivial. These two issues will be addressed in this paper via approximation and optimal transport techniques. Before ending this section, let us mention that thanks to the duality result in [START_REF] Hamza Ennaji | Augmented Lagrangian method for Hamilton-Jacobi equtions[END_REF], we have the following Corollary 2.6. The extremal values (BK) and

(OF D ) : inf φ∈DM p (Ω) Ω σ(x, φ(x))dx -φ • n, g : -div(φ) = ρ in D (Ω) , coincide.
The formulation of the problem (OF D ) as well as the definition of DM p (Ω) and further comments are recalled in Remark 3.10.

Proofs

T : L 1 (Ω) N → R ∪ {+∞}, φ → T (φ) =      Ω σ(x, φ(x))dx if -div(φ) = ρ -ν in D (R N ) +∞ otherwise, E : Lip(Ω) → R ∪ {-∞}, u → E(u) =      Ω ud(ρ -ν) if σ * (x, ∇u(x)) ≤ 1 a.e. x ∈ Ω -∞ otherwise.
Lemma 3.7. Assume that σ is a degenerate Finsler metric. Let T , E be defined as above and

ρ(Ω) = ν(D). Then inf φ∈L 1 (Ω) N T (φ) = sup u∈Lip(Ω) E(u).
Proof. The proof will be divided into two steps. We first prove for the case of non-degenerate Finsler metric σ, i.e., there exist two positive constants K 1 , K 2 such that

K 1 |p| ≤ σ(x, p) ≤ K 2 |p| for any x ∈ Ω, p ∈ R N .
In this setting, due to the non-degeneracy of σ, the qualification conditions are satisfied and the result follows directly from the Fenchel-Rockafellar duality (see e.g. [START_REF] Ekeland | Convex Analysis and Variational Problems[END_REF]Theorem III 4.1]).

For the general case, we check at once that sup E ≤ inf T by taking u as a test function in the divergence constraint -div(φ) = ρ -ν in D (R N ). Therefore, it remains to prove that

inf T ≤ sup E. (3.10)
We now proceed by an approximation via non-degenerate Finsler metrics. For n ∈ N * and x ∈ Ω, define

σ n (x, p) := max(σ(x, p), |p| n ) for every p ∈ R N ,
which establishes a sequence of non-degenerate Finsler metrics σ n on Ω satisfying

|p| n ≤ σ n (x, p) ≤ K|p| for some constant K > 0.
Thanks to [7, Thereom 5.1], we have that d σn → d σ in the space of Finsler distances endowed with the topology induced by uniform convergence on compact subsets of Ω × Ω. To prove the inverse inequality (3.10), let us introduce for n ∈ N * the functionals

T n : L 1 (Ω) N → R ∪ {+∞}, φ →      Ω σ n (x, φ(x))dx if -div(φ) = ρ -ν in D (R N ) +∞ otherwise, E n : Lip(Ω) → R ∪ {-∞}, u →      Ω ud(ρ -ν) if σ * n (x, ∇u(x)) ≤ 1 a.e. x ∈ Ω -∞ otherwise.
It follows from the non-degeneracy of σ n and the first step of the proof that inf T n = sup E n . We are now in a position to show that sup E n → sup E as n → ∞. Let u n be a maximizer for E n , i.e., sup

E n = Ω u n d(ρ -ν) and σ * n (x, ∇u n (x)) ≤ 1 a.e. in Ω. Fix x 0 ∈ Ω. Since ρ(Ω) = ν(D), we see that u n := u n -u n (x 0
) is still a maximizer. Thus we can assume that u n (x 0 ) = 0 for any n. Since

u n (y) -u n (x) ≤ d σn (x, y) ≤ K|x -y| in Ω × Ω,
{u n } is equi-Lipschitz continuous. By Ascoli-Arzelà's theorem, there exists a Lipschitz function u such that, up to a subsequence,

u n ⇒ u uniformly in Ω. Since d σn → d σ as n → ∞ we deduce that u is admissible for E, i.e. u(y) -u(x) ≤ d σ (x, y) in Ω × Ω. Consequently lim n (sup E n ) = lim n Ω u n d(ρ -ν) = Ω ud(ρ -ν) ≤ sup E,
and inf T ≤ lim n (inf T n ) ≤ sup E
as claimed in (3.10).

Before ending up this subsection, we recall the notion of disintegration of measures which will be useful in the proof of existence of optimal solution to the Beckmann-type problem.

Theorem 3.8. (cf. [START_REF] Santambrogio | Optimal transport for applied mathematicians[END_REF]) Let X, Y be locally compact metric spaces and π : X → Y a Borel map. For any η ∈ M + b (X) there exist a family of probability measures (η y ) y∈Y on X concentrated on π -1 ({y}) such that for any test function u ∈ C(X), the mapping y → X udη y is Borel measurable and inf

X u(x)dη(x) = Y X u(x)dη y (x)dπ η(y).
φ∈L 1 (Ω) N ,ν∈M b (D) Ω σ(x, φ(x))dx + D gdν : -div(φ) = ρ -ν in D (R N ) .
Proof. First observe that

max(M D ) ≤ inf(BK) ≤ inf (BK). (3.11)
Indeed, take u satisfying σ * (x, ∇u(x)) ≤ 1 a.e. x in Ω, u = g on D as a test function in the divergence constraint -div(φ) = ρ -ν in D (R N ), we get

Ω udρ = Ω φ |φ| • ∇ud|φ| + D gdν ≤ Ω σ x, φ |φ| (x) d|φ| + D gdν.
This implies (3.11). So, it is sufficient to show the duality between (M D ) and (BK), i.e. max(M D ) = inf (BK). We use a perturbation technique as in [START_REF] De | Duality theory and optimal transport for sand piles growing in a silos[END_REF]. Define on C(D) the following functional

F : v ∈ C(D) → -max u Ω udρ : u ∈ Lip(Ω), σ * (x, ∇u(x)) ≤ 1, u + v = g on D ,
which is well-defined. Let us show that F is convex and l.s.c.. Consider

v 1 , v 2 ∈ C(D) and set v = tv 1 + (1 -t)v 2 for t ∈ [0, 1]. Let u 1 , u 2 ∈ Lip(Ω) be two maximizers corresponding to v 1 and v 2 respectively, i.e. σ * (x, ∇u i (x)) ≤ 1, u i + v i = g on D and F (v i ) = - Ω u i dρ for i = 1, 2. Define u = tu 1 + (1 -t)u 2 . It is evident that u + v = g on D.
And using the homogeneity of σ * , we obtain σ * (x, ∇u) ≤ tσ * (x, ∇u 1 ) + (1 -t)σ * (x, ∇u 2 ) ≤ 1 so that u is admissible for v. Finally, we get

F (v) ≤ - Ω udρ = t - Ω u 1 dρ + (1 -t) - Ω u 2 dρ = tF (v 1 ) + (1 -t)F (v 2 ),
which proves the convexity. For the lower semicontinuity, take a sequence v n ⇒ v uniformly on D. For every n ∈ N, consider a maximizer u n corresponding to v n such that u n + v n = g on D, u n are 1-Lipschitz w.r.t. d σ (i.e., u n (y) -u n (x) ≤ d σ (x, y) or equivalently, σ * (x, ∇u n (x)) ≤ 1 a.e. x ∈ Ω), and

F (v n ) = - Ω u n dρ.
Since u n (y) -u n (x) ≤ d σ (x, y) ≤ K|x -y|, the functions u n are equi-Lipschitz in the Euclidean distance. Since the sequence {v n } is convergent, it is bounded on D, and so is the sequence {u n }. By Ascoli-Arzelà's theorem, there exists a Lipschitz function u such that u n ⇒ u uniformly in Ω as n → ∞. It is clear that u + v = g on D and u(y) -u(x) ≤ d σ (x, y), i.e. u is admissible for v. The lower semicontinuity is completed by

F (v) ≤ - Ω udρ = lim n→∞ - Ω u n dρ = lim inf n→∞ F (v n ).
Since F is convex and l.s.c., we have F = F * * , in particular F (0) = F * * (0). Let us finish the proof by computing F * and F * * . For any ν ∈ M b (D), we see that

F * (ν) = sup v∈C(D) D vdν -F (v) = sup v∈C(D), u∈Lip(Ω) D vdν + Ω udρ : σ * (x, ∇u(x))) ≤ 1, u + v = g on D = sup u∈Lip(Ω) Ω udρ + D (g -u)dν : σ * (x, ∇u(x))) ≤ 1, g -u ∈ C(D) .
For any constant c ∈ R and a Lipschitz extension g of g, one can see that u := g + c is an admissible test function in the definition of F * (ν) and

Ω udρ + D (g -u)dν = c(ρ(Ω) -ν(D)) + Ω g dρ. So, if ρ(Ω) = ν(D), then F * (ν) = +∞. This implies that F * (ν) =        sup u∈Lip(Ω) Ω udρ + D (g -u)dν : σ * (x, ∇u(x))) ≤ 1 if ν(D) = ρ(Ω) +∞ otherwise.
Following Lemma 3.7, for any ν ∈ M b (D) such that ν(D) = ρ(Ω), we have

F * (ν) = sup u∈Lip(Ω) Ω udρ + D (g -u)dν : σ * (x, ∇u(x))) ≤ 1 = D gdν + sup u∈Lip(Ω) Ω ud(ρ -ν) : σ * (x, ∇u(x))) ≤ 1 = D gdν + inf φ∈L 1 (Ω) N Ω σ(x, φ(x))dx : -div(φ) = ρ -ν in D (R N ) = inf φ∈L 1 (Ω) N Ω σ(x, φ(x)) dx + D gdν : -div(φ) = ρ -ν in D (R N ) .
Consequently,

max(M D ) = -F (0) = -F * * (0) = -sup ν∈M b (D) -F * (ν) = inf (BK).
Remark 3.10. Following our approach in [START_REF] Hamza Ennaji | Augmented Lagrangian method for Hamilton-Jacobi equtions[END_REF], it is possible to show that the optimal ν in (BK) is somehow related to the trace of the optimal flow φ. For completeness let us recall briefly some of the main ingredients we used, and for a simpler presentation, consider the case where D = ∂Ω. Then we prove in [START_REF] Hamza Ennaji | Augmented Lagrangian method for Hamilton-Jacobi equtions[END_REF] the duality between (M D ) and

(OF D ) : inf φ∈DM p (Ω) Ω σ(x, φ(x))dx -φ • n, g : -div (φ) = ρ in D (Ω) ,
where we define for 1 ≤ p ≤ ∞,

DM p (Ω) := φ ∈ L p (Ω) N : div φ =: µ ∈ M b (Ω) ,
endowed with the graph norm

φ DM p (Ω) := φ L p (Ω) + |div φ|(Ω).
The main interest of introducing such a space is to give a sense to the trace term φ•n which is not always defined for a general measure field φ. In particular, for any measure field φ ∈ DM p (Ω), one can define a trace φ • n on ∂Ω as a linear form on Lip(∂Ω) such that

φ • n, ξ /∂Ω = Ω ξ divφ + Ω ∇ξ • φ, for any ξ ∈ C 1 (Ω).
One can see that, at least formally, ν plays the role of -φ • n in (OF D ).

Our aim now is to use the optimal mass transportation interpretation to prove that the inf in (BK) is actually a min, i.e. the existence of optimal solution to the Beckmann-type problem (BK). In particular we see that min(BK) = max(M D ).

Proof. Take (γ, ν) be a solution of (MK) and define a vector measure Φ γ through

< Φ γ , v >= Ω×Ω 1 0 v(ξ(t)) ξ(t)dtdγ(x, y) ∀v ∈ C(Ω) N ,
with ξ being a geodesic joining x and y with respect to d σ . Let us check the feasibility of Φ γ for (BK), i.e.

-div

(Φ γ ) = ρ -ν in D (R N ). (3.12)
For any w ∈ C 1 (Ω), by definition, we have

< Φ γ , ∇w > = Ω×Ω 1 0 dw(ξ(t)) dt dtdγ(x, y) = Ω×Ω (w(y) -w(x))dγ(x, y) = Ω wd(ρ -ν),
which gives (3.12). The next task is to show that

Ω σ(x, Φ γ |Φ γ | (x))d|Φ γ | + D gdν ≤ Ω×Ω d σ (x, y)dγ(x, y) + D gdν.
To do so, for each t ∈ [0, 1], define vector measure E t by setting E t (v) := Ω×Ω v(ξ(t)) ξ(t)dγ(x, y)

for v ∈ C(Ω) N . We get Φ γ = 1 0
E t dt and, by Jensen's inequality,

Ω σ(x, Φ γ |Φ γ | (x))d|Φ γ | ≤ 1 0 Ω σ(x, E t |E t | (x))d|E t |dt. (3.13) Now define π t : (x, y) ∈ Ω × Ω → ξ(t) for t ∈ [0, 1],
where ξ is as before, a geodesic joining x and y. Consider η t = (π t ) γ. Using disintegration theorem (see Theorem 3.8) for γ with respect to π t , we can find probability measures γ z supported on π -1 t ({z}) such that

Ω×Ω u(x, y)dγ(x, y) = Ω×Ω×Ω u(x, y)dγ z (x, y)dη t (z). (3.14)
We check at once that E t << η t with the density dE t dη t (z) = Ω×Ω ξ(t)dγ z (x, y), which follows from the fact that, for any test function v,

Ω v(z)dE t = Ω v(ξ(t)) ξ(t)dγ(x, y) = Ω×Ω×Ω v(ξ(t)) ξ(t)dγ z (x, y)dη t (z) = Ω v(z) Ω×Ω ξ(t)dγ z (x, y)dη t (z).
On the other hand, by (3.14) and Jensen's inequality,

Ω×Ω σ(ξ(t), ξ(t))dγ(x, y) = Ω×Ω×Ω σ(ξ(t), ξ(t))dγ z (x, y)dη t (z) = Ω×Ω×Ω σ(z, ξ(t))dγ z (x, y)dη t (z) ≥ Ω σ z, Ω×Ω ξ(t)dγ z (x, y) dη t (z) = Ω σ(z, E t |E t | (z))d|E t |(z). (3.15)
Finally, we observe that since ξ is a geodesic

Ω×Ω d σ (x, y)dγ(x, y) + D gdν = Ω×Ω 1 0 σ(ξ(t), ξ(t))dγ(x, y)dt + D gdν ≥ 1 0 Ω σ(z, E t |E t | (z))d|E t |(z) + D gdν (by (3.15)) ≥ Ω σ(x, Φ γ |Φ γ | (x))d|Φ γ | + D gdν (by (3.13)). Consequently inf(BK) ≤ Ω σ(x, Φ γ |Φ γ | (x))d|Φ γ | + D gdν ≤ Ω×Ω d σ (x, y)dγ(x, y) + D gdν = min(MK).
Proof of Theorem 2.5. Following Theorem 1.1 on Kantorovich duality, one has max(M N ) = min(MK). Combining this with Propositions 3.9 and 3.11, we conclude that the Beckmann-type problem (BK) admits an optimal solution and max(M D ) = min(BK) = max(M N ) = min(MK).

Proof of Theorem 2.3. The proof of the duality between max(M D ) and inf(BK) is followed from Proposition 3.9 while the existence of optimal solution to the Beckmann-type problem (BK) is a consequence of Proposition 3.11 (see also Theorem 2.5). Let us now show the optimality conditions. Indeed, u and (φ, ν) are optimal solutions for (M D ) and (BK), respectively, if and only if

Ω σ(x, φ |φ| (x))d|φ| + D gdν = Ω udρ, or Ω σ(x, φ |φ| (x))d|φ| = Ω ud(ρ -ν) = Ω φ(x) • ∇ |φ| u(x)d|φ|.
This is equivalent to

σ(x, φ |φ| (x)) = φ(x) • ∇ |φ| u(x) for |φ|-a.e. x,
as desired.

Numerical results

As we pointed out in the previous sections, the maximization problem (M D ) is linked to the Monge-Kantorovich type problem (MK). The measure ρ needs not to satisfy the standard mass balance condition. However, transporting a part of the mass from/to the Dirichlet region D is allowed. In addition, taking ρ = L N |Ω the restriction of Lebesgue measure on Ω, (M D ) allows recovering the solution of HJ equation. This was essentially addressed in [START_REF] Hamza Ennaji | Augmented Lagrangian method for Hamilton-Jacobi equtions[END_REF] for several Hamiltonians. Here, we will focus essentially on the solution of (BK).

Formulation of the problem. We set

X = W 1,∞ (Ω) and Y = L ∞ (Ω) N × C(D), with F(u) = - Ω u dρ, Λu = (∇u, u |D ), for any u ∈ X and G(η, h) =    0 if σ * (x, η) ≤ 1 and h = g on D +∞ otherwise, for any (η, h) ∈ Y .
Thus, we can rewrite the problem (M D ) in the form

-inf u∈X F (u) + G(Λu).
Thanks to Theorem (2.5) and Proposition (3.9), we have

-min F(u)+G(Λu) : u ∈ X = -sup -F * (-Λ * (φ, ν))-G * (φ, ν) : φ ∈ L 1 (Ω) N , ν ∈ M b (D) = min φ∈L 1 (Ω) N ,ν∈M b (D) Ω σ(x, φ(x))dx + D gdν : -div (φ) = ρ -ν in D (R N ) .
Introducing a new primal variable (p, q) ∈ Y , (M D ) reads L(u, (p, q); φ, ν)

- inf (u,(p,q))∈X ×Y Λu=(p,q) F(u) + G(p, q).
where

L(u, (p, q); φ, ν) = F(u) + G(p, q) + Ω φ • (∇u -p) dx + D ν • (u -q)dx.
As usual, it is convenient to consider the following augmented Lagrangian L r (u, (p, q); φ, ν) = L(u, (p, q); φ, ν)

+ r 2 |∇u -p| 2 + r 2 |u |D -q| 2 , r > 0,
which has the same saddle points as L. Thus, the problem we will focus on is

(S r ) : min (u,q)∈X ×Y sup φ∈L 1 (Ω) N ,ν∈M b (D)
L r (u, (p, q); φ, ν).

The numerical approximation of (M D ) was addressed in [START_REF] Hamza Ennaji | Augmented Lagrangian method for Hamilton-Jacobi equtions[END_REF]. The main difference with the above formulation and the one in [START_REF] Hamza Ennaji | Augmented Lagrangian method for Hamilton-Jacobi equtions[END_REF] is the presence of an additional Lagrange multiplier ν corresponding to the boundary condition u = g on D. The resolution of (M D ) in [START_REF] Hamza Ennaji | Augmented Lagrangian method for Hamilton-Jacobi equtions[END_REF] was used via ALG2 algorithm [START_REF] Ekeland | Convex Analysis and Variational Problems[END_REF][START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF] and was implemented using finite element method. Up to our knowledge, it is not straightforward to define finite element functions ν on a closed subset D (typically D = ∂Ω). In the next section we explain how to tackle this difficulty. 4.2. Practical implementation. As we pointed out in the previous subsection, it is not clear how to solve the saddle point problem (S r ) via ALG2 algorithm due to the presence of the measure ν on the Dirichlet region. However, thanks to Corollary 2.6 and Remark 3.10 which shows that ν is somehow linked to the trace of the optimal flow φ on D, we can restrict ourselves to a formulation involving only the potential u and the flow φ. More precisely, for any u ∈ X and η ∈ Z = L ∞ (Ω), we take

F(u) =      - Ω u dρ if u = g on D +∞ otherwise , G(η) =    0 if σ * (x, η) ≤ 1 +∞ otherwise
, and Λu = ∇u.

Thus, following our approach in [START_REF] Hamza Ennaji | Augmented Lagrangian method for Hamilton-Jacobi equtions[END_REF], we can focus on the following saddle point problem inf (u,q)∈X ×Z sup φ∈DM p (Ω\D)

L(u, q; φ)

where

L(u, q; φ) = F(u) + G(q) + Ω φ • (Λ(u) -q) dx, for any (u, q, φ) ∈ X × Z × DM p (Ω \ D).
Hence, the augmented Lagrangian reads L r (u, q; φ) = F(u) + G(q) + φ, Λu -q + r 2 |Λu -q| 2 , for any (u, q, φ) ∈ X × Z × DM p (Ω \ D).

The definition of the space DM p (Ω \ D) is recalled in Remark 3.10 and the application of ALG2 algorithm for the augmented Lagrangian L r can be found in [START_REF] Hamza Ennaji | Augmented Lagrangian method for Hamilton-Jacobi equtions[END_REF]. The interested reader can check [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF][START_REF] Ekeland | Convex Analysis and Variational Problems[END_REF][START_REF] Glowinski | Numerical methods for nonlinear variational problems[END_REF][START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF][START_REF] Fortin | Augmented Lagrangian methods: applications to the numerical solution of boundary-value problems[END_REF][START_REF] Igbida | Augmented Lagrangian method for optimal partial transportation[END_REF][START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF][START_REF] Benamou | Augmented Lagrangian methods for transport optimization, mean field games and degenerate elliptic equations[END_REF] for more details on saddle-point problems and ALG2 algorithm. For the first test we take k(x, y) = 1 and ρ + = 2, ρ -= δ (0.5,0.5) . In the second test we take ρ + = 4χ and k(x, y) = 5 -3e -2 * ((x-0.5) 2 +(y-0. In the third test we take two Gaussian densities ρ + = e -40 * ((x-0.75) 2 +(y-0.3) 2 ) and ρ -= e -40 * ((x-0.3) 2 +(y-0.65) 2 ) . We change the distance d σ by taking a degenerate k. In particular, we choose

IsoValue

k(x, y) = (1 -2x) 2 (y -y 2 ) 2 + (1 -2y) 2 (x -x 2 ) 2 χ B where B = {(x, y) ∈ [0, 1] 2 :
(x -0.5) 2 + (y -0.5) 2 > 0.25}. We clearly see that the flux is concentrated essentially on the region where k vanishes, i.e, on B c , which in terms of optimal transport with respect to d σ represents a free transport region. In the last test we consider a Finsler metric of cystalline type, namely

σ(v) = max i=1,••• ,5 v • d i with d 1 = (1, -1), d 2 = (1, -0.8), d 3 = (-0.8, 1), d 4 = (-1, 1), d 5 = (-1, -1

Comments and extentions

A natural extension one can think of is the HJ equation with double obstacles on the boundary. More precisely one can consider the equation H(x, ∇u) = 0 on Ω, g 1 ≤ u ≤ g 2 on ∂Ω where g i : ∂Ω → R are continuous functions satisfying the compatibility condition g 1 (x)-g 2 (y) ≤ d σ (y, x) for any x, y ∈ ∂Ω.

In order to establish the link between this problem and a Bekmann-type problem, we consider as previously the following maximization problem (M D ) o : max Ω u dρ : u ∈ W 1,∞ (Ω), σ * (x, ∇u) ≤ 1 and g 1 ≤ u ≤ g 2 on ∂Ω .

Similarly, we can state the following result. i.e. the equality holds in (5.16), which is equivalent to the system of optimality conditions as desired.

Let us mention that problems of the form (M D ) o and (BK) o can arise when studying optimal transport problems with some import/export costs on the boundary. In the case where the transport cost is given by the Euclidean distance c(x, y) = |x -y| (in our case when considering the Eikonal equation |∇u| = 1) we refer the reader to [START_REF] José | An optimal transportation problem with a cost given by the Euclidean distance plus import/export taxes on the boundary[END_REF].

(1. 5 ) 1 . 2 .

 512 Hamilton-Jacobi equation and Beckmann's problem. The connection between HJ equation coupled with a Dirichlet condition (1.1) and Beckmann's problem is not straightforward. In the case where H(x, p) = |p| -k(x) (i.e., the Eikonal equation |∇u(x)| = k(x)) complemented with Dirichlet boundary condition u = g on ∂Ω, the viscosity solution of (1.1) can be characterized through the following optimization problem max Ω u dx : u(x) -u(y) ≤ d k (y, x) and u = g on ∂Ω ,

Theorem 2 . 3 .

 23 The optimization problem (M D ) coincides with the following Beckmann-type problem (BK) :

3. 1 .

 1 Preparatory results. Let ν ∈ M b (D) satisfy ρ(Ω) = ν(D). Define two functionals

3. 2 .

 2 Proofs of the main results. We get started with the proof of Theorem 2.3 by the following result. Proposition 3.9. We have max(M D ) = inf(BK) = inf (BK), where (BK) :

Proposition 3 . 11 .

 311 There exist ν ∈ M b (D) and a vector measure Φ such that -div(Φ) = ρ -ν in D (R N ) and inf(BK) ≤ Ω

  This allows us to rewrite (M D ) and (BK) in a saddle point form (S ) : min (u,(p,q))∈X ×Y sup φ∈L 1 (Ω) N ,ν∈M b (D)

4. 3 .

 3 Some examples. We take Ω = [0, 1] × [0, 1] and g = 0 on ∂Ω. The first three examples are performed with a Finsler distance d σ of Riemannian type d σ (x, y) t))| ζ(t)|dt.

Figure 1 .

 1 Figure 1. (A): the potential u, (B): the flow φ.

Figure 2 .

 2 Figure 2. (A): the potential u, (B): the flow φ.

Figure 3 .

 3 Figure 3. (A): the potential u, (B): the flow φ privileging the zero set of k .

Figure 4 .

 4 Figure 4. (A): the potential u, (B): the flow φ .

Theorem 5 . 12 .-

 512 The optimization problem (M D ) o coincides with the following Beckmann-type problem, denoted by (BK) o ,min φ∈M b (Ω) N ,ν∈M b (∂Ω) Ω σ(x, φ |φ| (x))d|φ| + ∂Ω g 2 dν + -∂Ω g 1 dν -: -div(φ) = ρ -ν in D (R N ) .Moreover, u and (φ, ν) are optimal solutions to (M D ) o and (BK) o , respectively if and only if div(φ) = ρ -ν in D (R N ) φ(x) • ∇ |φ| u(x) = σ x, φ |φ| (x)for |φ| -a.e. x g 1 ≤ u ≤ g 2 on ∂Ω and u = g 1 for ν --a.e. x and u = g 2 for ν + -a.e. x.Sketch of proof. We define on C(∂Ω) × C(∂Ω) the following functional on C(∂Ω) × C(∂Ω) by, for(v, w) ∈ C(∂Ω) × C(∂Ω), E(v, w) = -sup Ω udρ : u ∈ Lip(Ω), σ * (x, ∇u(x)) ≤ 1, g 1 ≤ u + v, u + w ≤ g 2 on ∂Ω .Most of the arguments of Section 3 can be reproduced to show that E is convex and lower semicontinuous, which gives E(0, 0) = E * * (0, 0). It follows that max(M D ) o = min(BK) o . We now turn to the optimality conditions. First observe that, for any feasible u and (φ, ν), we get used∇ |φ| u(x)φ(x) ≤ σ(x, φ |φ|(x)) for |φ| -a.e. x by the fact that u is 1-Lipschitz w.r.t. d σ . On the other hand, u and (φ, ν) are optimal for (M D ) o and (BK) o respectively, if and only if

  [(x-0.3) 2 +(y-0.2) 2 <0.03] and ρ -= 4χ [(x-0.7) 2 +(y-0.8) 2 <0.03]

  5) 2 ) .

	IsoValue	IsoValue
	-0.799859	-4.22222
	-0.684195	-3.77778
	-0.607086	-3.33333
	-0.529977	-2.88889
	-0.452868	-2.44444
	-0.375759	-2
	-0.29865	-1.55556
	-0.221541	-1.11111
	-0.144432	-0.666667
	-0.0673228	-0.222222
	0.00978628	0.222222
	0.0868954	0.666667
	0.164004	1.11111
	0.241113	1.55556
	0.318223	2
	0.395332	2.44444
	0.472441	2.88889
	0.54955	3.33333
	0.626659	3.77778
	0.819432	4.22222
		Vec Value
		0
		0.0578684
		0.115737
		0.173605
		0.231474
		0.289342
		0.34721
		0.405079
		0.462947
		0.520816
		0.578684
		0.636552
		0.694421
		0.752289
		0.810157
		0.868026
		0.925894
		0.983763
		1.04163
		1.0995
	(a)	(b)

  ).

	IsoValue	IsoValue
	-0.162111	-4.22222
	-0.138935	-3.77778
	-0.123484	-3.33333
	-0.108033	-2.88889
	-0.0925827	-2.44444
	-0.077132	-2
	-0.0616814	-1.55556
	-0.0462308	-1.11111
	-0.0307801	-0.666667
	-0.0153295	-0.222222
	0.000121148	0.222222
	0.0155718	0.666667
	0.0310224	1.11111
	0.0464731	1.55556
	0.0619237	2
	0.0773743	2.44444
	0.092825	2.88889
	0.108276	3.33333
	0.123726	3.77778
	0.162353	4.22222
		Vec Value
		0
		0.101396
		0.202792
		0.304187
		0.405583
		0.506979
		0.608375
		0.70977
		0.811166
		0.912562
		1.01396
		1.11535
		1.21675
		1.31814
		1.41954
		1.52094
		1.62233
		1.72373
		1.82512
		1.92652
	(a)	(b)