Semi-discrete optimal transport methods for the semi-geostrophic equations

We give a new and constructive proof of the existence of global-in-time weak solutions of the 3-dimensional incompressible semi-geostrophic equations (SG) in geostrophic coordinates, for arbitrary initial measures with compact support. This new proof, based on semi-discrete optimal transport techniques, works by characterising discrete solutions of SG in geostrophic coordinates in terms of trajectories satisfying an ordinary differential equation. It is advantageous in its simplicity and its explicit relation to Eulerian coordinates through the use of Laguerre tessellations. Using our method, we obtain improved time-regularity for a large class of discrete initial measures, and we compute explicitly two discrete solutions. The method naturally gives rise to an efficient numerical method, which we illustrate by presenting simulations of a 2-dimensional semi-geostrophic flow in geostrophic coordinates generated using a numerical solver for the semi-discrete optimal transport problem coupled with an ordinary differential equation solver.

In this paper we consider SG in geostrophic coordinates, associated to flows on an arbitrary convex bounded (physical) domain Ω ⊂ R 3 , which we interpret as the active transport equation for the time-dependent measure-valued map α, which is known as the potential vorticity. The connection between SG and optimal transport is contained in the nonlocal divergence-free velocity field W [α], which is often called the geostrophic velocity. Let B[α t ] denote the unique mean-zero convex function whose gradient is the optimal transport map between the Lebesgue measure on Ω and the Borel measure α t with respect to the quadratic cost, and let B[α t ] * denote its Legendre-Fenchel transform on Ω. At each time t, W [α] is given by where id R 3 denotes the identity on R 3 and J := Guided by the work of Cullen and Purser [14], this connection was first established rigorously by Benamou and Brenier. In [4], those authors proved the existence of global-in-time weak solutions of (1) for initial measures which are absolutely continuous with respect to the Lebesgue measure and have compactly-supported L p density for p > 3. This result was extended by Lopes Filho and Nussenzveig Lopes in [41] to the case where p ≥ 1, and by Loeper [40] and later Feldman and Tudorascu [23] to the case where the initial measure need only have compact support in R 3 . In [22,Proposition 4.14], Feldman and Tudorascu use Ambrosio and Gangbo's abstract techniques for Hamiltonian ODEs in Wasserstein space [3] to prove that when the initial measure is an arbitrary convex combination of Dirac masses there exists a global-in-time solution that maintains the discrete structure of the initial data. The two known results regarding the uniqueness of solutions of (1) are the local-in-time uniqueness of Hölder continuous periodic solutions proved in [40], and weak-strong uniqueness under uniform convexity proved in [24]. Otherwise, the problem of uniqueness of solutions remains open.
The main contribution of this paper is an alternative proof the existence of global-in-time weak solutions of (1) for arbitrary compactly-supported initial measures, which uses recently developed techniques from semi-discrete optimal transport to treat the case where the initial measure is discrete (see Sect. 3 and, in particular, Theorems 3.5 and 3.6). For a wide class of discrete initial measures our result recovers [22,Proposition 4.14] with improved time regularity (twice continuously differentiable rather than Lipschitz) and uniqueness. More significantly, our application of semi-discrete optimal transport to SG illuminates an explicit and intuitive connection between geostrophic coordinates and corresponding flows in the physical domain Ω. It also gives a constructive way of determining solutions explicitly, and it forms the basis of an effective numerical scheme, as we illustrate in Sect. 7.

SG in geostrophic coordinates and semi-discrete optimal transport
In this section, we describe our approach to studying (1) using semi-discrete optimal transport, which is the special case of optimal transport in which the source measure is absolutely continuous with respect to the Lebesgue measure and the target measure is discrete.
In recent years, semi-discrete optimal transport theory has seen significant expansion in its theoretical foundations (see [5,15,17,30,32,36,38,[42][43][44][45]). It has also been applied to many diverse problems in the sciences, both within fluid dynamics [28,37] and elsewhere such as materials science [6,7,33], economics [27,Chapter 5], crowd dynamics [34] and image interpolation [36]. Inspired by ideas in the original work of Cullen and Purser [14] on piecewise constant solutions of semi-geostrophic slice models, we use semi-discrete optimal transport to analyse (1) in the special case where the initial potential vorticity α 0 = α is a discrete measure, i.e., for some m i > 0 and z i ∈ R 3 . We show (Theorem 3.5) that for well-prepared discrete initial data (see Definition 3.4) there exists a corresponding discrete solution t → α t of (1) of the form where the trajectories z i are twice continuously differentiable. Denoting by L 3 the Lebesgue measure on R 3 and by L 3 ¬ Ω its restriction to Ω, the optimal transport map between L 3 ¬ Ω and α t given by (3) is a piecewise constant function is a tessellation of Ω by convex sets, known as the optimal Laguerre tesselation (see Definitions 2.4 and 2.5) generated by the seed vector z := (z 1 , . . . , z N ) subject to the mass constraint Letting x i (z) denote the centroid of the Laguerre cell C i (z), one can show (see Lemma 4.2) that the time-dependent measure-valued map defined by (3) is a weak solution of (1) if and only if the trajectories z 1 , . . . , z N satisfy the ODE initial value problem (IVP) for i ∈ {1, . . . , N }. At each time t, the seeds z i (t) in geostrophic space generate a Laguerre tessellation of the physical domain Ω (see Fig. 1 for a 2D illustration).
We obtain a solution of (1) for an arbitrary compactly-supported initial measure α by generating a sequence (α N ) N ∈N of well-prepared discrete measures converging to α in the Wasserstein 2-distance, evolving each of these discrete measures according to the corresponding ODE-IVP (4), and using compactness in the space of continuous measure-valued maps to pass to the limit as N → ∞. As such, our construction method can be thought of as a Fig. 1 Typical snapshot at a time t of a 2D discrete solution α of (1) and the corresponding tessellation of the physical domain Ω, which is taken here to be a square subset of R 2 . The right-hand plot shows a configuration of N = 8 seeds z(t) = (z 1 (t), . . . , z 8 (t)) in 2D geostrophic space, on which the measure α t is supported. The left-hand plot shows the (approximate) Laguerre tessellation of Ω generated by z(t) subject to the constraint that all cells have the same mass. Black lines represent Laguerre cell boundaries and the circle within each cell C i represents its centroid x i (z(t)) meshless or particle method. By comparison with the proof given in [4], and later generalised in [40,41], the discretisation occurs in the spatial domain rather than in the time domain. Note that in [13] Cullen, Gangbo and Pisante analyse a variant of SG using a spatial discretisation different from the one considered in this paper. Analytically, the essential benefit of the discretisation of the initial measure α is that the study of the active transport equation (1), whose velocity field is only in general of class BV loc , is replaced by the study of the ODE-IVP (4), whose right hand side is continuously differentiable L 3N -almost everywhere. Mollifications of the vector-field and related quantities used in [4,40,41], as well as the abstract techniques for Hamiltonian ODEs in Wasserstein space used in [22], are therefore avoided, resulting in a more direct solution procedure.

Background on the semi-geostrophic equations
In this section, we briefly describe how Eq. (1) is derived. In their traditional Eulerian formulation in a fixed spatial domain Ω ⊂ R 3 , the (non-dimensionalised) semi-geostrophic equations are given by the coupled system ⎧ ⎨ where u g := (u g,1 , u g,2 , 0) T is the geostrophic velocity field, u is the Eulerian velocity field of the fluid, u a := u − u g is the ageostrophic velocity field, θ is the potential temperature, and the matrix J , defined by (2), encodes planetary rotation. Importantly, the hydrodynamic and thermodynamic fields u g and θ are linked through the fluid pressure p by the identity When posed on a suitably smooth bounded domain Ω ⊂ R 3 , the system (5) is typically supplemented with a no-slip boundary condition u · n = 0 on ∂Ω, where n is the outward unit normal field on the boundary, which is sufficient to ensure that fluid points remain in the domain Ω for all times. Define the geopotential P pointwise by for x ∈ Ω and t ≥ 0. The system (5) can then be expressed as where id Ω denotes the identity map on Ω, and U : ∇ P → u is the formal solution operator associated to the (time-independent) div-curl boundary-value problem By way of this simple change of dependent variable, SG can be viewed as an inhomogeneous active transport equation (7) whose unknown ∇ P is a time-dependent conservative vector field on Ω. The Eulerian velocity field is then formally defined through the action of the solution operator U . This change of dependent variable also highlights a substantial mathematical difficulty one faces when constructing solutions of (7): for the boundary-value problem (8) to be of elliptic type at each time t, P(·, t) must be strictly convex.
The state-of-the-art regarding the existence of solutions of SG in Eulerian coordinates is due to Ambrosio et al. [1]. Using the W 2,1 loc -regularity of Alexandrov solutions of a class of Dirichlet boundary-value problems for the Monge-Ampère equation established in [16], the authors proved the existence of global-in-time distributional solutions of SG in Eulerian coordinates posed on smooth convex domains Ω ⊂ R 3 for a class of initial geopotentials P 0 satisfying However, for such solutions, the support of the pushforward measure ∇ P(·, t) # L 3 ¬ Ω is the whole space R 3 at each time t. The relation (6) then implies that the temperature field θ satisfies θ(·, t) / ∈ L ∞ (Ω). Interpreted physically, this means that the atmospheric fluid is arbitrarily hot on sets of positive measure at all times. At the time of writing, the existence of either local-in-time or global-in-time distributional solutions of (7) for physical initial data P 0 satisfying ∇ P 0 (Ω) ⊂⊂ R 3 remains open.
Since the pioneering work of Hoskins [31], Cullen and Purser [14], and Benamou and Brenier [4] on the semi-geostrophic equations, it has become customary to regard ∇ P(·, t) formally as a diffeomorphism between Ω and its image ∇ P(Ω, t) for each time t. The system (5) is then transformed to the time-dependent coordinate system determined by ∇ P, known as geostrophic coordinates. It is a remarkable property of SG that, as shown in [4], this formal change of coordinates yields a closed equation which is free of the field u. Indeed, under the assumption that ∇ P is a smooth solution of (7) and P(·, t) is strictly convex at each time t, it can be shown that the time-dependent pushforward measure is a distributional solution of the active transport equation (1).

Outline of the paper
We begin in Sect. 2 with a brief introduction to semi-discrete optimal transport theory. Section 3 contains the statement of the following existence results whose novel proofs are the main contribution of this paper: 1. discrete geostrophic solutions with well-prepared discrete initial data exist, are unique, and are defined by trajectories that are twice continuously differentiable in time (see Definition 3.4 and Theorem 3.5); 2. Lipschitz-in-time solutions of SG in geostrophic coordinates with arbitrary compactlysupported initial measure can be constructed as the uniform limit of a sequence of discrete geostrophic solutions that are twice continuously differentiable in time (Theorem 3.6).
These results are proved in Sects. 4 and 5 respectively. Section 6 contains the explicit calculation of two exact solutions of SG in geostrophic coordinates, as well as a brief discussion on equilibrium solutions. Finally, in Sect. 7, we illustrate the theory developed in the paper by simulating a 2D semi-geostrophic flow in geostrophic coordinates, and we plot the corresponding Laguerre tessellations of the physical domain Ω.

Notation
Let d ∈ N. We denote by R d > the subset of R d consisting of all vectors whose components are positive. For i ∈ {1, . . . , d}, the i th canonical basis vector in R d is denoted by e i . Let A ⊆ R d be a Borel set. We denote the identity map on A by id A , and the characteristic function of A by 1 A . We denote the interior of A by Int(A) and the boundary of A by ∂ A.
Measures We denote by L d the Lebesgue measure on R d and by L d ¬ A its restriction to A. The set of Borel probability measures on A is denoted by P( A). Given a Borel map T : A → R d and a measure μ ∈ P( A), the pushforward of μ by T is denoted by T # μ and is defined by T # μ(B) = μ(T −1 (B)) for all Borel sets B ⊆ R d . The set of Borel probability measures on R d with compact support is denoted by P c (R d ). For any p ∈ [1, +∞), P p (R d ) denotes the set of all Borel probability measures on R d with finite moments of order p, equipped with the Wasserstein p-distance W p . This is defined for μ, ν ∈ P p (R d ) by where π x and π y denote the projections onto the first and second variables, respectively. Throughout this paper spaces of probability measures are understood to be equipped with the Wasserstein 2-distance unless otherwise stated.

Convex functions
Given a convex function f : A → R, the subdifferential of f is the set-valued function, mapping from A into the set of subsets of R d , defined by The Legendre-Fenchel transform of f is the function f * :

Test functions
We denote by D(R d ) the space of test functions C ∞ c (R d ) equipped with the standard semi-norm topology (see, for example, [26]).
Physical domain Throughout this paper (with the exception of Sect. 7) Ω is taken to be an arbitrary convex open bounded subset of R 3 , and, without loss of generality, we use the normalisation convention that L 3 (Ω) = 1 so that all measures under consideration are probability measures.

Semi-discrete optimal transport
In this section we review some basic aspects of semi-discrete optimal transport theory. For further information on semi-discrete optimal transport see [44,Section 4] and, for more on optimal transport theory in greater generality, see [48,51,52].
Given a target measure ν ∈ P 2 (R 3 ), a Borel map T : Ω → R 3 is said to be an optimal transport map between L 3 ¬ Ω and ν with respect to the quadratic cost c : subject to the constraint that The problem of finding an optimal transport map given source and target measures is known as the Monge problem. For any ν ∈ P 2 (R 3 ) such a map T exists, is unique, and can be expressed as the gradient of a convex function Φ belonging to the Sobolev space H 1 (Ω) (see for example [48,Theorem 1.22] or [51, Theorem 2.12]).

Definition 2.1
We define the operator B : P 2 (R 3 ) → H 1 (Ω) to be that which sends any given ν ∈ P 2 (R 3 ) to the unique mean-zero convex function in H 1 (Ω) whose gradient is the unique optimal transport map between L 3 ¬ Ω and ν with respect to the quadratic cost.
The stability of optimal transport, including the continuity of B, has been studied in, for example, [5,38,43] and [52,Theorem 5.20 and Corollary 5.23]. In [43] it is shown that for any bounded set A ⊂ R 3 , the restriction of B to P( A) is Hölder continuous.
Semi-discrete optimal transport (see, for example, [32], [44,Section 4], [48, Section 6.4.2]) refers to the special case where, for some N ∈ N, the target measure ν belongs to the class where

Definition 2.3 We call a vector
Hence the optimal transport problem between L 3 ¬ Ω and ν ∈ Q N (R 3 ) is reduced to an optimal partitioning problem. Moreover, the unique optimal partition, and corresponding transport map, can be characterised using the notion of Laguerre tessellations. Definition 2.4 (Laguerre tessellation) Given a seed vector z = (z 1 , . . . , z N ) ∈ R 3N and a weight vector w = (w 1 , . . . , w n ) ∈ R N , the Laguerre tessellation of Ω generated by the pair (z, w) is defined to be the family Note that any Laguerre cell is convex since it is the intersection of finitely many half-spaces with the convex set Ω. In particular, Laguerre cells that do not intersect ∂Ω are polyhedra. Moreover, for any seed vector z, weight vector w and indices i = j, the intersection Using the Kantorovich Duality Theorem (see, for example, [48, Section 1.2]), one can show that the optimal transport cost between L 3 ¬ Ω and ν is the supremum over all weight vectors w = (w 1 , . . . , w N ) ∈ R N of the Kantorovich functional g : R N → R defined by The Kantorovich functional g is concave, and maximisers w ∈ R N of g exist and satisfy (See, for example, [44,Theorem 40].) We call such w optimal weight vectors. Note that where ψ is an optimal Kantorovich potential. Given an optimal weight vector w ∈ R N , the unique optimal transport map from L 3 ¬ Ω to ν is given by In particular, if we define the function Φ : Ω → R by then T = ∇Φ and, using the notation introduced in Definition 2.1, This follows from the Gangbo-McCann Theorem (see, for example, [48,Theorem 1.17]). The Monge problem is a non-convex optimisation problem. As outlined above, in the case of semi-discrete optimal transport this can be replaced by an unconstrained, finite dimensional optimisation problem (maximising g), which is numerically tractable. As we demonstrate in Sect. 7, this is one motivation for using semi-discrete optimal transport to construct solutions of SG in geostrophic coordinates.

Optimal weight map
Optimal weight vectors are not unique. Indeed, let ν ∈ Q N (R 3 ) have mass vector m and seed vector z, and let e = (1, . . . , 1) ∈ R N . Using (11) it is easy to see that for any λ ∈ R, Hence, by the characterisation of optimal weight vectors (13), w ∈ R N is optimal if and only if every vector in w + span{e} is optimal. Conversely, if w, w ∈ R N and w − w / ∈ span{e}, then the pairs (z, w) and (z, w) define distinct Laguerre tessellations, so at least one of w and w is not optimal. In particular, it is easy to deduce that there is a unique optimal weight vector whose N th component is zero. This leads to the following definition, where e i denotes the i th canonical basis vector in R N . Definition 2.5 Given a fixed mass vector m ∈ R N , we define the optimal weight map w * : D → R N to be that which sends each seed vector z ∈ D to the unique optimal weight vector w * (z) ∈ span{e 1 , . . . , e N −1 }. We refer to the family {C i (z, w * (z))} N i=1 as the optimal Laguerre tessellation of Ω generated by z ∈ R N . Note that there are many possible definitions of the optimal weight map, which all yield the same definition of an optimal Laguerre tessellation. For instance, a natural choice of range space would be span{e} ⊥ . We choose the range space span{e 1 , . . . , e N −1 } so that subsequent arguments concerning the regularity of w * can be carried out using only the canonical basis of Euclidean space.

Statement of existence theorems
After stating some preliminary definitions, we state the two existence results (Theorems 3.5 and 3.6) for which we give novel proofs. See Definition 2.1 and Eq. (9) for the definitions of the operator B and the space Q N (R 3 ), respectively.
Note that the geostrophic energy is traditionally written in terms of the geostrophic velocity field u g = (u g,1 , u g,2 , 0) T and the potential temperature θ as

Remark 3.2
The geostrophic energy functional is continuous on P(K ) for any compact set K ⊂ R 3 . To see this first note that, for any ν ∈ P c (R 3 ), and, by continuity of B on P(K ) (Theorem 2.2), for all ϕ ∈ D(R 3 × R), where the matrix J is defined by (2). In what follows, we will refer to such a map α as a geostrophic solution. In particular, if there exists N ∈ N, a (k-times continuously differentiable) map z = (z 1 , . . . , z N ) : [0, T ] → R 3N and a mass vector m ∈ R N such that for all t ∈ [0, T ], we will refer to α as a (k-times continuously differentiable) discrete geostrophic solution. If the corresponding geostrophic energy E[α t ] is constant in time, then we say that α is energy-conserving. for all t ∈ [0, T ] then, since by a change of variables, Hence α is a distributional solution of the active transport equation This is the setting considered in [4].
is well-prepared for SG in geostrophic coordinates if there exists r > 0 such that In other words, β is well prepared if the seeds z i lie in distinct horizontal planes.
Well-preparedness of the initial data ensures that the seeds do not collide and therefore do not enter the set on which the right-hand side of the ODE in (4) is discontinuous. This allows us to prove the following theorem.

this solution is energy-conserving.
In Example 6.1, we show that Theorem 3.5 extends easily to the case N = 1 by deriving an explicit expression for the solution.
The uniqueness and regularity of solutions attained in Theorem 3.5 are significant in the context of numerical analysis since these factors determine the convergence of the corresponding numerical method. As mentioned in the introduction, the result of Theorem 3.5 is analogous to [22,Proposition 4.14], but with the following differences. In [22,Proposition 4.14], the physical domain Ω need not be convex, and the initial measure can be any convex combination of Dirac masses. However, with these slightly weaker hypotheses, the seed trajectories z are only known to be Lipschitz in time and are not known to be unique. For any T ∈ (0, ∞) and any α ∈ P c (R 3 ), there exists an energy-conserving geostrophic solution α ∈ C 0,1 ([0, T ]; P c (R 3 )) with initial measure α and a sequence (α N ) N ∈N of twice continuously differentiable discrete geostrophic solutions which converges uniformly in C([0, T ]; P c (R 3 )) to α: The existence of geostrophic solutions with arbitrary compactly supported initial measure was first proved by Loeper [40,Theorem 2.3] and later by Feldman and Tudorascu in the context of weak Lagrangian solutions (see [23,Theorem 3.2]). Being uniform in time, as opposed to pointwise in time, the convergence obtained in Theorem 3.6 is stronger than that obtained in these previous works. Since we work only with compactly supported measures, the spatial convergence obtained in Theorem 3.6 is equivalent to that obtained in [23]. As in [23], the limit point obtained in Theorem 3.5 is Lipschitz in time. Note that the corresponding theorems in [23,40] do not include the hypothesis that Ω is convex. We include this hypothesis for technical reasons relating to the regularity of the optimal centroid map (see Definition 4.1 and Remark 4.12).

Remark 3.7 (Conservation of transport cost)
From the proof of Lemma 4.14 it is easy to deduce that each discrete solution α N conserves not only the geostrophic energy but also the transport cost W 2 (α N t , L 3 ¬ Ω). Therefore any solution α constructed as the uniform limit of discrete solutions also conserves the transport cost W 2 (α t , L 3 ¬ Ω) by continuity of the Wasserstein distance. This is essentially due to the Hamiltonian structure of Eq. (1) (see [3,Example 8.1(c)]) and can be seen as a special case of [3,Theorem 5.2]. Conservation of the transport cost for weak Lagrangian solutions is proved in [22,Corollary 5.2]. For clarity, we include an independent proof based on our semi-discrete solution procedure.

Proof of Theorem 3.5
Fix N ∈ N, N ≥ 2. (The case N = 1 is discussed in Example 6.1.) To construct a twice continuously differentiable discrete geostrophic solution with well-prepared initial measure α ∈ Q N (R 3 ), we substitute the expression (18) into the transport equation (17) and, by appropriate choice of test functions, derive an ODE-IVP for the paths z = (z 1 , . . . , z N ). We then prove in Proposition 4.4 that this ODE-IVP has a unique C 2 -solution. Conversely, we show that any C 1 -solution of the ODE-IVP gives rise to an energy-conserving discrete geostrophic solution via the formula (18), from which Theorem 3.5 follows. The most involved part of the proof of Theorem 3.5 is the proof of the regularity of the optimal centroid map, which we now define.
where D is the set of seed vectors defined by (10). D is the set of generators of Laguerre tessellations of Ω with no zero-mass cells. We define the centroid map x : D → Ω N , x = (x 1 , . . . , x N ) by Moreover, given a fixed mass vector m, we define the optimal centroid map x = (x 1 , . . . , where w * is the optimal weight map (Definition 2.5).
We now characterise k-times continuously differentiable discrete geostrophic solutions in terms of solutions of an ODE-IVP involving the optimal centroid map.
is a k-times continuously differentiable discrete geostrophic solution with initial measure α, where z = (z 1 , . . . , z N ), and J N ∈ R 3N ×3N is the block diagonal matrix Proof Suppose that z = (z 1 , . . . , z N ) : [0, T ] → R 3N is a k-times continuously differentiable solution of (23) and let α be given by (22). Since z is continuous, α ∈ C([0, T ]; P c (R 3 )). We must check that α satisfies (17). Letting ϕ ∈ D(R 3 × R) be arbitrary, we have as required. Conversely, suppose that α : [0, T ] → Q N (R 3 ) given by (22) is a k-times continuously differentiable discrete geostrophic solution with initial measure α, in the sense of Definition 3.3. Then substitution of (22) into (17) for all ϕ ∈ D(R 3 × R). We now show that the paths z i can be separated, and deduce that the ODE-IVP (23) Hence For a fixed coordinate index k ∈ {1, 2, 3}, consider a test function By the Chain Rule, for all t ∈ I , Hence, with this particular choice of test function, Eq. (25) becomes Since φ ∈ C ∞ c (I ), the first integral is zero. By varying over all k ∈ {1, 2, 3} and all φ ∈ C ∞ c (I ) we obtain the pointwise statement thaṫ for all t ∈ I , in particular for t = t 0 . In a similar way, it can be shown that the inital condition z i (0) = z i is satisfied. Repeating this for all t 0 ∈ (0, T ), and all i ∈ {1, . . . , N }, we deduce that z = (z 1 , . . . , z N ) satisfies the ODE-IVP (23). (23) is precisely the discrete analogue of the active transport equation (1). Indeed, let ν ∈ Q N (R 3 ) have seed vector z = (z 1 , . . . , z N ) and consider the extension by +∞ of the convex function B[ν] to R 3 , which we will again denote by B
Having characterised k-times continuously differentiable discrete geostrophic solutions in terms of solutions of the ODE-IVP (23), we now aim to show that solutions of (23) exist and that the corresponding geostrophic solutions are energy-conserving. After proving some preliminary results (Lemma 4.13 and Lemma 4.10), we prove the following Proposition. z = (z 1 , . . . , z N
By Lemma 4.2, this gives rise to a unique discrete geostrophic solution with initial measure (21) via the formula (22). As we will see in Sect. 5, the condition (26) on the initial data is not restrictive for the purpose of constructing a geostrophic solution with arbitrary initial measure in P c (R 3 ). We begin by explaining its relevance.
We show below that the map W , defined by (24), is continuously differentiable, and therefore locally Lipschitz on D. However, as we demonstrate in Remark 4.11, W does not, in general, admit continuous extension to R 3N . In order to apply the Picard-Lindelöf Existence Theorem to obtain a unique solution of (23) on a given time interval [0, T ], it is therefore necessary to ensure that any solution trajectory is bounded away from the boundary of D. This is guaranteed a priori if z satisfies (26). Indeed, since the third row of the matrix J is zero, if z ∈ R 3N satisfies (26) and z is a solution of the corresponding ODE-IVP (23), then z(t) satisfies (26) for all times t.
We now prove the claimed regularity of the map W in the case where the domain Ω is open, bounded and convex (see Remark 4.12 for a discussion of the case where Ω is nonconvex). To do so, we use the following results regarding the regularity of the centroid map x (see Definition 4.1) and the volume map V (see Definition 4.8), as well as the structure of the matrix of partial derivatives of V with respect to w, which can be described using the notions of a Laplacian matrix and the dual graph of a Laguerre tesselation.

Definition 4.5 (Dual graph) Given a Laguerre tessellation {C
. . , N } the set of neighbours of i is defined to be the set The degree matrix of G, D G ∈ R N ×N , is the diagonal matrix such that The Laplacian matrix of G, L G ∈ R N ×N , is then given by f (x) dx.

Then F is continuously differentiable. In particular, for (z, w) ∈ D let G = (V , E) denote the dual graph of the Laguerre tessellation
Then the matrix D w F(z, w) of partial derivatives of F with respect to w evaluated at (z, w) is the Laplacian matrix of the weighted graph G = (V , E, h).
Note that the expressions for the partial derivatives of F with respect to the weights given above differs by a factor of 2 from those given in [15] due to our choice of quadratic cost function. Note also that the assumption that the seed locations are generic with respect to the cost, which is used in the proof of [44,Theorem 45], is not needed for the quadratic cost if the cells are non-empty. z, w)) .

Definition 4.8 We define the volume map
By combining Theorem 4.7 with the quotient rule for derivatives, we immediately obtain the following corollary.
where H 2 denotes the 2-dimensional Hausdorff measure on R 3 . Then the matrix D w V (z, w) of partial derivatives of F with respect to w evaluated at (z, w) is the Laplacian matrix of the weighted graph G = (V , E, h).
The optimal centroid map x (see Definition 4.1) is the composition of the centroid map x with the optimal weight map w * (see Definition 2.5). To show that x is continuously differentiable we now prove that w * is continuously differentiable. At the time of writing, we believe this result to be novel. Observe that A A T w * = w * so for all z ∈ D r(z, A T w * (z)) = 0. Now fix an arbitrary seed vector z ∈ D. Since m i > 0 for each i ∈ {1, . . . , N }, (z, w * (z)) ∈ D. By Corolloary 4.9, V is therefore continuously differentiable in a neighbourhood of (z, w * (z)). The residual map r is therefore continuously differentiable in a neighbourhood of (z, A T w * (z)) and The weighted Laplacian matrix D w V (z, w * (z)) has positive entries and corresponds to a connected graph. Hence it is symmetric, positive semi-definite and its kernel is span{e}, where e = (1, . . . , 1) ∈ R N (see for example [46,Section 2.4]). The symmetric matrix L = D w r(z, A T w * (z)) is therefore invertible. Indeed, if y ∈ R N −1 , By the Implicit Function Theorem (see, for example, [18, Theorem 10.2.1, p.270]) applied to r at the point (z, A T w * (z)), the function A T w * : D → R N −1 is therefore continuously differentiable in a neighbourhood of z. Since A A T w * = w * , it follows that the optimal weight map w * is continuously differentiable in a neighbourhood of z, as required.
Letting x * denote the centroid of the spherical cap it can easily be shown that the corresponding optimal centroids are given by does not admit a continuous extension to the whole interval (−1, 1). Since z is continuous on (−1, 1), this means that the optimal centroid map x does not admit a continuous extension at the point (0, 0) ∈ R 6 , precisely where the two seeds are the same. For s ∈ (−1/2, 1/2), the cell C 1 (z(s), w) is then given by and has volume z(s), w) is not differentiable at s = 0 which, since z is differentiable, implies that V 1 is not differentiable at the point (z(0), w).
If Ω is open, bounded and convex, we can use Corollary 4.9 to deduce that the optimal centroid map x is continuously differentiable and, therefore, locally Lipschitz. While it is known that the volume map V is Lipschitz in w even when the domain Ω is non-convex (see, for example [44,Proposition 41]), to prove that x is locally Lipschitz for non-convex domains would require a finer analysis of the regularity of the centroid map x and the volume map V with respect to z.
We now prove that for any finite time horizon T , solutions of the ODE-IVP (23) on the interval [0, T ] are bounded and have bounded first derivatives. Our immediate application of these a priori estimates is to prove that such solutions exist (Proposition 4.4). As we show in Sect. 5, these estimates also yield the necessary compactness of the corresponding sequence of measure-valued maps to pass to the limit as N → ∞.

Lemma 4.13 (A priori estimates)
Proof Since the matrix J is skew symmetric and has operator norm 1, for any i ∈ {1, . . . , N } we have on the interval [0, T ] is given by (Note that this is the unique solution unless z i = 0.) Applying a comparison lemma such as [49, Theorem 13.2] then establishes (27). Hence, , which establishes (28).
Using the regularity results and the a priori estimates established above, we now prove Proposition 4.4.

Proof of Proposition 4.4
Combining Corollary 4.9 with Lemma 4.10 we see that the optimal centroid map x is continuously differentiable on its domain D. The map W is, therefore, locally Lipschitz on D. Now, let z = (z 1 , . . . , z N ) ∈ D satisfy (26) By the Picard-Lindelöf Existence Theorem (see, for example, [50, Theorem 2.2]) there exists a unique C 1 -solution of (23) on the interval [0, T 0 ].
Since the third row of the matrix J is zero, by (26) Hence B r 2 (z(T 0 )) ⊂ D. By the a priori estimate (27) proved in Lemma 4.13, for all z ∈ B r 2 (z(T 0 )), Applying the Picard-Lindelöf Existence Theorem again we conclude that there exists a unique C 1 -solution of (23) on the interval [0, min{T , 2T 0 }]. Repeating this process, we obtain a unique C 1 -solution z of (23) on the interval [0, T ]. Since the function W is To conclude this section we prove Theorem 3.5. We begin by proving that any continuously differentiable discrete geostrophic solution is energy-conserving.

Lemma 4.14 Any continuously differentiable discrete geostrophic solution is energyconserving.
Proof Suppose that α ∈ C([0, T ]; P c (R 3 )) is a continuously differentiable discrete geostrophic solution with initial measure

By Lemma 4.2, α is then given by
where, for each i ∈ {1, . . . , N }, Moreover, letting Then, by Eq. (16), We now show that the function t → E[α t ] is differentiable on the interval (0, T ) and that its derivative is zero, from which it follows that α is energy-conserving. Trivially, the third term of (30) is constant in time. By (29), since the third row of the matrix J is zero, Therefore the time derivative of the second term of (30) is also zero on (0, T ). Now define the function ζ : so that the first term of (30) is precisely Ω ζ(x, t) dx. Observe that for any fixed s ∈ (0, T ), the function t → ζ(x, t) is continuously differentiable at s for L 3 -almost every x in Ω.
For x in the interior of the cell C i (t) By hypothesis, z(s) ∈ D, so continuity at s of the map t → ∂ t ζ(x, t) follows from Corollary 4.9 and Lemma 4.10. By Lemma 4.13, ∂ t ζ(·, s) is uniformly bounded on Ω. Hence, by the Mean Value Theorem and the Dominated Convergence Theorem, we may exchange the order of differentiation and integration and, using (29) and the skew-symmetry of the matrix J , we obtain which completes the proof.

Proof of Theorem 3.5
Since α is well-prepared in the sense of Definition 3.4, by Proposition 4.4 the corresponding ODE-IVP (23) has a unique C 2 -solution. By Lemma 4.2, this gives rise to a unique twice continuously differentiable discrete geostrophic solution with initial measure α via the formula (22). This solution is energy-conserving by Lemma 4.14.

Proof of Theorem 3.6
We now construct a geostrophic solution with arbitrary initial measure α ∈ P c (R 3 ). We begin by proving the existence of a sequence (α N ) N ∈N of well-prepared discrete probability measures (see Definition 3.4) converging to α with respect to the Wasserstein 2-distance (Lemma 5.1). The existence of a sequence of discrete geostrophic solutions α N with initial measures α N , respectively, is then guaranteed by Theorem 3.5. To obtain a geostrophic solution with initial measure α, we then apply the Arzelà-Ascoli Compactness Theorem combined with the continuity of B (Theorem 2.2) and pass to the limit as N → ∞.
There exists a compact set U ⊂ R 3 , a sequence of discrete probability measures β N ∈ Q N (R 3 ) given by and a sequence of positive real numbers r N > 0 such that and Proof Take a compact set U ⊂ R 3 such that supp(β) ⊂ U and a sequence of discrete probability measuresβ whose support is contained in U , such that lim N →∞ Such a sequence exists by for example [44,Lemma 10].
For each N ∈ N, let z N i ∈ U , i ∈ {1, . . . , N }, satisfy the following: Combining (33) and (34) completes the proof. Now fix T ∈ (0, ∞) and α ∈ P c (R 3 ). By Lemma 5.1 there exists a sequence of wellprepared discrete probability measures α N ∈ Q N (R 3 ), given by By Theorem 3.5, for each N ∈ N there exists a unique twice continuously differentiable discrete geostrophic solution α N ∈ C([0, T ]; P c (R 3 )) with initial measure α N 0 = α N , given by Proof The a priori estimates (27) on the paths z N i mean that the metric space (P c (R 3 ), W 2 ) can be replaced by the compact metric space (P(K ), W 2 ), where K ⊂ R 3 is the closed ball of radius R 1 := M + R(Ω)T centred at the origin. Since K is compact, W 1 and W 2 are equivalent metrics on P(K ) (see [48, p.179]), so it is sufficient to prove that (α N ) N ∈N has a uniformly convergent subsequence in C ([0, T ]; (P(K ), W 1 )) whose limit point is Lipschitz with respect to the W 1 metric. Moerover, since the arrival space (P(K ), W 1 ) is a compact metric space, by the Ascoli-Arzelá Theorem (see for example [48,p.10,Box 1.7]), such a sequence exists if and only if for every ε > 0, there exists δ(ε) > 0 such that, whenever s, t ∈ [0, T ] and |t − s| < δ(ε), W 1 (α N t , α N s ) < ε for all N . For μ, ν ∈ P(K ), (See for example [51,Theorem 1.14].) Let φ : K → R be 1-Lipschitz and let t, s ∈ [0, T ].
By the a priori estimate (28) we have Using the characterisation of the Wasserstein 1-distance given by (35) we obtain Choosing δ(ε) = ε/L, the Ascoli-Arzelá Theorem guarantees the existence of a uniformly convergent subsequence (α k(N ) ) N ∈N . Denote by α its limit point. Combining the Lipschitz estimate (36) and using the triangle inequality in (P(K ), W 1 ), for all N ∈ N we have Passing to the limit as N → ∞, we see that α is Lipschitz with respect to W 1 as required.
We conclude this section with the proof of Theorem 3.6.

Proof of Theorem 3.6
Let α ∈ C 0,1 ([0, T ]; P c (R 3 )) be the limit point of the sequence (α k(N ) ) N ∈N obtained in Lemma 5.2. First, we prove that α satisfies the transport equation (17). For clarity, given β ∈ C([0, T ]; P(K )) and ϕ ∈ D(R 3 × R) we define Since the space is a dense subspace of D(R 3 × R) (see [26]), to show that α satisfies (17) it is enough to check that F [α, ϕ] = 0 for any ϕ ∈ D. Moreover, for each N ∈ N we have By the triangle inequality, it is therefore sufficient to show that Since K is compact, W 1 and W 2 are equivalent metrics on K . The sequence (α k(N ) ) N ∈N therefore converges to α uniformly in C([0, T ]; (P(K ), W 1 )). By the characterisation of W 1 given by (35), this implies that for any Lipschitz function η : K → R, For Therefore, by (38), Moreover, since uniform convergence implies pointwise convergence, we have Finally, let By combining (39), (40) and (41), we see that (37) holds.
To complete the proof we note that each discrete solution α k(N ) is energy-conserving so, by continuity of the geostrophic energy functional E on P(K ) (see Remark 3.2), for any t ∈ [0, T ] which means that α is also energy-conserving.

Explicit solutions
Here we consider two special cases (Examples 6.1 and 6.2) where we can obtain explicit expressions for discrete geostrophic solutions. In Example 6.3 we show that L 3 ¬ Ω is an equilibrium solution of (1) and that its optimal quantisers are equilibrium solutions of the corresponding ODE-IVP (23).

Example 6.1 (A single mass)
First we consider the case of a single Dirac mass. This example has been discussed in both [21,Section 5] and [23,Section 2.2] in the context of the 2dimensional semi-geostrophic equations on the physical domain B 1 (0) ⊂ R 2 . In contrast with previous approaches, which use approximations of the Dirac mass by characteristic functions on balls, our solution follows immediately from the characterisation of discrete geostrophic solutions in terms of the ODE-IVP (23).
Let Ω ⊂ R 3 be open, bounded and convex, denote by x Ω the centroid of Ω, and let z ∈ R 3 . By an argument analogous to the proof of Lemma 4.2, a map α ∈ C([0, T ]; P c (R 3 )) given by is a k-times continuously differentiable discrete geostrophic solution with initial measure α = δ z if and only if z : [0, T ] → R 3 is a k-times continuously differentiable solution of the ODEż satisfying the initial condition z(0) = z. Such a map z is unique and is given by Moreover, Let z 3 = z · e 3 . By Eq. (16), the corresponding geostrophic energy satisfies be an optimal quantiser of L 3 ¬ Ω. By, e.g., [19,Proposition 3.1], [29,Corollary 4.3], the seeds z = (z 1 , . . . , z N ) generate a centroidal Voronoi tessellation of Ω, i.e, This means that W (z) = 0 so z is an equilibrium solution of the ODE-IVP (23). The behaviour of L 3 ¬ Ω under the dynamics of Eq. (1) therefore agrees exactly with that of its optimal discrete approximants.

Numerical simulations
One advantage of the constructive existence proof given in this paper is that it naturally leads to a numerical method (a meshfree method) and moreover it tells us something about the convergence of this method. To be precise, given an initial measure α ∈ P c (R 3 ), we can construct a numerical approximation of a solution of the semi-geostrophic equations (17) as follows: 1. Approximate α by a discrete measure α N = N i=1 m i δ z i . This leads to the semi-discrete numerical scheme (continuous in time, discrete in space) used in the proof of Theorem 3.6, where an exact solution α N = N i=1 m i δ z i (t) of (17) with initial condition α N is constructed by solving the system of ODEsż = W (z) given in Eq. (23). To turn this into a bona fide numerical method we require a further discretisation: 2. Use a time-stepping scheme to approximately solve the ODEż = W (z). Every evaluation of the vector field W requires a further numerical approximation; to evaluate W (z(t)) we must solve the semi-discrete transport problem W 2 (L 3 ¬ Ω, α t ) in order to compute the centroids x(z(t)).
We demonstrate the viability of this numerical method by giving an example in two dimensions, implemented in MATLAB. Due to a lack of space we postpone the three-dimensional implementation and a more thorough numerical study to a further paper. Before giving the example we briefly discuss some implementation issues and convergence.

ODE solver
We used the explicit Runge-Kutta scheme RK4 [35,Example 5.13] to solve the ODEż = W (z). For larger simulations it may be better to use a linear multistep method [35,Section 5.9] since each evaluation of the vector field W is expensive (linear multistep methods only require one new vector field evaluation per time step, whereas RK4 requires four per time step).
Semi-discrete transport solver The semi-discrete transport problem W 2 (L 3 ¬ Ω, α t ) can be solved by maximising the concave function g (see Eq. (12)) as described in Sect. 2. We did this in MATLAB using a quasi-Newton method (the MATLAB function fminunc). Every evaluation of g and its gradient requires a Laguerre diagram to be computed, which we did using the MATLAB function power_bounded [25]. This simple method, which is described in more detail in [6, Algorithm 1 and Section 4], was sufficient for our proof of concept simulations here. In general, however, to maximise g it would be much faster to use the damped Newton method from [32], especially for large 3D simulations.
Convergence By Theorem 3.6 the sequence (α N ) N ∈N generated by the semi-discrete scheme has a subsequence that converges (uniformly with respect to the Wasserstein metric) to a solution of the semi-geostrophic equations (17) with initial condition α. In particular, if Eq. (17) has a unique weak solution with initial measure α, then the whole sequence of approximations α N converges to the true solution. Local-in-time uniqueness of Hölder continuous periodic solutions of (17) was proved in [40], but is not known in general. By the conservation of the transport cost W 2 (L 3 ¬ Ω, α t ) (see Remark 3.7) it is easy to see that the whole sequence (α N ) N ∈N also converges in the very special case where the initial measure α is the Lebesgue measure on Ω. (Note that the Lebesgue measure is an equilibrium solution of (17) as discussed in Example 6.3.) In general, proving convergence of the whole sequence (α N ) N ∈N is beyond the scope of this paper, as is proving convergence of the fully discrete scheme. We will study these in a future paper. Note that we have dropped the previous normalisation convention that α(R 2 ) = L 2 (Ω) = 1, which is not necessary for the analysis; we simply require that α(R 2 ) = L 2 (Ω). We approximated α by a discrete measure of the form α N = N i=1 m i δ z i with N = 2000 seeds using 1000 iterations of Lloyd's algorithm [19]. We approximated the solution of the ODE-IVP (23)  Note that all the seeds start off in Ω but they are not confined there.
As an accuracy check, we repeated these simulations with a smaller time step size of h = 0.005 and a finer optimal transport tolerance of ε = 0.05. We found that the x-and y-coordinates of the seeds z i at the final time step T = 5 agreed with our previous results to within 10 −3 . Recall from Remark 3.7 that the exact dynamics (23) preserves the transport cost: In our simulations (with h = 0.01, ε = 0.1) the transport cost was conserved to within 7.5 × 10 −7 : max n W 2 2 (L 2 ¬ Ω, α t n ) − 1 501 500 m=0 W 2 2 (L 2 ¬ Ω, α t m ) < 7.5 × 10 −7 where t n = nh, n ∈ {0, 1, . . . , 500}.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.