A new implementation of the geometric method for solving the Eady slice equations

We present a new implementation of the geometric method of Cullen&Purser (1984) for solving the semi-geostrophic Eady slice equations which model large scale atmospheric flows and frontogenesis. The geometric method is a Lagrangian discretisation, where the PDE is approximated by a particle system. An important property of the discretisation is that it is energy conserving. We restate the geometric method in the language of semi-discrete optimal transport theory and exploit this to develop a fast implementation that combines the latest results from numerical optimal transport theory with a novel adaptive time-stepping scheme. Our results enable a controlled comparison between the Eady-Boussinesq vertical slice equations and their semi-geostrophic approximation. We provide further evidence that weak solutions of the Eady-Boussinesq vertical slice equations converge to weak solutions of the semi-geostrophic Eady slice equations as the Rossby number tends to zero.


Introduction
In this paper, we study a meshfree method for solving the semi-geostrophic Eady slice equations, which are stated in equations (8), (9) in Eulerian coordinates, and in equation (16) in Lagrangian coordinates. The Eady slice model is a simplified model of large-scale (small Rossby number) atmospheric flow, and is capable of predicting the formation of atmospheric fronts [9]. We develop an efficient implementation of the geometric method [10] for solving this PDE using an adaptive time-stepping algorithm (Algorithm 2). As an application, we use it to test the validity of the semi-geostrophic approximation of the Euler-Boussinesq equations.
The diagnostics extracted in [36,Section 3.3] show that they systematically violate energy conservation properties after front formation, resulting in a significant loss of total energy over multiple lifecycles.
The main application of our new algorithm is to obtain numerical evidence that the semigeostrophic Eady slice equations (8), (9), in which geostrophic balance of the out-of-slice wind and hydrostatic balance are enforced, are the correct limit of the Eady-Boussinesq vertical slice equations (7) in the large-scale limit (i.e. as the Rossby number Ro tends to 0). In the classical setting, solutions of the Eady-Boussinesq vertical slice equations with periodic boundary conditions converge strongly in L 2 to solutions of the corresponding semi-geostrophic equations as Ro → 0 [7,Theorem 4.5]. Correspondingly, the numerical solutions obtained in [34] satisfy geostrophic and hydrostatic balance, up to an error of order O(Ro 2 ), while the solutions are smooth, and order O(Ro) after front formation. On the other hand, in both [34] and [36] it was remarked that there is a significant difference between the Eulerian solutions presented in those papers and the semigeostrophic solution presented in [11], which cannot be accounted for by the dissipative nature of the former. Subsequent investigations revealed that the problems studied in the two papers are not identical. Our implementation of the geometric method enables a controlled comparison between the two approaches. Using the physical parameters from [34], we obtain numerical solutions of the semi-geostrophic Eady slice equations (8), (9) and confirm that the differences between these and the Eulerian solutions obtained in [34] are consistent with the loss of Lagrangian conservation in the Eulerian solutions. This provides numerical evidence that the semi-geostrophic approximation is the correct small Rossby number approximation of the Eady-Boussinesq vertical slice equations, and suggests that a result similar to [7,Theorem 4.5] may also hold for weak solutions which, unlike classical solutions, can possess non-dissipative singularities that represent the evolution of weather fronts.
Outline of the paper. In Section 2 we introduce the semi-geostrophic Eady slice equations in Eulerian coordinates (Section 2.1) and Lagrangian coordinates (Section 2.3), as well as an important steady shear flow. In Section 3 we discretise the Eady slice equations using the geometric method. This leads us to the system of ODEs (32). We give an algorithm for solving these ODEs in Section 4, which involves discretising the initial data (Section 4.1), evaluating the right-hand side of the ODE by solving a semi-discrete optimal transport problem (Section 4.2), and using an adaptive timestepping scheme (Section 4.4). Section 5 includes some numerical experiments, where we study the stability of a steady shear flow. We illustrate front formation in Section 5.2 and study the validity of the semi-geostrophic approximation in Section 5.4. This numerical study is complemented by an analytical linear instability analysis in Appendix C.

Governing equations
We study the semi-geostrophic approximation of the Eady-Boussinesq vertical slice model in an infinite channel of height H. We will consider solutions that are 2L-periodic in the first coordinate direction. We start by summarising the derivation of this model from the Euler-Boussinesq equations and we then draw comparisons with the Eady-Boussinesq vertical slice model considered in [34].

Eady slice model
Our starting point is the Euler-Boussinesq equations for an incompressible fluid: ∇ · u = 0.
Here s < 0 is a constant, soθ represents the decrease in potential temperature when moving away from the equator and towards the north pole. Observe that φ 0 +φ is in hydrostatic balance with θ 0 +θ, which means that ∂ z (φ 0 +φ) = g θ 0 (θ 0 +θ).
We seek vertical slice solutions of (1)- (3), where the velocity depends only on x, z and t. We consider the fluid equations in the (x, z)-domain we impose 2L-periodic boundary conditions in x, and we require that w(x, z, t) = 0 for z ∈ {−H/2, H/2}. Substituting u = u(x, z, t), (4), and (5) into (1)-(3) yields the following form of Eady-Boussinesq vertical slice equations: where D t := ∂ t + u∂ x + w∂ z is the in-slice material derivative operator.
The semi-geostrophic approximation is obtained from (7) by neglecting the derivatives D t u and D t w. This results in the system Using (9) we obtain the relation . (12) The first component of ∇P is commonly referred to as the absolute momentum. The semigeostrophic (SG) Eady slice equations (8) then become where id x (x, t) = x, e 1 := (1, 0), and The global-in-time existence of solutions of (13) and its 3-dimensional analogue for physicallyrelevant initial data remains an open problem. For future reference, we record that the steady state (11) corresponds to We now introduce the Lagrangian form of the equations. Let F denote the in-slice flow corresponding to u so that F(x, t) ∈ Ω is the position at time t of a particle that started at position x, that is Since u(·, t) is divergence free, the flow F is mass-preserving. Written in terms of the Lagrangian variable the transport equation (13) becomes In contrast to the Eulerian setting, global-in-time solutions of the equations in Lagrangian coordinates are known to exist for a wide and physically-relevant class of initial data [17].

Energy
We define the total geostrophic energy at time t to be where is the total kinetic energy due to v and is the total potential energy. (The final integral in the definition of P(t) ensures that the total potential energy of the steady state (11) is zero. We include it for consistency with [36].) The total geostrophic energy E is conserved by the Eulerian equations (13). Note that in contrast to the energy defined in [36, equations (23)- (26)], the total geostrophic energy E does not include the kinetic energy due to the in-slice velocity u.

The geometric method
The geometric method is a spatial discretisation of equation (16), which yields a system of ordinary differential equations. It is derived using an energy minimisation principle known as Cullen's stability principle. It was first described and implemented in [10], and later used in the context of the SG Eady slice problem in [11]. We describe the original implementation in Section 4.5. In this section we recast the geometric method in the language of semi-discrete optimal transport theory, which underpins our novel implementation and will aid its description.

The stability principle and semi-discrete optimal transport
In the present context, the stability principle can be stated as follows: Stable solutions of (16) are those which, at each time, minimise the total geostrophic energy (17) over all periodic mass-preserving rearrangements of fluid particles that conserve the absolute momentum and potential temperature.
This can be shown to be equivalent to assuming that P is convex.
In the geometric method, we approximate the modified geopotential P at each time by a piecewise affine function, and apply the stability principle. The gradient of any piecewise affine function on the fluid domain Ω is uniquely identified by a tessellation of Ω by cells S i , i ∈ {1, . . . , n}, and corresponding points z i ∈ R 2 , where z i is the gradient of the piecewise affine function on S i . Each periodic rearrangement of fluid particles corresponds to a different tessellation of Ω by cells S i . Such a rearrangement is mass-preserving if S i has the same area as S i for all i ∈ {1, . . . , n}. The absolute momentum and potential temperature are conserved when the corresponding image points z i are fixed.
By writing the geostrophic energy (17) in terms of the points z i and sets S i , accounting for the periodic boundary conditions on φ, and neglecting terms that are constant over all periodic mass-preserving rearrangements, the stability principle can be rephrased as a semi-discrete optimal transport problem (Definition 3.1); see Appendix A for the derivation. This is a type of optimal partitioning problem. We refer to its solution as an optimal partition.
In what follows, for A ⊂ R 2 we denote by |A| the area of A.
Definition 3.1 (Optimal partition). Given seeds z = (z 1 , . . . , z n ) ∈ R 2n with z i = z j if i = j and target masses m = (m 1 , . . . , m n )∈ R n with m i > 0 and n i=1 m i = |Ω|, a partition of Ω is said to be optimal if it minimises the transport cost of Ω that satisfy the mass constraint Here |x − y| per is the distance between x, y ∈ R 2 taking into account 2L-periodicity in the first component: Recall that the characteristic function 1 A : Any modified geopotential P that is piecewise affine at time t and satisfies the stability principle has the form where {S i (t)} n i=1 is an optimal partition corresponding to the seeds z i (t) in the sense of Definition 3.1, and accounts for the periodic boundary condition on the geopotential φ (cf. [33,Theorem 1.25]). Note that k * (x, z i ) is well defined for almost-every x ∈ Ω. Moreover, without loss of generality, Optimal partitions can be described in terms of periodic Laguerre diagrams. These are partitions of the domain Ω into cells parametrised by the seeds z = (z 1 , . . . , z n ) and a set of weights w = (w 1 , . . . , w n ) ∈ R n .
This is the i th periodic Laguerre cell generated by (z, w), and the collection of all cells is the periodic Laguerre diagram generated by (z, w). It is well known that given z = (z 1 , . . . , z n ) ∈ R 2n , with z i = z j whenever i = j, and given m = (m 1 , . . . , m n ) ∈ R n with m i > 0, n i=1 m i = |Ω|, there exists a unique weight vector w * (z) ∈ R n with final entry 0 that satisfies the mass constraint m(z, w * (z)) = m (see for example [28,Corollary 39]). The optimal partition of Ω is then the periodic Laguerre diagram generated by (z, w * (z)). We call w * (z) the optimal weight vector for the seeds z. Observe that, for all λ ∈ R, where e = (1, . . . , 1) ∈ R n . Choosing the last entry of w * (z) to be 0 ensures uniqueness of the optimal weight vector, without loss of generality.
To solve the semi-discrete optimal transport problem (Definition 3.1), it is therefore sufficient to compute the optimal weight vector. To do this, we use the well-known fact that w * (z) ∈ R n is the unique maximum (with final entry 0) of the Kantorovich functional, which is the concave function K : R n → R defined by See for example [28,Theorem 40]. In effect, the constrained optimal partitioning problem is transformed into an unconstrained, finite-dimensional, concave maximisation problem, which is numerically tractable.

Spatial discretisation
In summary, we seek solutions of (16) for which the associated geopotential P is piecewise affine in space at each time t. By the stability principle, ∇P must have the form (22) with for some time-dependent map z = (z 1 , . . . , z n ), where w * is the maximum of K and the chosen target masses m = (m 1 , . . . , m n ) do not change over time.
The assumption that S i has the form (27) is equivalent to assuming that P is convex. We now sketch the derivation of an ODE for z, further details of which can be found in Appendix B. We begin by making the following definition, which arises naturally in the derivation due to the periodic boundary conditions. Definition 3.4. Let z = (z 1 , . . . , z n ) ∈ R 2n with z i = z j if i = j. Let w = (w 1 , . . . , w n ) ∈ R n . For i ∈ {1, . . . , n}, we define the convex polygon C i (z, w) by Note the difference between Definitions 3.2 and 3.4: Definition 3.2 defines a periodic Laguerre tessellation of [−L, L) × [−H/2, H/2] generated by a finite set of seeds and weights (z, w), while Definition 3.4 defines a non-periodic Laguerre tessellation of R × [−H/2, H/2] generated by all periodic copies of (z, w), namely by (z i + k, w i ) for all i ∈ {1, . . . , n}, k ∈ K. In general, the cells C i (z, w) are not contained in Ω, but if C i (z, w) ⊂ Ω then C i (z, w) = C i,per (z, w). In any case, both C i (z, w) and C i,per (z, w) have the same area for each i ∈ {1, . . . , n}.
By definition of the optimal weight vector, at each time t the mass constraint m z(t), w * z(t) = m is satisfied. For brevity, we define Assume that there exists a mass-preserving flow F such that at all times t. For each i ∈ {1, . . . , n}, we define c i (z) to be the centroid of C i (z, w * (z)), namely and we define c(z) := (c 1 (z), . . . , c n (z)).
Formally, after substituting the ansatz (22), (27) into (16), integrating over the cell C i,per (0), and using (29) to apply the change of variables x → F(x, t), one finds that the vector of seed trajectories z satisfies the ODE for all i ∈ {1, . . . , n}, where z = (z 1 , . . . , z n ) denotes the initial seed positions: see Appendix B for the details of the derivation. The ambient space containing individual seed positions z i (t) is referred to as geostrophic space. The ODE (31) is the discrete analogue of the Lagrangian equation (16), where the centroid map c plays the role of the flow F. The full system can be written compactly as where J and Π are the 2n × 2n matrices J = diag(J, . . . , J) and Π = diag(e 1 ⊗ e 1 , ..., e 1 ⊗ e 1 ).
The matrix Π acts on each seed z i by projection onto the horizontal coordinate. We will use this form of the system when we describe the time discretisation of the system: see Algorithm 2 below.

Structure preservation and recovery of physical variables
If z satisfies the ODE (32) and w(t) = w * (z(t)), then the total geostrophic energy (17) corresponding to a modified geopotential P of the form (33) is constant in time. The geometric method therefore inherits the energy conservation property possessed by the Eulerian SG Eady slice equation (13). Moreover, the one-parameter family of Laguerre tessellations t → {C i (z(t), w(t))} i∈{1,...,n} corresponds to an area-preserving flow of the fluid since the mass of each cell is conserved by definition of w * . Given seed trajectories z = (z 1 , . . . , z n ) satisfying the ODE (32), let w(t) = w * (z(t)). The corresponding modified geopotential P is and k * (x, z i ) was defined in (23). This is a piecewise affine convex function whose gradient is defined almost-everywhere and is given by Expressions for the approximate out-of-slice velocity and potential temperature can then be recovered using (12).

Implementation
We now give details of our numerical implementation of the geometric method. Our software is available at the following GitHub repositories.

Generating discrete initial data
For fixed n ∈ N, we approximate given initial data ∇P (x, 0) by a piecewise constant function ∇P n 0 of the form (34) defined by initial seeds z = (z 1 , . . . , z N ) and cell areas m = (m 1 , . . . , m n ). We now describe how we choose z and m. In Section 5 we will take ∇P (x, 0) to be a small perturbation of the steady state ∇P , which was defined in (14). Therefore the seeds z i are taken from the image domain ∇P (Ω, 0), which is a small perturbation of the rectangle We define z i as follows: 1. use Lloyd's algorithm [13] to generate points y = (y 1 , . . . , y n ) that are approximately uniformly distributed in R; 2. for each i ∈ {1, . . . , n}, map y i into Ω by inverting ∇P : We briefly describe Lloyd's algorithm, which is an iterative method for quantising measures (see, for example, [13, Section 5.2]). Let y = (y 1 , . . . , For each i ∈ {1, . . . , n}, define the i th periodic Voronoi cell generated by y as The collection of all cells {V i,per (y)} n i=1 is called the periodic Voronoi tessellation of R. At each iteration of the Lloyd algorithm, the value of y i is updated by moving it to the centroid of its Voronoi cell: The centroids can be computed exactly (without numerical integration): see for example [6, online supplementary material, equation (4)]. In our implementation of Lloyd's algorithm, we start with points y = (y 1 , . . . , y n ) on a regular triangular lattice in R and perform 100 iterations to obtain y = (y 1 , . . . , y n ). The target areas are defined for each i ∈ {1, . . . , n} as The method above generates points z i that are approximately optimally sampled from the distribution on ∇P (Ω, 0) that is obtained by pushing forward the uniform distribution on Ω by ∇P (·, 0). It is numerically cheaper than directly sampling this distribution using Lloyd's algorithm since it avoids numerical integration.
Since the discrete dynamics (32) conserves the total geostrophic energy E, it is desirable that the initial condition starts on the correct energy surface, i.e. that the discrete initial data ∇P n 0 has the same total geostrophic energy as ∇P (x, 0). While there exists discrete initial data with this property [15], the choice above does not. Nevertheless, it is easy to generate and it approximates the initial energy well enough if n is sufficiently large.

Generating Laguerre tessellations
In two dimensions the worst-case complexity of computing a Laguerre tessellation with n seeds is O(n log n), where n is the number of seeds. This can be achieved for example by the lifting method of Aurenhammer [2]. We computed Laguerre and Voronoi tessellations using the C++ library Voro++ [32] and our own mex file to interface with MATLAB. While this does not achieve the optimal scaling in n, it is sufficiently fast for the values of n that we use.

Solving the semi-discrete optimal transport problem
In any numerical scheme for solving the ODE (32), it is necessary to evaluate the right-hand side at each time step. Given seeds z, this involves solving the corresponding semi-discrete optimal transport problem to compute c(z). This is the most expensive part of the algorithm. Here we do this by finding the maximum of the Kantorovich functional K (equation (26)) using the damped Newton method developed in [23]; see Algorithm 1. This solves the nonlinear algebraic equation First, a guess w for the optimal weight vector is proposed. The Newton direction d is then determined by solving the linear system where e n = (0, . . . , 0, 1) ∈ R n . (Formulas for D 2 K and ∇K are given in Appendix D.) For the linear system (36) to have a unique solution, it is necessary and sufficient for w to satisfy the mass-positivity condition for all i ∈ {1, . . . , n}. To ensure that this condition is met by the subsequent iterate, backtracking is used to determine the length of the Newton step. The algorithm terminates once |∇K| is less than a given tolerance. Convergence of the damped Newton algorithm is guaranteed if and only if the initial guess for the weights satisfies (37), in which case it converges globally with linear speed, and locally with quadratic speed, as the number of iterations diverges (see [23, Proposition 6.1]). Our adaptive time-stepping scheme (Algorithm 2) for solving the ODE (32) provides a robust way to generate a good initial guess for the weights at each time step given the optimal weights from the previous time step.
It remains to generate a good initial guess for the weights at time t = 0. By [25, Theorems 3 and 4], each cell C i,per (z, w) is non-empty if and only if w is c-concave in that sense that there exists ϕ : Ω → R such that This condition ensures that the cells C i,per (z, w) are non-empty but not necessarily that they satisfy the mass-positivity condition (37). If ϕ = 0, then w defined by (38) satisfies (37) if in addition the horizontal components of the seeds z i are distinct (or if z i ∈ Ω for all i, but this is not the case for our simulations in Section 5). At time t = 0, we therefore first randomly perturb the horizontal components of the seeds: where ζ i ∈ R are appropriately scaled, randomly generated numbers. Then we apply Algorithm 1 with seeds z i and initial weights to obtain w * ( z). Since w * is continuous, for sufficiently small perturbations ζ i , w * ( z) is a good initial guess for w * (z) (which is computed in the initial step of Algorithm 2). Note that (39) is the square of the K-periodic distance from z i to Ω which can be simply computed as where e 2 = (0, 1).

Algorithm 1 Damped Newton algorithm of Kitagawa, Mérigot and Thibert [23]
Input: Seeds z = (z 1 , . . . , z n ), target masses m = (m 1 , . . . , m n ), a percentage mass tolerance η, and an initial guess for the weights w such that Initialisation: Set w 0 = w and k = 0. Convert η to an absolute mass tolerance: Step 1: Solve the following linear system for the Newton direction d (k) : Step 2: Determine the length of the Newton step using backtracking: find the minimum ∈ N∪{0} such that w (k, Step 3: Define the damped Newton update w (k+1) := w (k) + 2 − d (k) and k ← k + 1.

Solving the ODE using an adaptive time-stepping method
The most computationally intensive part of solving the ODE (32) numerically is evaluating the function c. In addition, c is continuously differentiable (on the set of distinct seed positions), but not twice continuously differentiable in general: see [4]. As such, we can only expect convergence of the ODE solver up to second order in the time step size. We therefore use an Adams-Bashforth 2-step method (AB2) with adaptive time-stepping. This requires only one function evaluation at each time step. Our adaptive time-stepping scheme (Algorithm 2) speeds up the function evaluations by using the ODE (32) to generate a good initial guess for the weights for Algorithm 1, as we now describe. Suppose that at a given time step we have seeds z and optimal weights w. Following the notation of Algorithm 2, for a proposed time step size h l , the AB2 scheme determines an increment z l inc for the seed positions. The seeds at the subsequent time step are then defined as We generate an initial guess w l for the weights corresponding to the seeds z l by taking a first order Taylor expansion of w * (z l ) around z: If the seeds and weights (z l , w l ) generate a periodic Laguerre tessellation of Ω with no zero-area cells, then the proposed step size is accepted. Otherwise, the seeds and weights (z l , w l ) are discarded, the proposed step size is halved and the updates are recalculated. This ensures the convergence of the damped Newton algorithm when used at the subsequent time step.
From numerical experiments, we observed that this method is up to forty times faster than either of the following methods for generating the initial guess for the weights: (i) using the weights from the previous time step, and a small enough time step to ensure that the cells have positive area; (ii) using (39) with a suitable perturbation of the seeds, without using any information from the previous time step or adaptive time stepping. Moreover, the sensitivity analysis in Section 5.5 shows that the adaptively chosen times step sizes are suitable, i.e. not unreasonably small.
The expression for D z w * (z) can be obtained by implicit differentiation of the mass constraint. Indeed, by definition of the optimal weight map w * , m(z, w * (z)) = m.
By differentiating (41) with respect to z, we see that the derivative D z w * (z) satisfies the equation The matrix D w m(z, w * (z)) is symmetric and singular with kernel spanned by e := (1, . . . , 1) ∈ R n . Each column of D z m(z, w * (z)) is orthogonal to e. Hence, the linear system (42) has a solution D z w * (z). We choose the solution where 0 ∈ R 2n is the zero vector and the matrices A ∈ R (n−1)×(n−1) and B ∈ R (n−1)×2n are defined by Algorithm 2 Adaptive time-stepping scheme Input: Initial seeds z and masses m, a final time T , a default time-step size h def (seconds), and a percentage mass tolerance η. Initial step: Compute the optimal weight vector w * (z) using Algorithm 1, and do a forward Euler step to determine seed positions at time h def . Set While: t ≤ T : Step 1: Compute the optimal weight vector w * (z) using Algorithm 1 with w as the initial guess for the weights. Set Step 2: Determine the minimum l ∈ N ∪ {0} such that min i∈{1,...,N } (m i (z l , w l )) > 0, where the updated seeds z l and weights w l are defined for l ∈ N as follows: Propose step size: Compute increments of seeds and weghts: Update seeds and guess for weights: Step 3: Set Output: Seed positions at each time step.
for i, j ∈ {1, . . . , n − 1} and for k ∈ {1, . . . , n − 1} and l ∈ {1, . . . , n}. Expressions for the derivatives of m with respect to z and w in the non-periodic setting are given for example in [12] and [6,Lemma 2.4]. In Appendix D, we state the analogous expressions in the periodic setting.

Comparison with the original implementation
In the original presentation of the geometric method [10], a convex piecewise affine modified geopotential P is constructed directly at each time t. The projection of its graph onto the (x 1 , x 2 )-plane gives a tessellation {C i } n i=1 of Ω. While the term 'Laguerre diagram' is not used in [10], this tessellation is equivalent to the optimal partition described above in terms of Laguerre cells (up to the inclusion of periodic boundary conditions). The existence and uniqueness of such a map P was also proved formally.
As described in [9, Section 5.1.2], for some unknown scalars p j , j ∈ {1, . . . , n}, each cell C i consists of all x ∈ Ω such that for all j ∈ {1, . . . , n}. The unknowns p j play the role of the weights used in Definition 3.2. (In the non-periodic case, the weights w i are related to the p i by p i = 1 2 w i − 1 2 |z i | 2 , and the modified geopotential P is given by P (x) = max i (x · z i + p i ).) The task is to find (p 1 , . . . , p n ) such that C i has given area m i for each i ∈ {1, . . . , n}. This is achieved by applying the nonlinear conjugate gradient method to the objective function At the first iteration, p i are defined such that all cells have positive area (see [9,Equation 5.10]). This initial guess is comparable to the initial guess for the weights used at the first time step in our algorithm. At each iteration, the cells C i are constructed using the divide-and-conquer algorithm of Preparata and Hong [31].

Results
In this section we use the geometric method to simulate a shear instability and the formation of an atmospheric front. In what follows, the root mean square of the meridional velocity v (RMSv) at time t is given by where |Ω| is the area of Ω.

Parameters and initial conditions
In all simulations we used the following physical parameters, taken from [34]: We considered four initial conditions ∇P (x, 0) of the form where G is a small perturbation of the steady shear flow (14). The perturbations G are related to perturbations θ of the potential temperature and v of the meridional velocity by In each case, a corresponding domain height H was chosen, either in line with the linear instability analysis, or to enable comparison with the literature: see Table 1.
The perturbations G u and G s , listed in Table 1 and defined below, are normal modes of equation (13) linearised around the steady state (14). Full details of the linearisation are contained in Appendix C. Define the Burger number by As shown in Section C.3, exponentially growing normal modes exist only when the Burger number is less than a critical value Bu crit , otherwise all normal modes are oscillatory. The domain height H in Table 1 corresponding to G u is chosen so that Bu < Bu crit , the maximum growth rate is achieved by the lowest frequency normal mode G u , and all other normal modes are either decaying or oscillatory. Conversely, the domain height in Table 1 corresponding to G s is chosen so that Bu > Bu crit , and the lowest frequency normal mode G s is oscillatory and has a wave speed of one domain length (i.e. 2L) every 16 days. Define the constants where and σ : R → R is given by The unstable normal mode G u is where and where The perturbation used in [34] and [36] is Finally, where B = 0.25 K.
These, and the corresponding parameter values in Table 1, are taken directly from the stated sources. Note that neither G V nor G C are normal modes of (13) linearised around the steady state (14). They do, however, have the same horizontal wavelength as θ u . Finally, in each simulation we specify three discretisation parameters: n, the number of cells in the spatial discretisation of Ω; h def , the default time step size (see Algorithm 2); η, the percentage mass tolerance (see Algorithm 1).

Unstable normal mode
The results reported in this subsection were obtained using the initial data given in Table 1, Row 1, and the simulation parameters η = 0.01, h def = 30 and n = 2678, unless otherwise stated. The dimensional growth rate of the unstable mode (53), and the corresponding RMSv, under the dynamics of the linearised equations (78)-(82) is This is derived in the appendix: see equation (102).
The time-evolution of the RMSv calculated from the simulation data is presented in Figure 1a. Using equations (12) and (34), the RMSv at time t corresponding to a solution z = (z 1 , . . . , z n ) of (32) is given by where k * (x, z i (t)) is defined by (23), and C i,per (t) is defined by (28). After a short decline, the RMSv grows at the predicted rate ω: see Figure 1b. From around time t = 5 days, this growth slows until a peak RMSv value is reached at time t = 7.5 days. After a period of decline, the RMSv reaches a trough and proceeds to oscillate between the peak and trough values with a frequency of around 7 days. The initial decline in the RMSv is due to numerical errors incurred by the spatial discretisation: see Section 5.5.
The initial peak in the RMSv curve at time t = 7.5573 days coincides with the formation of a frontal discontinuity, as seen in Figure 2. Subsequent peaks/troughs of the RMSv occur when the frontal discontinuity is strongest/weakest, respectively. Note that the meridional velocity v displayed in Figure 2 is computed from the numerical solution z using (12) and (34). It is a piecewise affine function with respect to the Laguerre cells C i,per , which are shown in Figure 2, Rows 1 and 3. They are coloured according to the meridional velocity at their centroids.
Up to first onset of frontogenesis, the distribution of seeds in geostrophic space appears to roughly approximate a two-dimensional subset of R 2 : see Figure 2, t = 4.7031 days. At the onset of frontogenesis, however, it is clear that a small subset of the seeds is distributed along a onedimensional curve: see Figure 2, t = 7.5573 days. In other words, the frontal discontinuity in physical space appears to correspond to a singular part of the potential vorticity 1 measure α, which is defined at time t to be the push-forward measure α t = ∇P (·, t)#L 2 ¬ Ω. For the case where α is non-singular (has a two-dimensional support), it is defined by The time evolution of the total geostrophic energy (17), the total kinetic energy (18) and the total potential energy (19) is shown in Figure 3a  the total geostrophic energy is conserved by a numerical solution, we define the energy conservation error at time t as where E n (t) is the total geostrophic energy at time t calculated from the numerical solution, E n is the temporal mean of this quantity, and n is the number of seeds used in the simulation. For the values of n listed in Table 2, |ε n (t)| < 2 × 10 −5 for all t. Figure 3b demonstrates that |ε n (t)| is at a local maximum at times t when the RMSv is at a peak or trough. By definition, peaks and troughs of the kinetic energy align with those of the RMSv. Since the total energy is conserved, the total potential energy has the same behaviour as the total kinetic energy but with opposite sign, and energy is transferred from kinetic to potential, and vice-versa, over multiple lifecycles. Numerical solutions of the Eady-Boussinesq vertical slice equations (7) obtained in previous works using Eulerian [36] or semi-Lagrangian [34] methods exhibit significantly larger energy conservation errors. These losses in the total energy occur each time a frontal discontinuity forms. The energetic losses incurred in the results of [36] are attributed to the approximation of the advection of the meridional velocity v in the numerical method. The geometric method has two advantages in this regard. First, it is a Lagrangian method so there is no need to approximate an advection term. Second, at each time step the Laguerre tessellation of the fluid domain is defined by the values of θ and v, and so it is adapted to the location of the front, regardless of the chosen resolution. This  (18) and potential energy (19) for the numerical solution of (16) with initial data from Table 1, Row 1, and simulation parameters η = 0.001, h def = 30 and n = 2678. This shows that our numerical scheme is energy-conserving to a high accuracy.
(b) Energy conservation error ε n (defined in (59)) obtained using simulation parameters η = 0.001 and h def = 30. The t coordinate of each orange or blue vertical line is a time when the RMSv curve corresponding to n = 2678 is at a peak or trough, respectively. We see that the energy conservation of the numerical method is worst at peaks and troughs of the RMSv. is reflected by the fact that the order of the energy conservation error was observed to be similar when using different numbers of seeds n.
The qualitative behaviour of our numerical solutions in physical space coincide with those from previous works [8,11,34,36]. To the best of our knowledge, the behaviour of the corresponding seeds in geostrophic space has not been previously illustrated.

Stable normal mode
The results reported in this subsection were obtained using the initial condition and physical parameters from row 2 of Table 1, and the discretisation parameters η = 0.001, h def = 30 and n = 990.
Under the dynamics of the linearised equations (78)-(82), disturbances of the steady shear flow (14) corresponding to the stable perturbation (55) propagate west with constant wave speed where κ and σ are defined by (51) and (52): see Appendix C.4. With our choice of physical parameters this equates to one domain length (i.e. 2L) every 16 days. This supports our numerical results, as can be seen in Figure 4, where we plot the perturbation of the steady potential temperature. Figure 4: Plots of the potential temperature perturbation in the stable parameter regime with initial data from Table 1, Row 2, and simulation parameters η = 0.001, h def = 30 and n = 990. We see that the initial potential temperature perturbation propagates left with roughly constant amplitude. The numerically observed wave speed is in good agreement with the theoretical prediction of the linear instability analysis (Section C.4), which predicts that the wave takes 16 days to cross the domain. In particular, the yellow patch in the bottom right corner at t = 0 crosses half the domain by day 8.
Indeed, Figure 4 shows that the large scale disturbance initially located at the right-most boundary of Ω (see the yellow patch in the bottom right corner of Figure

Numerical evidence for the convergence of the Eady-Boussinesq vertical slice equations to the semi-geostrophic Eady slice equations
We now compare our numerical solutions of the SG Eady slice equations (16) to the numerical solutions of the Eady-Boussinesq vertical slice equations (7) obtained in [29] and [34]. As outlined in [34, Section 2.3], the SG Eady slice equations (8)-(9) can be understood formally as a small Rossby number approximation of (7). This is made concrete by a rescaling argument using a scaling parameter β such that the limit Ro → 0 corresponds to the limit β → 0. In [34], numerical solutions of (7) are obtained for a sequence of decreasing values of β using a semi-Lagrangian method. The resulting RMSv curves are then compared to that from [8, Figure 4.6], which was obtained by solving the SG Eady slice equations (16) numerically using the geometric method. A similar programme is followed in [36] using a compatible finite element method to solve (7) numerically. It is observed that the maximum amplitude of the RMSv obtained in [8] is much greater than that obtained in both subsequent studies, even for small β and high-resolution  (16) was solved with initial data from Table 1, Row 1, and simulation parameters η = 0.01, h def = 30 and n = 2678. The limit β → 0 corresponds to the large scale limit Ro → 0. As β → 0, the blue curves tend towards the solid black curve, at least for small times. This supports the hypothesis that the SG Eady slice equations are the small Rossby number limit of the Eady-Boussinesq vertical slice equations.
simulations. This would not support the hypothesis that the SG Eady slice equations are the limit of the Eady-Boussinesq vertical slice equations as Ro → 0. However, the reason for the discrepancy between the RMSv curves is that the physical parameters used in [8] are not the same as those used in the two subsequent papers, despite being reported as such.
To correct the comparisons made in [34] and [36], we implement the geometric method using the physical parameters listed in those papers. Note that the initial condition used in both [34] and [36] is that of Table 1, Row 3, not the most unstable mode discussed in Section 5.2. When using this initial condition, it is therefore some time before the fastest growing unstable normal mode dominates the numerical solution and the expected growth rate of the RMSv is achieved. To account for this, in [34] and [36] the numerical solutions are translated backwards in time so that the maximum value of v at t = 0 days matches that of [29]. We use the normal mode initial condition from Table 1, Row 1, and do not translate in time. The discretisation parameters we use are η = 0.01, h def = 30 and n = 2678. The resulting RMSv curve is plotted against those from [34] in Figure 5 2 . We see that as β → 0, the RMSv curves of [34] tend towards the RMSv curve of our numerical solution of the SG Eady slice equations (16) (solid black curve in Figure 5), at least for small times. This supports the hypothesis that (8)-(9) are the small Rossby number limit of (7).  Table 2: Relative errors of the numerical solution of (16) with initial data from Table 1 In particular, except at very early times where the SG solution is affected by discretisation errors, the initial growth rate of the RMSv curves from [34] approaches that of the SG numerical solution as β → 0. This supports the initialisation procedure employed in [34] and [36]. For a fair and more detailed comparison, however, it would be desirable to perform further experiments using the semi-Lagrangian method of [34] and the compatible finite-element method of [36] with the normal mode initial condition so that every simulation uses the same initial condition and there is no need to time-translate the solutions.

Sensitivity analysis
In this section we investigate the sensitivity of our numerical solution to the simulation parameters. These are η, the percentage area tolerance supplied to the semi-discrete optimal transport solver; h def , the default time step size in seconds; n, the number of seeds. The initial condition and physical parameters used to generate the results in this section are those from Table 1, Row 1. Figure 6: Plots of the relative errors from Table 2.
Accuracy tests. To evaluate the accuracy of our numerical solutions of (16), we ran three sets of simulations, varying each parameter in turn. For the first set of simulations in which η was varied, we used h def = 30, and n = 1470; in the second set of simulations in which h def was varied, we used η = 0.001, and n = 1470; in the third set of simulations in which n was varied, we used η = 0.01 and h def = 30. The remaining parameter values are listed in the first column of Table 2. For each set of simulations, the numerical solutions were compared to that obtained using the discretisation parameter giving the finest discretisation, which is listed in the final row of the corresponding section of Table 2.
Next we state how we compare solutions with different discretisation parameters. Associated to each solution z of the ODE (32) with corresponding target masses m = (m 1 , . . . , m n ) is a time-dependent family of discrete distributions known as the potential vorticity. A natural way to quantify the discrepancy between two discrete distributions n i=1 m i δ z i and n j=1 m j δ z j is by using the Wasserstein 2-distance (see, for example, [33,Chapter 5]). The Sinkhorn loss introduced in [22, Theorem 1] provides an approximation of the Wasserstein 2-distance that can be computed easily and quickly by solving a regularised optimal transport problem. Indeed, its square-root converges to the Wasserstein 2-distance as the strength of the regularisation goes to zero [19]. If n = n and m i = m i for all i ∈ {1, . . . , n}, one can also consider the weighted Euclidean distance between z and z given by which provides an upper bound on the Wasserstein 2-distance and is much simpler to compute. For each set of simulations, we computed both the square root of the Sinkhorn loss with regularisation parameter ε = 0.001 and the weighted Euclidean distance at several times t and normalised them by the factor where z is the numerical solution of (32) obtained using the simulation parameter giving the finest discretisation, and z i (t) ∞ is the maximum absolute value of the components of z i (t). These values are presented in Table 2 and Figure 6 and, because of the choice of normalisation, can be interpreted as relative errors. The relative error decreases with η, and tends to increase over time, as expected (see Figure  6, top left and bottom left). When quantified using the weighted Euclidean distance, the relative error also decreases with the default time step size h def , but when quantified using the Sinkhorn loss, this is not true at later times (see Figure 6, top middle and bottom middle). The relative error decreases with n at t = 2 days and t = 4 days, but appears to increase slightly for large values of n at t = 6 days and t = 8 days (see Figure 6, top right). Nevertheless, in all cases, the relative errors are small. The unexpected trends may occur because the reference solutions are numerical rather than exact, or because the time step sizes chosen by the adaptive method depend on the simulation parameters.
Behaviour of RMSv. Figure 7a demonstrates that for sufficiently large n the RMSv is well approximated over multiple life-cycles by our numerical results. We deduce from Figures 7b and 7c that the initial decline in the RMSv is due to errors incurred by the spatial discretisation of the initial condition (see Section 4.1), which decrease at a rate faster than n − 1 2 , and do not significantly affect the behaviour of the RMSv after t = 2 days.
Adaptive time step size. Figure 8 demonstrates that the number of times that the time step size is halved per iteration by the adaptive time stepping algorithm correlates with the magnitude of the RMSv, and hence the strength of the frontal discontinuity. This is to be expected because the vertical component ofż i (t) is proportional to the meridional velocity v(c i (z(t)), t) so, roughly speaking, as the RMSv increases, so does the magnitude ofż i (t). As such, the first-order Taylor expansion of w * becomes less accurate as the RMSv increases, so a smaller time step is needed in order for it to generate a good initial guess for the weights for Algorithm 1. (c) RMSv curves up to t = 2 days with large n. The light-blue curve is the RMSv curve of the exact solution of the linearised equations (78)-(82). We see that the initial dip in the RMSv is an artefact of the discretisation. Figure 7: Sensitivity of the RMSv to the number of seeds n. The initial data was taken from Table  1, Row 1, and the discretisation parameters were η = 0.01 and h def = 30.  Table  1, Row 1. The t coordinate of each orange or blue vertical line is a time when the corresponding RMSv curve is at a peak or trough, respectively. The number of time step refinements is maximum at peaks of the RMSv curve, which correspond to peaks in the strength of the frontal discontinuity. Therefore the adaptive time stepping scheme successfully identifies the frontal discontinuities.

Conclusion
In this paper, we recast the geometric method for solving the SG Eady slice equations in the language of semi-discrete optimal transport theory (Section 3), and develop a new implementation using the latest results from semi-discrete optimal transport theory and a novel adaptive time-stepping algorithm that is tailored to the ODE (32) (Algorithm 2). The numerical solutions that we obtain via our implementation support the conjecture that weak solutions of the Eady-Boussinesq vertical slice equations (7) converge to weak solutions of the semi-geostrophic Eady slice equations (8), (9) as Ro → 0 (Section 5.4). Rigorous numerical tests in Section 5.5 demonstrate the sensitivity of the algorithm with respect to the discretisation parameters. To clarify the use of different initial conditions in the literature on the Eady slice problem [8,[34][35][36], we include a linear instability analysis of the steady shear flow (11), validating and extending the work of Eady [14] (see Appendix C). We perform simulations in different physical parameter regimes, which verify the linear instability analysis (Sections 5.2 and 5.3).
The linear instability analysis provides benchmark initial conditions in both stable and unstable parameter regimes. Along with our implementation of the geometric method, this could be used in future work to carry out a more rigorous numerical study of the convergence of weak solutions of the Eady-Boussinesq vertical slice equations to weak solutions of the SG Eady slice equations as Ro → 0, and to explore the behaviour of solutions in different physical parameter regimes.

A Derivation of the semi-discrete transport problem
In this section we derive the semi-discrete optimal transport problem (Definition 3.1) from the stability principle. Consider a geopotential φ : R × [−H/2, H/2] → R that is 2L-periodic in the x 1 direction. The modified geopotential P : R × [−H/2, H/2] → R defined by such that The corresponding total geostrophic energy is The stability principle says that at each point in time this energy is minimised over all periodic rearrangements of particles that preserve potential temperature and absolute momentum. Each such rearrangement corresponds to a collection of sets S i ⊂ R × [−H/2, H/2], i ∈ {1, . . . , n}, satisfying (62) with specified masses Necessarily, The terms (64) are constant over all such collections. This leads to the following minimisation problem: |x − z i − k| 2 dx : (62) and (65) hold .
We show that the semi-discrete optimal transport problem (Definition 3.1) is equivalent to (66). (62) and (65). For each i ∈ {1, . . . , n}, define Note that |C i | = m i for all i ∈ {1, . . . , n}, and {C i } n i=1 is a tessellation of Ω. Then By taking the minimum over {S i } n i=1 , we see that the minimum (66) is greater than or equal to the minimum attained in the semi-discrete optimal transport problem (Definition 3.1).
On the other hand, consider a tessellation {C i } n i=1 of Ω satisfying the mass constraint |C i | = m i for all i ∈ {1, . . . , n}. For each i ∈ {i, . . . , n} define Elementary calculations show that {S i } n i=1 satisfies (62) and (65), and that By taking the minimum over {C i } n i=1 , we see that the minimum (66) is less than or equal to the minimum attained in the semi-discrete transport problem (Definition 3.1). By combining this with the opposite inequality above, we see that the minimum values are equal, and that the stability principle is equivalent to the semi-discrete optimal transport problem, as claimed. Minimisers are related via (67) and (68).
Note that k * (F(x, t), z i (t)) and 1 C i,per (t) (F(x, t)) are piecewise constant in time. Therefore provided that the derivative exists. By (29) and the mass-preserving property of the flow F, On the other hand, integrating the right-hand side of (16) over C i,per (0) gives where c i was defined in (30), and where we used the fact that Here C i (t) := C i (z(t), w * (z(t))) is a non-periodic Laguerre cell (see Definition 3.4). By combining (16), (69) and (70), we derive (31) as desired.

C Linear instability analysis
In this section we study the linear instability of the following steady solution (planar Couette flow) of equations (8) and (9): This was first done by Eady [14], and the unstable perturbations found by Eady were used as initial conditions for the simulations of Williams [35]. We reproduce Eady's results here to make the paper self-contained, to clarify the connection between the initial conditions used here and by [34][35][36], and to derive the stable modes, which are not given in [14]. Note that (71) is not the only steady solution of (8), (9). In fact so is for any φ(z).
For k = 0, ω = 0 is an eigenvalue of (85)-(90) with eigenfunctionsû =v =ŵ = 0,φ arbitrary, and θ determined by (89). This eigenvalue corresponds to the family of steady planar Couette flows (72). From now on we assume that k = 0. We now eliminateû,v,φ andθ. This can be done by noting thatû Taking the z-derivative of (88), multiplying (89) by ikπg/(θ 0 f L), and adding the resulting equations gives Define the change of variables It follows from (93) that W satisfies the ODE The general solution of this ODE is given by This can be used to write the general solution of (93). First we use the boundary conditions (90) to find the eigenvalues ω.
Introduce the dimensionless growth rate Thenŵ(z) = W (Z(z)), wherẽ The boundary conditions (90) mean that Therefore, from (95), we read off that where M is the 2-by-2 matrix with components Non-trivial solutions are then those for which the vector (A, B) belongs to the kernel of M . Imposing the condition that det(M ) = 0 yields 1 + κ 2 − ω 2 sinh 2κ = 2κ cosh 2κ.
Since κ ∈ R, any ω that satisfies the dispersion relation (97) is either real or purely imaginary. Correspondingly the growth rate ω is either purely imaginary or real.
It is worth observing that there are no exponentially growing modes when where κ crit = 1.19968 (to 6 s.f.) is the smallest positive root of (99). We study this parameter regime in the following section. It can be shown that 2κ coth 2κ − 1 − κ 2 is maximised when κ = ±κ * with κ * = 0.803058 (to 6 s.f.). This is achievable by an integer mode number k if the Burger number satisfies for some k ∈ Z. Take k = 1 and the values of N , L, f , s, θ 0 , g given in Section 5.1. We choose H so that (101) is satisfied. This gives H = 10224.85 m.
The fastest growing mode can be rewritten as ω = u z Bu H L σ(κ * ).
In particular, we can read off that it is proportional to u z and the aspect ratio H/L of the box, and inversely proportional to Bu. Next we compute the eigenfunctions for the unstable modes. By (96) we can writê w(z) = W (z(z) + iσ),z = 2κ H z.

C.4 Stable modes
In this section we study the parameter regime Bu > Bu crit , where the eigenvalues ω are purely imaginary for all k ∈ Z. This means that the steady Couette solution (71) is linearly stable with respect to normal mode perturbations. It turns out, however, that the normal modes do not form a complete basis, in the following sense. For each mode number k, there are only two solutions ±σ(κ k ) of (99), and therefore only two eigenvalues ±ω k . Consequently the boundary value problem (93), (90) only has two eigenvalues, and not a countable set of eigenvalues as one might expect. This means that it is not possible to represent every initial condition of the linear perturbation equations (78)-(82) as a sum of normal modes (eigenfunctions). Therefore stability with respect to normal mode perturbations does not guarantee stability with respect to all perturbations. The origin of this problem is the denominator ω − ikπ L u in equation (93), which can vanish. For a closely related but simpler baroclinic instability problem, Pedlosky [30] showed that any initial perturbation can be represented by supplementing the discrete spectrum of normal modes with a continuous spectrum. By doing this he showed that the linear stability of the steady solution is in fact determined by the incomplete set of normal modes. Our simulations suggest that the same is true for our problem. We plan to explore this in a future paper.
It follows that Substituting these expressions into (103) and taking k = 1 gives the initial condition that we use in Section 5.3. In particular, v s (x, z) = λv 1 (x, z, 0) and θ s (x, z) = λθ 1 (x, z, 0), where λ was defined in (104). Observe that this perturbation is a travelling wave with wave speed

D Derivatives of K
Recall that K was defined in equation (26). The gradient and Hessian of K are given by ∇K = m − m, where m is the cell-area map from Definition 3.3. For a proof, see for example [23,28]. The system (40) defining the k th Newton direction d (k) can therefore be written as −D w m z, w (k) d (k) = m z, w (k) − m, d (k) · e n = 0.
For i = j, ∂m i ∂w j = − 1 2 k∈K length(e(i, j, k)) where e(i, j, k) is the edge between the non-periodic Laguerre cells with seeds z i + k and z j , which is defined by e(i, j, k) = x ∈ R × − H 2 , H 2 : ∀ m ∈ {1, . . . , n}, l ∈ K, Note that this set may be empty, in which case ∂m i /∂w j = 0. The diagonal entries of D w m are The expressions for the derivatives of m are proved for example in [23,28] for non-periodic Laguerre tessellations and in [5] for periodic Laguerre tessellations.