Rotating shallow water flow under location uncertainty with a structure-preserving discretization
RRotating shallow water flow under location uncertaintywith a structure-preserving discretization
R¨udiger Brecht † , Long Li § , Werner Bauer ‡ , Etienne M´emin §† Department of Mathematics and Statistics, Memorial University of Newfoundland,St. John’s (NL) A1C 5S7, Canada § Inria/IRMAR, Campus universitaire de Beaulieu, Rennes, France ‡ Imperial College London, Department of Mathematics, 180 Queen’s Gate, London SW7 2AZ,United Kingdom.
Abstract
We introduce a new representation of the rotating shallow water equations based on a stochastictransport principle. The derivation relies on a decomposition of the fluid flow into a large-scalecomponent and a noise term that models the unresolved small-scale flow. The total energy ofsuch a random model is demonstrated to be preserved along time for any realization. To preservethis structure, we combine an energy (in space) preserving discretization of the underlyingdeterministic model with approximations of the stochastic terms that are based on standard finitevolume/difference operators. This way, our method can directly be used in existing dynamicalcores of global numerical weather prediction and climate models. For an inviscid test case on thef-plane we use a homogenous noise and illustrate that the spatial part of the stochastic schemepreserves the total energy of the system. Moreover, using an inhomogenous noise, we show for abarotropically unstable jet on the sphere that the proposed random model better captures thestructure of a large-scale flow than a comparable deterministic model.
Plain summary
The motion of geophysical fluids on the globe needs to be modelled to get some insights oftomorrow’s climate. These forecasts must be precise enough while remaining computationallyaffordable. An ideal system should also deliver, across time, an accurate measurement of the un-certainties introduced through physical or numerical approximations. To address these issues, weuse the rotating shallow water equations, which provide a simplified version of the dynamics, anda stochastic representation of the unresolved small-scale processes. The former is approximatedwith a structure preserving numerical model enabling the conservation of physical quantitiessuch as mass and energy and the latter is modelled by the location uncertainty framework thatrelies on stochastic transport and has the great advantage to be energy conserving. Our methodcan directly be used in existing dynamical cores of global numerical weather prediction andclimate models. Numerical results illustrate the energy conservation of the numerical model.Simulating a barotropically unstable jet on the sphere, we demonstrate that the random modelbetter captures the structure of a large-scale flow than a comparable deterministic model. Therandom dynamical system is also shown to be associated with good uncertainty representation.1 a r X i v : . [ phy s i c s . f l u - dyn ] F e b Introduction
Numerical simulations of the Earth’s atmosphere and ocean plays an important role in developingour understanding of weather forecasting. A major focus lies in determining the large scale flowcorrectly, which is strongly related to the parameterizations of sub-grid processes (Frederiksenet al., 2013). The non-linear and non-local nature of the dynamical system make the large-scaleflow structures interact with the smaller components. The computational expense for solvingthe Kolmogorov scales (Pope, 2000) of a geophysical flows is fare beyond reach today and likelyin the future. Thus, the effect of unresolved scales has to be modeled or parametrized.For several years, there is a growing interest in geophysical sciences to incorporate a stochas-tic representation (Franzke and Majda, 2006; Majda et al., 2008; Grooms and Majda, 2014;Gottwald et al., 2017) of the small-scale processes. In this study, we propose to stick to a spe-cific stochastic model, the so-called
Location Uncertainty (LU) derived by M´emin (2014), whichemerges from a decomposition of the Lagrangian velocity into a time-smooth drift and a highlyoscillating uncertainty term. Such random model allows us to develop by stochastic calculus anew stochastic transport operator (Resseguier et al., 2017a) for the extensive scalars. In partic-ular, this transport operator involves a multiplicative random forcing, a heterogeneous diffusionand a corrected advection resulting from the inhomogeneity of the random flow. This stochastictransport principle has been used as a fundamental tool to derive stochastic representations oflarge-scale geophysical dynamics (Resseguier et al., 2017a; Chapron et al., 2018; Bauer et al.,2020a). In the present work, we use this mathematical principle together with some physicalconservation laws to derive a stochastic version of the rotating shallow water (RSW) system.One strong property of this random model is that it preserves the total energy of the resolvedflow in time for each realization.Recently, the LU model performed very well in Resseguier et al. (2017b,c); Bauer et al.(2020a,b) when studying oceanic quasi-geostrophic flows. It was found to be more accurate inpredicting the extreme events, in diagnosing the frontogenesis and filamentogenesis, in struc-turing the large-scale flow and in reproducing the long-term statistics. Besides, Chapron et al.(2018) investigated the Lorentz-63 test case and demonstrated that the LU model was moreeffective in exploring the range of the strange attractors compared to classical models.In this work, the performance of the LU model is assessed for the numerical simulation of theRSW system, which can be considered as the first step towards developing numerical randomglobal climate models. This is the first time that the LU model is implemented for the dynamicsevolving on the sphere.We propose to combine the discrete variational integrator for RSW fluids as introduced in(Bauer and Gay-Balmaz, 2019) and (Brecht et al., 2019) with the numerical LU model in or-der to mimic the continuous conservation properties. Variational integrators are designed byfirst discretizing the given Lagrangian, and then by deriving a discrete system of associatedEuler-Lagrange equations from the discretized Lagrangian (see Marsden and West (2001)). Theadvantage of this approach is that the resulting discrete system inherits several important prop-erties of the underlying continuous system, notably a discrete version of Noether’s theorem thatguarantees the preservation of conserved quantities associated to the symmetries of the discreteLagrangian (see Hairer et al. (2006)). Variational integrators also exhibit superior long-termstability properties. Therefore, they typically outperform traditional integrators if one is inter-ested in long-time integration or the statistical properties of a given dynamical system. Thebenefit of the proposed method that relies on a combination of a variational integrator witha potentially differently approximated LU model is that it can directly be applied to existingdynamical cores of numerical weather prediction and climate models.Apart from taking into account the unresolved processes, it is paramount in uncertaintyquantification and ensemble forecasting to model the uncertainties along time (Resseguier et al.,2020). For a long time, operational weather forecast centres had relied on random perturbations2f initial conditions (PIC) to spread the ensemble forecasts. However, in the application ofdata assimilation to geophysical fluid dynamics, such PIC model is known to under-estimatethe true uncertainty compared to the observations (Gottwald and Harlim, 2013; Franzke et al.,2015). Hence, an assimilation system is overconfident for such a random model. To overcomethis issue, the covariance inflation method (Anderson and Anderson, 1999) is often adopted,in which the ensemble covariance is increased by a carefully tuned parameter. In the presentwork, we compare the reliability of the ensemble spread of such a PIC model with our RSW-LUsystem, under the same strength of uncertainty.The remainder of this paper is structured as follows. Section 2 describes the basic principles ofthe LU model and the derivation of the rotating shallow water system under LU associated withthe energy conservation property. Section 3 explains the parameterizations of the uncertaintyand the numerical discretization of the stochastic dynamical system. Section 4 discusses thenumerical results for an inviscid test case with homogeneous noise and a viscous test case withheterogeneous noise. Finally, in Section 5 we draw some conclusions and provide an outlook forfuture work.
In this section, we first review the LU model introduced by M´emin (2014), then we derivethe rotating shallow water equations under LU, denoted as RSW–LU, following the classicalstrategy as shown in Vallis (2017). In particular, we demonstrate one important characteristicof the RSW–LU, namely it preserves the total energy of the large-scale flow.
The LU model is based on a temporal-scale-separation assumption of the following stochasticflow: d X t = w ( X t , t ) d t + σ ( X t , t ) d B t , (2.1)where X is the Lagrangian displacement defined within the bounded domain Ω ⊂ R d ( d =2 or 3), w is the large-scale velocity that is both spatially and temporally correlated, and σ d B t is the small-scale uncertainty (also called noise) term that is only correlated in space. The spatialstructure of such noise is specified through a deterministic integral operator σ : ( L ( Ω )) d → ( L ( Ω )) d , acting on square integrable vector-valued functions f ∈ ( L ( Ω )) d , with a boundedkernel ˘ σ such that σ [ f ]( x , t ) = (cid:90) Ω ˘ σ ( x , y , t ) f ( y ) d y , ∀ f ∈ ( L ( Ω )) d . (2.2)The randomness of such noise is driven by the cylindrical I d -Wiener process B t (Da Prato andZabczyk, 2014). The fact that the kernel is bounded, i.e. sup ( x , y ) ∈ Ω | ˘ σ ( x , y ) | < + ∞ , impliesthat the operator σ is Hilbert-Schmidt on ( L ( Ω )) d . Therefore, the resulting small-scale flow σ d B t is a centered (of null ensemble mean) Gaussian process with the following covariancetensor , denoted as Q , being well-defined: Q ( x , y , t, s ) = E (cid:104)(cid:0) σ ( x , t ) d B t (cid:1)(cid:0) σ ( y , s ) d B s (cid:1) T (cid:105) = δ ( t − s ) d t (cid:90) Ω ˘ σ ( x , z , t ) ˘ σ T ( y , z , s ) d z , (2.3)where E stands for the expectation and δ is the Kronecker symbol. The strength of the noiseis measured by its variance , denoted as a , which is given by the diagonal components of thecovariance per unit of time: a ( x , t ) (cid:52) = Q ( x , x , t, t ) / d t = σσ T ( x , t ) . (2.4)3e remark that such variance tensor a has the same unit as a diffusion tensor (m · s − ) andthat the density of the turbulent kinetic energy (TKE) can be specified by tr( a ) / d t .The previous representation (2.2) is a general way to define the noise in LU models. Inparticular, the fact that σ is Hilbert-Schmidt ensures that the covariance operator per unit oftime, Q / d t , admits an orthogonal eigenfunction basis { Φ n ( • , t ) } n ∈ N weighted by the eigenvaluesΛ n ≥ (cid:80) n ∈ N Λ n < ∞ . Therefore, one may equivalently define the noise and itsvariance, based on the following spectral decomposition: σ ( x , t ) d B t = (cid:88) n ∈ N Φ n ( x , t ) d β nt , (2.5a) a ( x , t ) = (cid:88) n ∈ N Φ n ( x , t ) Φ T n ( x , t ) , (2.5b)where β n denotes n independent and identically distributed (i.i.d.) one-dimensional standardBrownian motions.The core of LU models is based on a stochastic Reynolds transport theorem (SRTT), intro-duced by M´emin (2014), which describes the rate of change of a random scalar q transported bythe stochastic flow (2.1) within a flow volume V . In particular, for incompressible small-scaleflows, ∇· σ = 0, the SRTT can be written asd t (cid:16) (cid:90) V ( t ) q ( x , t ) d x (cid:17) = (cid:90) V ( t ) (cid:0) D t q + q ∇· ( w − w s ) (cid:1) d x , (2.6a) D t q (cid:52) = d t q + ( w − w s ) ·∇ q d t + σ d B t ·∇ q − ∇· ( a ∇ q ) d t, (2.6b)d t q (cid:52) = q t +d t − q t , w s (cid:52) = 12 ∇· a , (2.6c)in which the stochastic transport operator D t (Resseguier et al., 2017a) and the Itˆo-Stokes drift(ISD) w s (Bauer et al., 2020a) are included. The latter term arises from the effect of statisticalinhomogeneity of the small-scale flow on the large-scale component, which can be consideredas a generalization of the Stokes drift in ocean circulations. In the definition of the stochastictransport operator in (2.6b), the first term on the right-hand side (RHS), defined in (2.6c),stands for a forward time-increment of q at a fixed point x , and the last two terms describe,respectively, a backscattering from the small-scales to the large-scales and an inhomogeneousdiffusion at the small-scales. In particular, for an isochoric flow with ∇· ( w − w s ) = 0, one mayimmediately deduce from (2.6a) the following transport equation of an extensive scalar: D t q = 0 , (2.7)where the energy of such random scalar q is globally conserved, as shown in Resseguier et al.(2017a):d t (cid:16) (cid:90) Ω q d x (cid:17) = (cid:16) (cid:90) Ω q ∇· ( a ∇ q ) d x (cid:124) (cid:123)(cid:122) (cid:125) Energy loss by diffusion + 12 (cid:90) Ω ( ∇ q ) T a ∇ q d x (cid:124) (cid:123)(cid:122) (cid:125) Energy intake by noise (cid:17) d t = 0 . (2.8)Indeed, this can be interpreted as a process where the energy brought by the noise is exactlycounter-balanced by that dissipated from the diffusion term.4 .2 Derivation of RSW–LU This section describes in detail the derivation of the RSW–LU system. We remark that aformulation of the shallow water equations under LU in a non-rotating frame is outlined byM´emin (2014), whereas the new model that we present in this work is fully stochastic andincludes rotation such that it is suited for simulations of geophysical flows.The above SRTT (2.6a) and Newton’s second principle allow us to derive the following (three-dimensional) stochastic equations of motions in a rotating frame (Resseguier et al., 2017a; Baueret al., 2020a):
Horizontal momentum equation : D t u + f × (cid:0) u d t + σ H d B t (cid:1) = − ρ ∇ H (cid:0) p d t + d p σt (cid:1) + ν ∇ (cid:0) u d t + σ H d B t (cid:1) , (2.9a) Vertical momentum equation : D t w = − ρ ∂ z (cid:0) p d t + d p σt (cid:1) − g d t + ν ∇ (cid:0) w d t + σ z d B t (cid:1) , (2.9b) Mass equation : D t ρ = 0 , (2.9c) Continuity equation : ∇· (cid:0) w − w s (cid:1) = 0 , ∇· σ = 0 , (2.9d)where u (resp. σ H d B t ) and w (resp. σ z d B t ) are the horizontal and vertical components ofthe three-dimensional large-scale flow w (resp. the small-scale flow σ d B t ); f = (2 ˜ Ω sin Θ) k isthe Coriolis parameter varying in latitude Θ, with the Earth’s angular rotation rate ˜ Ω and thevertical unit vector k = [0 , , T ; ρ is the fluid density; ∇ H = [ ∂ x , ∂ y ] T denotes the horizontalgradient; p and ˙ p σt (cid:52) = d p σt / d t (informal definition) are the time-smooth and time-uncorrelatedcomponents of the pressure field, respectively; g is the Earth’s gravity value and ν is the kine-matic viscosity. For the following derivation of the shallow water equations we drop the viscousterms.In order to model the large-scale circulations in the atmosphere and ocean, the hydrostaticbalance approximation is widely adopted (Vallis, 2017). Under a small aspect ratio, H / L (cid:28) L and H the horizontal and vertical scales of the motion, the acceleration term D t w on theleft-hand side (LHS) of Equation (2.9b) has a lower order of magnitude than the RHS terms,hence the vertical momentum equation reduces to ∂ z (cid:0) p d t + d p σt (cid:1) = − g d t. (2.10a)According to the Doob’s theorem – unique decomposition of a semimartingale process (Kunita,1997), the previous equation is equivalent to ∂ z p = − ρg, ∂ z d p σt = 0 . (2.10b)Integrating vertically these hydrostatic balances (2.10b) from 0 to z (see Figure 1) under aconstant density ρ , we have p ( x, y, z, t ) = p ( x, y, t ) − ρ gz, (2.10c)d p σt ( x, y, z, t ) = d p σt ( x, y, , t ) , (2.10d)where p denotes the pressure at the bottom of the basin ( z = 0). Following Vallis (2017), weassume that the weight of the overlying fluid is negligible, i.e. p ( x, y, η, t ) ≈ η the heightof the free surface, leading to p = ρ gη . This allows us to rewrite Equation (2.10c) such thatfor any z ∈ [0 , η ] we have p ( x, y, z, t ) = ρ g (cid:0) η ( x, y, t ) − z (cid:1) . (2.10e)5ubsequently, the pressure gradient forces in the horizontal momentum equation (2.9a) reduceto − ρ ∇ H (cid:0) p d t + d p σt (cid:1) = − g ∇ H η − ρ ∇ H d p σt , (2.10f)which do not depend on z according to Equations (2.10e) and (2.10d). Therefore, the accelerationterms on the LHS of Equation (2.9a) must not depend on z , hence the shallow water momentumequation can be written as D H t u + f × (cid:0) u d t + σ H d B t (cid:1) = − g ∇ H η d t − ρ ∇ H d p σt , (2.11a) D H t u (cid:52) = d t u + (cid:0) ( u − u s ) d t + σ H d B t (cid:1) · ∇ H u − ∇ H · (cid:0) a H ∇ H u (cid:1) d t, (2.11b) u s (cid:52) = 12 ∇ H · a H , a = (cid:18) a H a Hz a Hz a z (cid:19) , (2.11c)where D H t is the horizontal stochastic transport operator, u s is the two-dimensional ISD, a H , a z and a Hz are the horizontal, vertical and cross components of the three-dimensional variancetensor a . Note that Equation (2.11a) is valid only when the cross component a Hz is verticallyindependent, i.e. ∂ z a Hz = 0. For instance, one may consider that the horizontal small-scale flow σ H d B t is spatially uncorrelated with the vertical small-scale flow σ z d B t , i.e. a Hz = 0.In order to derive the shallow water mass equation, let us first integrate vertically the con-tinuity equation (2.9d) from the bottom topography η b to the free surface η (see Figure 1):( w − w s ) | z = η − ( w − w s ) | z = η b = − h ∇ H · ( u − u s ) , (2.12a) σ d B t | z = η − σ d B t | z = η b = − h ∇ H · σ H d B t , (2.12b)where h = η − η b denotes the thickness of the water column. On the other hand, a small vertical(Eulerian) displacement at the top and the bottom of the fluid leads to a variation of the positionof a particular fluid element (Vallis, 2017): (cid:0) ( w − w s ) d t + σ d B t (cid:1)(cid:12)(cid:12) z = η = D H t η, (2.12c) (cid:0) ( w − w s ) d t + σ d B t (cid:1)(cid:12)(cid:12) z = η b = D H t η b . (2.12d)Combining Equations (2.12), we deduce the following stochastic mass equation: D H t h + h ∇ H · (cid:0) ( u − u s ) d t + σ H d B t (cid:1) = 0 . (2.13)The above two equations (2.11a) and (2.13) constitute a general formulation of the RSW–LUsystem. Using different levels of noise strength in such a stochastic system allows us to modeldifferent physical regimes of the large-scale flow. To characterise these regimes, Resseguier et al.(2017b) introduced the following scaling number (cid:15) = T σ T T KEM KE , (2.14a)6 igure 1. Illustration of a single-layered shallow water system (inspired by Vallis (2017)). h is the thickness of awater column, η is the height of the free surface and η b is the height of the bottom topography. As a result, wehave h = η − η b . where T and T σ are the correlation time scales of the large-scale flow and the small-scale com-ponent, respectively. The mean kinetic energy scale ( M KE ) is given by U with U = L / T thetypical velocity scale, and the turbulent kinetic energy scale ( T KE ) is defined by A / T σ with A the magnitude of the variance tensor a . As such, the dimensional noise associated with itsdimensional variance can be specified by σ H d B t = √ (cid:15) L ( σ H d B t ) (cid:48) , a = (cid:15) U L a (cid:48) , (2.14b)where • (cid:48) denotes adimensional variables. From expressions (2.14), one may easily concludethat the greater the scaling number (cid:15) , the stronger the noise σ H d B t (with higher variance a ).Furthermore, as interpreted in Resseguier et al. (2017c), a strong noise ( (cid:15) (cid:29)
1) modifies theclassical geostrophic equilibrium of the large-scale flow by including some correction terms to theisobaric velocities. In the present work, only moderate noise ( (cid:15) ∼
1) is adopted for the RSW–LUsystem. Under such assumption, the small-scale flow becomes approximately geostrophic andincompressible, i.e. f × σ H d B t ≈ − ρ ∇ H d p σt and ∇ H · σ H d B t = 0. As a result, the RSW–LUsystem simplifies to D H t u + f × u d t = − g ∇ H η d t, (2.15a) D H t h + h ∇ H · ( u − u s ) d t = 0 , (2.15b) ∇ H · σ H d B t = 0 . (2.15c)We remark that an additional incompressible constraint must be imposed on the horizontal ISD, i.e. ∇· u s = 0, so that the previous system preserves the total energy of the large-scale flow.This will be shown in the subsequent section. For the sake of readability, in the following wedrop the symbol H for all horizontal variables. This section demonstrates the energy conservation of the RSW–LU system (2.15). Let us recallthat the density of the kinetic energy (KE) and of the potential energy (PE) of the large-scaleflow in the shallow water system (Vallis, 2017) is, respectively, given byKE = (cid:90) h ρ | u | d z = ρ h | u | , (2.16a)7E = (cid:90) h ρ gz d z = ρ gh . (2.16b)The density of total energy is defined as the sum of them:E = KE + PE (2.16c)where | u | = u · u and we assume that ρ = 1 and the bottom is flat, i.e. η b = 0 for algebraicsimplicity.In order to explain the conservation of energy more concisely, we adopt an equivalent Stratonovichrepresentation of the RSW–LU system (2.15), namely D t ◦ u + f × u d t = − g ∇ h d t, (2.17a) D t ◦ h + h ∇· ( u − u s ) d t = 0 , (2.17b) f × σ ◦ d B t = − ∇ d t ◦ p σ , ∇· σ ◦ d B t = 0 , (2.17c) D t ◦ u (cid:52) = d t ◦ u + (cid:0) ( u − u s ) d t + σ ◦ d B t (cid:1) ·∇ u, (2.17d)where d t ◦ (cid:52) = u t + dt/ − u t − dt/ stands for a central time-increment based on the definition ofStratonovich integrals and D t ◦ denotes the stochastic transport operator under Stratonovichnotations. We remark that the equivalence between the Itˆo form (2.11b) and the Stratonovichform (2.17d) are fully detailed in Appendix C of Bauer et al. (2020a). As shown by Kunita(1997), Stratonovich integrals are defined such that the chain rule and the integration-by-partformula of ordinary calculus holds. In particular, for two random tracers f and g , we haved t ◦ ( f g ) = f d t ◦ g + g d t ◦ f. (2.18a)Therefore, from the definition of the Stratonovich transport operator (2.17d), we deduce thefollowing product rule: D t ◦ ( f g ) = g D t ◦ f + f D t ◦ g. (2.18b)Applying this rule on the definition of PE (2.16b) together with the mass equation (2.15b), D t ◦ PE = gh D t ◦ h = − gh ∇· ( u − u s ) d t, (2.19a)or D t ◦ PE + 2PE ∇· ( u − u s ) d t = 0 . (2.19b)Similarly, from both mass equation and momentum equation in (2.15), we derive the evolutionof KE (2.16a): D t ◦ KE = h u · D t u + 12 | u | D t ◦ h = − u ·∇ (cid:0) gh (cid:1) d t − h | u | ∇· ( u − u s ) d t, (2.19c)8oting that u · ( f × u d t ) = 0 and recalling that η b = 0, which yields D t ◦ KE + u ·∇ PE d t + KE ∇· ( u − u s ) d t = 0 . (2.19d)Subsequently, we deduce the evolution of the density of total energy: D t ◦ E + ∇· ( u PE) d t − PE ∇· u s d t + E ∇· ( u − u s ) d t = 0 . (2.20a)Expanding the Stratonovich transport operator (2.17d), the previous equation can be re-writtenas d t ◦ E + ∇· (cid:0) F d t + F ◦ d B t (cid:1) = PE ∇· u s d t, (2.20b)where F (cid:52) = ( u − u s ) E + u PE and F ◦ d B t (cid:52) = E σ ◦ d B t are the total energy flux due to thecorrected large-scale drift u − u s and the noise component, respectively. The additional termPE ∇· u s stands for sources or sinks of the potential energy due to the compressibility of theISD. In particular, if we assume that the ISD is incompressible, i.e. ∇· u s = 0, the evolution ofthe energy density reduces tod t ◦ E + ∇· (cid:0) F d t + F ◦ d B t (cid:1) = 0 . (2.21a)If the fluid domain has zero boundary conditions (e.g. the normal velocities vanish on each wallor there are no boundaries at all as on the sphere), then one can show that the total energy,E (cid:52) = (cid:82) Ω E( x , t )d x , is invariant in time:d t ◦ E = (cid:90) Ω d t ◦ E d x = − (cid:90) ∂Ω ( F d t + F ◦ d B t ) · n d l = 0 , (2.21b)where ∂Ω and n denote the domain’s boundaries and the unit normal vector, respectively.In sum, in this work we propose the following RSW–LU system that preserves the globalenergy of the large-scale flow in time for any realization of a random noise: Conservation of momentum : D t u + f × u d t = − g ∇ η d t, (2.22a) Conservation of mass : D t h + h ∇· u d t = 0 , (2.22b) Random geostrophic constraint : f × σ d B t = − ρ ∇ d p σt , (2.22c) Incompressible constraints : ∇· σ d B t = 0 , ∇· u s = 0 , (2.22d) Conservation of energy :d t (cid:90) Ω ρ (cid:0) h | u | + gh (cid:1) d x = 0 . (2.22e)Note that for a sufficiently small noise ( σ ≈ D t / d t , reducesto the material derivative. In order to perform a numerical simulation of the RSW–LU (2.22), the noise term σ d B t andthe variance tensor a have to be a priori parametrized. Then an adequate discretization inspace-time have to be specified for solving the dynamical system. This section describes thesetwo aspects. 9 .1 Parameterizations of noise In the following, we present two different kinds of spatial structure for the noise – homogeneousand heterogeneous. The first one is easy-to-implement, in particular when considering noise thatrespects the incompressible constraints (2.22d). We use such homogeneous noise to study thenumerical energy behaviour, as shown in Section 4.1. On the other hand, because heterogeneousnoise has more physical meaning, we will use the latter when studying realistic complex flows.As shown in Bauer et al. (2020a), heterogeneous noise induces a structuration of the large-scaleflow through the inhomogeneity of the small-scale flow. In Section 4.2, such heterogeneous noiseis adopted for identifying the barotropic instability of a mid-latitude jet.
From the general definitions (2.2) and (2.4), a homogeneous noise means that its correlationoperator σ is a convolution operator and the variance tensor a reduces to a constant matrix(independent of any position in the fluid domain). Furthermore, to ensure that a two-dimensionalnoise is incompressible, Resseguier et al. (2017b) proposed an isotropic model defined througha random stream function σ ( x ) d B t = ∇ ⊥ (cid:0) ˘ ϕ (cid:63) d B t (cid:1) ( x ) , (3.1)where ∇ ⊥ = [ − ∂ y , ∂ x ] T denotes the perpendicular gradient and ˘ ϕ (cid:63) d B t stands for the randomstream function with a convolution kernel ˘ ϕ (and the symbol (cid:63) denotes a convolution). Asshown in Resseguier et al. (2017b, 2020), both isotropy and incompressibility of the noise (3.1)result in a (constant) diagonal variance tensor a I with the eddy-viscosity-like coefficient a and the two-dimensional identity matrix I . In fact, the divergence-free constraint of the ISD inEquation (2.22d) is naturally satisfied (since ∇· u s = ∇· ∇· ( a I ) = 0). As discussed at the endof Section 2.2, for the RSW–LU system (2.22) under geostrophic noise , f × σ H d B t ≈ − ∇ H d p σt ,one can identify, for a constant Coriolis parameter f , the random pressure d p σt with the proposedrandom stream function by d p σt = f ˘ ϕ (cid:63) d B t .In practice, the convolution kernel ˘ ϕ is specified by three parameters: a fixed omni-directionalspectrum slope s , a band-pass filter f BP with support in the range of two wavenumbers κ m and κ M , and an eddy-viscosity-like coefficient a . In fact, the Fourier transform of the random streamfunction ˘ ϕ (cid:63) d B t can be defined as: (cid:92) ˘ ϕ (cid:63) d B t ( k ) (cid:52) = A √ ∆ t f BP ( (cid:107) k (cid:107) ) (cid:107) k (cid:107) − α (cid:98) ξ t ( k ) with α = (3 + s ) / , (3.2)where (cid:98) • denotes the Fourier transform coefficient, ξ t is a space-time white noise, and A is aconstant to ensure E (cid:13)(cid:13) σ d B t (cid:13)(cid:13) = 2 a ∆ t (see Equations (2.3) and (2.4)) with ∆ t the size of onetime stepping. In the simulations, the maximal wavenumber k M of the noise can usually bechosen as the effective resolution cutoff, the minimal wavenumber can be set to k m = k M /
2, andthe theoretical spectrum slope of a two-dimensional flow is given by s = −
3. The noise strengthparameter a will be specified in Section 4.1. The homogeneous noise defined in Section 3.1.1 is quite simple to construct and to interpret,however, it lacks to represent physically important contributions of the small-scale to the largescale flow, which is crucial in order to accurately model realistic scenarios in geophysical fluiddynamics. For this reason, two parameterizations of the heterogeneous noise are presented inthe following.These approaches result from the spectral decomposition (2.5) used to construct the eigen-function basis of the spatial covariance. In practice, we work with a finite set of Empirical10rthogonal Functions (EOFs) of the small-scale Eulerian velocity rather than with the La-grangian displacement. The first method for estimating the EOFs is an off-line procedure basedon the Proper Orthogonal Decomposition (POD) technique of high-dimensional data in whichthe EOFs are assumed to be time-independent, whereas the second one is an on-line estimationfrom a coarse-grid simulation where the EOFs are time-dependent. As will be shown in Section4.2, the former allows for incorporating data into the dynamical model and is more suitable formid-term simulations, yet the latter is independent from observations and is more adequate forlong-term simulations.
Off-line learning of EOFs
Let us consider a set of velocity snapshots { u o ( x , t i ) } i =1 ,...,N t ,that have been a priori coarse-grained from high-dimensional data using a low-pass filter (suchas the sharp spectral filter of Pope (2000) often used in large eddy simulations). Applying thesnapshot POD procedure (Sirovich, 1987) for the fluctuations u (cid:48) o = u o − u o (where • denotesa temporal average) enables us to build a set of EOFs { φ i } i =1 ,...,N t . In addition, we supposethat the fluctuations of the large-scale flow live in a subspace spanned by { φ i } i =1 ,...,m − (with m < N t ) and that the small-scale random drift σ d B t / ∆ t lives in the complemented subspacespanned by { φ i } i = m,...,N t such that1∆ t σ ( x ) d B t = N t (cid:88) i = m (cid:112) λ i φ i ( x ) ξ i , t a ( x ) = N t (cid:88) i = m λ i φ i ( x ) φ T i ( x ) , (3.3)where λ i is the eigenvalue associated to the spatial mode φ i and ξ i is a standard Gaussianvariable. In practice, there exists an opening question in (3.3), that is how to adequately choosethe “splitting mode” φ m . Recently, Bauer et al. (2020b) proposed to fix it by comparing thetime-averaged energy spectrum of the observations and the one from a coarse-grid deterministicsimulation. On-line learning of EOFs
The previously described data-driven estimation of EOFs is aquite efficient procedure. However, such observation data, either from direct measurements orfrom high-dimensional simulations, are not always available. Therefore, Bauer et al. (2020a);Resseguier et al. (2020) proposed an alternative approach in which some local fluctuations,called pseudo-observations (PSO), are generated directly from a coarse-grid simulation. Then,the singular value decomposition (SVD) is applied on those PSO to estimate a set of EOFssuch that the noise associated with its variance tensor will be built in the same way as in (3.3).Finally, the magnitude of the noise and variance should be scaled down to smaller scales basedon a similarity analysis (Kadri Harouna and M´emin, 2017).In the following, we describe in more details both the generation of PSO and the scalingtechnique. The approach proposed here defines N o PSO (denoted as u (cid:48) ) at each grid point. Fora given time t and a current coarse velocity u , we build the PSO by sliding a local windowof size N w × N w over the spatial grid (with N w the grid number in one direction of the localwindow). We denote the spatial scale of the window by L = N w l , where l is the smallest scale ofthe simulation. At every grid point x i,j , we list the N w velocity values contained in the windowcentered at that point: I ( x i,j , t ) (cid:52) = (cid:26) u ( x p,q , t ) (cid:12)(cid:12)(cid:12)(cid:12) | p − i | ≤ N w − , | q − j | ≤ N w − (cid:27) . (3.4)Note that appropriate boundary conditions (replication, periodicity, etc.) are adopted whenlooking at a point on the border. Then, independently for each n ∈ { , . . . , N o } and for eachpoint x i,j , we set the value of the PSO u (cid:48) ( x i,j , t, n ) by randomly choosing a value in the set I ( x i,j , t ). After this, we average over the realization index n . Then, from the SVD we obtain a11et of EOFs { φ ( L ) i } i =1 ,...,N o , and a spectral representation of the small-scale velocity:1∆ t σ ( L ) ( x , t ) d B t = N o (cid:88) i =1 φ ( L ) i ( x , t ) ξ i . (3.5a)Since the PSO u (cid:48) have been generated at a spatial scale of the window L = N w l , they must bescaled down to the “simulation scale” l . As such, the variance tensor a of the small-scale flow isrescaled according to a turbulence-power-law coefficient(Kadri Harouna and M´emin, 2017) suchthat a ( l ) = (cid:18) lL (cid:19) / a ( L ) , (3.5b)where a ( L ) and a ( l ) are the variance tensors at the scales L and l respectively. Finally, thesmall-scale flow can be simulated at the “simulation scale” l as σ ( l ) d B t = (cid:18) lL (cid:19) / σ ( L ) d B t . (3.5c)As will be shown in Section 4.2, such flow-dependent noise has a good performance in long-termsimulation, yet the drawback is that the computational costs are significantly higher comparedto the previous off-line procedure, as the SVD is computed at each time step. In this subsection, we introduce an energy conserving (in space) approximation of the abovederived stochastic system. Considering the definition of the stochastic transport operator D t defined in (2.6b) with a time increment d t q (cid:52) = q t +d t − q t defined in (2.6c), the RSW–LU systemin Eqn. (2.22a)–(2.22b) can be explicitly written asd t u = (cid:16) − u ·∇ u − f × u − g ∇ η (cid:17) d t + (cid:16) ∇· ∇· ( au ) d t − σ d B t ·∇ u (cid:17) , (3.6a)d t h = − ∇· ( u h ) d t + (cid:16) ∇· ∇· ( a h ) d t − σ d B t ·∇ h (cid:17) . (3.6b)We suggest to develop an approximation of the stochastic RSW–LU model (3.6a)–(3.6b) byfirst discretizing the deterministic model underlying this system with a structure-preservingdiscretization method (that preserves energy in space) and, then, to approximate (with a po-tentially different discretization method) the stochastic terms. Here, we use for the former avariational discretization approach on a triangular C–grid while for the latter we apply a stan-dard finite difference method. The deterministic dynamical core of our stochastic system resultsfrom simply setting σ ≈ As mentioned above, the deterministic model (or deterministic dynamical core) of the abovestochastic system results from setting σ ≈
0, which leads via (2.4) to a ≈
0. Hence, Equations(3.6a)–(3.6b) reduce to the deterministic RSW equationsd t u = (cid:16) − ( ∇ × u + f ) × u − ∇ ( 12 u ) − g ∇ η (cid:17) d t, d t h = − ∇· ( u h ) d t, (3.7)12 i T i + T i − T j T j + T j − ζ − ζ + e ij ˜ e ii − ˜ e ii + ˜ e jj − ˜ e jj + Figure 2. Notation and indexing conventions for the 2D simplicial mesh. where we used the vector calculus identity u ·∇ u = ( ∇ × u ) × u + u . Note that in thedeterministic case d t / d t agrees (in the limit d t →
0) with the partial derivative ∂/∂t . Variational discretizations.
In the following we present an energy conserving (in space)approximation of these equations using a variational discretization approach. While detailsabout the derivation can be found in Bauer and Gay-Balmaz (2019); Brecht et al. (2019), herewe only give the final, fully discrete scheme.To do so, we start with introducing the mesh and some notation. The variational discretiza-tion of (3.7) results in a scheme that corresponds to a C-grid staggering of the variables on aquasi uniform triangular grid with hexagonal/pentagonal dual mesh. Let N denote the numberof triangles used to discretize the domain. As shown in Fig. 2, we use the following notation: T denotes the primal triangle, ζ the dual hexagon/pentagon, e ij = T i ∩ T j the primal edge and˜ e ij = ζ + ∩ ζ − the associated dual edge. Furthermore, we have n ij and t ij as the normalizednormal and tangential vector relative to edge e ij at its midpoint. Moreover, D i is the discretewater depth at the circumcentre of T i , η bi the discrete bottom topography at the circumcentreof T i , and V ij = ( u · n ) ij the normal velocity at the triangle edge midpoints in the directionfrom triangle T i to T j . We denote D ij = ( D i + D j ) as the water depth averaged to the edgemidpoints.The variational discretization method does not require to define explicitly approximations ofthe differential operators because they directly result from the discrete variational principle. Itturns out that on the given mesh, these operators agree with the following definitions of standardfinite difference and finite volume operators:(Grad n F ) ij (cid:52) = F T j − F T i | ˜ e ij | , (Grad t F ) ij (cid:52) = F ζ − − F ζ + | e ij | , (Div V ) i (cid:52) = 1 | T i | (cid:88) k ∈{ j,i − ,i + } | e ik | V ik , (Curl V ) ζ (cid:52) = 1 | ζ | (cid:88) ˜ e nm ∈ ∂ζ | ˜ e nm | V nm , (3.8)for the normal velocity V ij and a scalar function F either sampled as F T i at the circumcentreof the triangle T i or sampled as F ζ ± at the centre of the dual cell ζ ± . The operators Grad n andGrad t correspond to the gradient in the normal and tangential direction, respectively, and Divto the divergence of a vector field:( ∇ F ) ij ≈ (Grad n F ) n ij + (Grad t F ) t ij , (3.9)13 ∇ · u ) i ≈ (Div V ) i , (3.10)( ∇ × u ) ζ ≈ (Curl V ) ζ . (3.11)The last Equation (3.11) defines the discrete vorticity and for later use, we also discretize thepotential vorticity as ∇ × u + fh ≈ (Curl V ) ζ + f ζ D ζ , D ζ = (cid:88) ˜ e ij ∈ ∂ζ | ζ ∩ T i || ζ | D i . (3.12) Semi-discrete RSW scheme.
With the above notation, the deterministic semi-discrete RSWequations read:d t V ij = L V ij ( V, D ) ∆ t, for all edges e ij , (3.13a)d t D i = L D i ( V, D ) ∆ t, for all cells T i , (3.13b)where L V ij and L D i denote the deterministic spatial operators, and ∆ t stands for the discrete timestep. The RHS of the momentum equation (3.13a) is given by L V ij ( V, D ) (cid:52) = − Adv(
V, D ) ij − K( V ) ij − G( D ) ij , (3.14)where Adv denotes the discretization of the advection term ( ∇ × u + f ) × u of (3.7), K theapproximation of the gradient of the kinetic energy ∇ ( u ) and G of the gradient of the heightfield g ∇ η . Explicitly, the advection term is given byAdv( V, D ) ij (cid:52) = − D ij | ˜ e ij | (cid:16) (Curl V ) ζ − + f ζ − (cid:17) (cid:18) | ζ − ∩ T i | | T i | D ji − | e ii − | V ii − + | ζ − ∩ T j | | T j | D ij − | e jj − | V jj − (cid:19) + 1 D ij | ˜ e ij | (cid:16) (Curl V ) ζ + + f ζ + (cid:17) (cid:18) | ζ + ∩ T i | | T i | D ji + | e ii + | V ii + + | ζ + ∩ T j | | T j | D ij + | e jj + | V jj + (cid:19) , (3.15)where f ζ ± is the Coriolis term evaluated at the centre of ζ ± . Moreover, the two gradient termsread: K( V ) ij (cid:52) = 12 (Grad n F ) ij , F T i = (cid:88) k ∈{ j,i − ,i + } | ˜ e ik | | e ik | ( V ik ) | T k | , (3.16)G( D ) ij (cid:52) = g (Grad n ( D + η b )) ij . (3.17)The RHS of the continuity equation (3.13b) is given by L D i ( V, D ) (cid:52) = − (cid:0) Div ( DV ) (cid:1) i , (3.18)which approximates the divergence term − ∇· ( u h ). Stabilization.
In addition, as often used in the simulations of large-scale atmospheric andoceanic flows, in order to stabilize the numerical solution (which will be important for thestochastic model), we include a biharmonic eddy viscosity with uniform coefficient µ (of unit m /s ) in the momentum equation:d t V = (cid:16) − Adv(
V, D ) ij − K( V ) ij − G( D ) ij − µL ( V ) ij (cid:17) ∆ t, (3.19)where: L ( V ) ij = (cid:0) Grad n (Div V ) ij − Grad t (Curl V ) ij (cid:1) . (3.20)14 ime scheme. For the time integrator we use a Crank-Nicolson-type scheme where we solvethe system of fully discretized non-linear momentum and continuity equations by a fixed-pointiterative method. The corresponding algorithm coincides for σ = 0 and µ = 0 with the onegiven in Section 3.2.3. The fully stochastic system has additional terms on the RHS of Equations (3.6a) and (3.6b).With these terms the discrete equations read:d t V ij = L V ij ( V, D ) ∆ t + ∆ G V ij , (3.21a)d t D i = L D i ( V, D ) ∆ t + ∆ G D i , (3.21b)where the stochastic LU-terms are given by∆ G V ij (cid:52) = (cid:16) − ∆ t (cid:0) ∇ · ∇· ( au ) (cid:1) ij + ( σ d B t ·∇ u ) ij (cid:17) · n ij , (3.21c)∆ G D i (cid:52) = − ∆ t (cid:0) ∇ · ∇· ( a D ) (cid:1) i + ( σ d B t ·∇ D ) i . (3.21d)Note that the two terms within the large bracket in (3.21c) comprise two Cartesian componentsof a vector which is then projected onto the triangle edge’s normal direction via n ij . The twoterms in (3.21d) are scalar valued at the cell circumcenters i .The parametrization of the noise described in Section 3.1 is formulated in Cartesian coordi-nates, because this allows using standard algorithms to calculate e.g. EOFs and POD. Likewise,we represent the stochastic LU-terms in Cartesian coordinates but to connect both determinis-tic and stochastic terms, we will calculate the occurring differentials with operators as providedby the deterministic dynamical core (see interface description below). Therefore, we write thesecond term in (3.21c) as( σ d B t ·∇ F ) ij = (cid:88) l =1 ( σ d B t ) lij ( ∇ F ) lij , (3.22)in which ( σ d B t ) ij denotes the discrete noise vector with two Cartesian components, constructedas described in Section 3.1 and evaluated at the edge midpoint ij . The scalar function F is aplaceholder for the Cartesian components of the velocity field u = ( u , u ). Likewise, the firstterm in (3.21c) can be written component-wise as( ∇ · ∇· ( a F )) ij = (cid:88) k,l =1 (cid:16) ∂ x k ( ∂ x l ( a kl F )) ij (cid:17) ij , (3.23)where a kl denotes the matrix elements of the variance tensor which will be evaluated, similarlyto the discrete noise vector, at the edge midpoints. For a concrete realization of the differentialson the RHS of both stochastic terms, we will use the gradient operator (3.9) as introduced next.To calculate the terms in (3.21d) we also use the representations (3.22) and (3.23) for a scalarfunction F = D describing the water depth. However, as our proposed procedure will result interms at the edge midpoint ij , we have to average them to the cell centers i .15 nterface between dynamical core and LU terms. As mentioned above, the constructionof the noise is done on a Cartesian mesh while the discretization of the deterministic dynamicalcore, corresponding to a triangular C-grid staggering, predicts the values for velocity normalto the triangle edges and for water depth at the triangle centers. We propose to exchangeinformation between the noise generation module and the dynamical core via the midpoints ofthe triangle edges where on such C-grid staggered discretizations the velocity values naturallyreside.Starting with a given predicted velocity vector with edge values V ij , we first have to recon-struct the full velocity vector field from these normal values. We use the reconstruction of thevector field in the interior of each triangle proposed by Perot et al. (2006): u i = 1 | T i | (cid:88) k = j,i − ,i + | e ik | ( x e ik − x T i ) V ik , (3.24)where x e ik are the coordinates of the edge midpoint and x T i are the coordinates of the trianglecircumcentre. By averaging values from neighboring triangles, we obtain the correspondingvalues at the edge midpoints or vertices (see Bauer (2013) for details).This reconstructed velocity vector field will be used to generate the noise as described inSection 3.1. After the noise has been constructed on the Cartesian mesh, we evaluate the discretenoise vector ( σ d B t ) ij and the discrete variance tensor ( a ) ij at the triangle edge midpoints. Thisinformation will then be used to calculate the LU noise terms in (3.21c) and (3.21d).To calculate the derivatives in these stochastic terms, we use the normal and tangentialgradient operators, i.e. the gradient operator of (3.9). To use it, we have to average values, e.g.the term ( a kl F ), to cell centers and vertices and the resulting differential will be an expressionlocated at the edge midpoint. In more detail, we can represent the partial derivative in Cartesiancoordinates by( ∂ x l F ) ij = (Grad n F ) n lij + (Grad t F ) t lij , l = 1 , . (3.25)Concretely, to discretize (3.23), we first compute ( ∂ x l ( a kl F )) ij using Equation (3.25). Thesubindex ij indicates that the resulting term is associated to the edge midpoint. To applythe second derivative in (3.23), i.e. (cid:16) ∂ x k ( ∂ x l ( a kl F )) ij (cid:17) ij , we proceed analogously, i.e. we firstaverage the terms describing the first derivative to cells and vertices and then apply once moreEquation (3.25). We proceed similarly to represent the term ∇ F in (3.22).As mentioned above, the terms in (3.21d) are calculated similarly to (3.21c) with the onlydifference that the former have to be averaged to the cell centers after the proposed discretizationprocedure has been applied. The iterated Crank-Nicolson method presented in Brecht et al. (2019) is adopted for the tem-poral discretization. Keeping the iterative solver and adding the LU terms results in an Euler–Maruyama scheme, which decrease the order of convergence of the deterministic iterative solver(see Kloeden and Platen (1992) for details). To enhance readability, we denote V t as the arrayover all edges e ij of the velocity V ij and D t as the array over all cells T i of the water depth D i at time t . The governing algorithm reads: Time-stepping algorithm
1. Start loop over k = 0 with initial guess at t : V ∗ k =0 = V t and ( D ∗ k =0 ) i = D ti + ∆ G D ij ( D t ) . Besides, we compute ∆ G V ij ( V t ). 16. Update water depth D ∗ k +1 and velocity V ∗ k +1 using explicit equation: D ∗ k +1 − D t ∆ t = − Div ( D ∗ k V ∗ k ) + Div ( D t V t )2 V ∗ k +1 − V t ∆ t = − Adv( V ∗ k , D ∗ k +1 ) + Adv( V t , D t )2 − K( V ∗ k ) + K( V t )2 − G( D ∗ k +1 ) − µL ( V ∗ k ) + ∆ G V ij ( V t )and set k + 1 = k .3. Stop loop if (cid:107) V ∗ k +1 − V ∗ k (cid:107) + (cid:107) D ∗ k +1 − D ∗ k (cid:107) < tolerance.For all simulations in this manuscript, we used a tolerance of 10 − for simulations on the planeand 10 − for simulation on the sphere.This algorithm will be used in the next section to evolve the fluid flow in time. In this section, we first study the energy behaviour of the numerical RSW–LU scheme fromabove for an inviscid test flow. Then, we show that for a viscous test flow, the stochasticmodel captures more accurately the referent structure of the large-scale flow when comparedto the deterministic model under the same coarse resolution. In addition, we demonstrate thatthe proposed RSW–LU system provides a more reliable ensemble forecast with larger spread,compared to a classical random model based on the perturbations of initial condition (PIC).
This first test case consists of two co-rotating vortices on the f -plane without viscosity (i.e. µ = 0). To illustrate the energy conservation of the spatial discretization of the RSW–LUsystem (2.22), we use the homogeneous stationary noise defined in Section 3.1.1 since the twoincompressible constraints ∇· σ d B t = 0 and ∇· ∇· a = 0 in (2.22d) are naturally satisfied.Then, no extra steps are required to satisfy the incompressible constraints. Initial conditions
The simulation is performed on a rectangular double periodic domain Ω = [0 , L x ] × [0 , L y ] with L x = 5000 km and L y = 4330 km, which is discretized into N = 32768 triangles. The large-scaleflow is assumed to be under a geostrophic regime at the initial state, i.e. f k × u = − g ∇ h . Weuse an initial height field elevation (as e.g. in Bauer and Gay-Balmaz (2019)) of the form h (cid:0) x, y, t = 0 (cid:1) = H − H (cid:48) (cid:18) exp (cid:16) − x (cid:48) + y (cid:48) (cid:17) + exp (cid:16) − x (cid:48) + y (cid:48) (cid:17) − πs x s y L x L y (cid:19) , (4.1a)where the background height H is set to 10 km, the magnitude of the small perturbed height H (cid:48) is set to 75 m and the periodic extensions are given by x (cid:48) i = L x πs x sin (cid:0) πL x ( x − x c i ) (cid:1) , y (cid:48) i = L y πs y sin (cid:0) πL y ( y − y c i ) (cid:1) , i = 1 , x c , y c ) = ( L x , L y ), ( x c , y c ) = ( L x , L y ) withparameters ( s x , s y ) = ( L x , L y ). To obtain the discrete initial water depth D i , we sample the17 igure 3. Contour plots of the potential vorticity fields after 2 days for (left) one realization of a LU simulationwith homogeneous noise and (right) a deterministic run. The contour interval is 0.4 days − km − . analytical function h at each cell centre. Subsequently, the discrete geostrophic velocities ateach triangle edge ij at the initial state can be deduced via V ij = − gf (Grad t D ) ij , (4.2)where the Coriolis parameter f is set to 5 . − . For the LU simulations, the magnitude ofthe homogeneous noise remains moderate with its constant variance a set to be 169 . · s − . Analysis of energy conservation
To analyze the energy conservation properties of our stochastic integrator, we use the aboveinitial conditions to simulate the two co-rotating vortices for 2 days. In Figure 3, we showcontour plots of the potential vorticity (as defined in (3.12)) fields of the deterministic andstochastic models. We observe that under the moderate noise with a as chosen above, thelarge-scale structure of the stochastic system is similar to that of the deterministic run.On the specific staggered grid as shown in Figure 2, the total energy of the shallow waterequations (2.16) for both deterministic and stochastic case is approximated byE( t ) ≈ N (cid:88) i =1 D i ( t ) | T i | (cid:88) k = j,i − ,i + | T i | h ik f ik (cid:0) V ik ( t ) (cid:1) + 12 g (cid:0) D i ( t ) (cid:1) | T i | . (4.3)As shown in Bauer and Gay-Balmaz (2019), the proposed discrete variational integrator (seeSection 3.2.1) together with an iterative Crank-Nicolson time stepping method exhibits a 1storder convergence rate of the energy error with smaller time step size. This will allows usimmediately to simply include the stochastic terms to result in an Euler-Maruyama type timeintegrator for stochastic systems (cf. Section 3.2.2).In the present work, we consider the energy behavior of the deterministic scheme (i.e. the vari-ational integrator) as reference, which is denoted as E REF ( t ) in the following. For the stochasticRSW model, the Euler-Maruyama time scheme might lead to a different behavior with respectto energy conservation when compared to the deterministic model. In order to quantify numeri-cally the energy conservation of the RSW–LU, we propose to measure the relative errors between18 days -16 -14 -12 -10 -8 -6 r e l a t i v e e rr o r (a) Evolution of the relative L errors betweenthe energy of the mean RSW–LU and the refer-ence, using ∆ t (blue line), ∆ t/
10 (red line) and∆ t/
100 (yellow line) respectively. -5 -4 t -9 -8 -7 (b) Convergence of the energy path of theRSW–LU to that of the reference w.r.t. timestep sizes. The blue line shows the global er-rors of the ensemble mean energy, the blue areadescribes the 68% confident interval of the en-semble errors and the dashed line stands for the1st order convergence rate. Figure 4. Analysis of the numerical energy conservation of the RSW–LU. the mean stochastic energy, denoted as E LU ( t ), and the reference E REF ( t ) by E LU ( t ) / E REF ( t ) − t results in smaller relative errors.To determine more quantitatively the convergence rate of the stochastic scheme (relative tothe reference) with respect to different time step sizes, we defined the following global (in spaceand time) error measure: ε (E LU ) (cid:52) = (cid:107) E LU ( t ) − E REF ( t ) (cid:107) L ,T ]) (cid:107) E REF ( t ) (cid:107) L ,T ]) , (4.4)where (cid:107) f ( t ) (cid:107) L ,T ]) = ( (cid:82) T | f ( t ) | d t ) / and T is set to 2 days. We determine for an ensemblewith 10 members such global errors in order to illustrate the convergence rate of each ensemblemember and the spread between those rates. This spread is illustrated as blue shaded area inFigure 4b. The area centre is determined by the mean of the errors, and the dispersion of thisarea is given by one standard derivation ( i.e.
68% confident interval of the ensemble of ε (E LU )).Besides, the minimal and maximal values of the errors of the ensemble are represented by thevertical bar-plots. The blue line of Figure 4b shows that the convergence rate (w.r.t. various∆ t ) of the ensemble mean energy is of 1st order. This is consistent with the weak convergencerate of order O (∆ t ) of the Euler-Maruyama scheme, cf. Section 3.2.3. Next, we want to show that our stochastic system better captures the structure of a large-scale flow than a comparable deterministic model. To this end, we use a viscous test case andheterogeneous noise.The viscous test case we use is proposed by Galewsky et al. (2004) and it consists of abarotropically unstable jet at the mid-latitude on the sphere. This strongly non-linear flow willbe destabilized by a small perturbation of the initial field, which induces decaying turbulenceafter a few days. However, the development of the barotropic instability in numerical simulations19ighly depends on accurately resolving the small-scale flow, which is particularly challenging forcoarse-grid simulations. For the same reason, the performance of an ensemble forecast systemin this test case is quite sensible to the numerical resolution. In the following, we demonstratethat the RSW–LU simulation on a coarse mesh under heterogeneous noises, provides betterprediction of the barotropic instability compared to the deterministic coarse simulation, andproduces more reliable ensemble spread than the classical PIC simulation.
Initial conditions
The values of the principle parameters for the simulations are specified in Table 1. Under thegeostrophic regime, the initial zonal velocity and height is respectively given by u (Θ , t = 0) = U e n exp (cid:16) − Θ )(Θ − Θ ) (cid:17) , for Θ < Θ < Θ , (4.5a) h (Θ , t = 0) = H − Rg (cid:90) Θ u ( θ, t = 0) (cid:16) Ω sin θ + tan θR u ( θ, t = 0) (cid:17) d θ, (4.5b)where e n = exp (cid:0) − / (Θ − Θ ) (cid:1) is used to rescale the jet magnitude to the maximal value U at the jet’s mid-point Θ = π/
4. As introduced by Galewsky et al. (2004), in order to initiatethe barotropic instability, the following localized bump is included in the height field: h (cid:48) (Υ , Θ) = H (cid:48) cos Θ exp (cid:16) − (3Υ) − (cid:0) π − Θ) (cid:1) (cid:17) , (4.5c)where Υ denotes the longitude. Analogously to the previous inviscid test case, we then use theseanalytic functions (4.5) to sample the discrete velocity at the edge mid-point and the height fieldat the cell centre on the staggered mesh (See Figure 2).Parameters Value Description(Θ , Θ ) (2 π, π ) /
14 rad Initial latitude limits H .
158 km Background height H (cid:48)
120 m Initial perturbation amplitude R . × km Mean radius of Earth g .
806 m · s − Gravity of Earth˜ Ω . × − s − Angular rotation rate of Earth U
80 m · s − Maximum zonal velocity µ l . × m · s − Fine-grid biharmonic viscosity µ L . × m · s − Coarse-grid biharmonic viscosity∆ t l
12 s Fine-grid time step∆ t L
50 s Coarse-grid time step N l N L Table 1. Parameter list for simulations of the barotropic instability.
For the LU simulations, we use the two heterogeneous noises described in Section 3.1.2, basedon either the off-line learning of EOFs from the high-resolution simulation data, denoted as LUoff-line, or on the on-line estimation of EOFs from the coarse-grid simulation, denoted as
LUon-line . To allow for comparisons, the strength of these two noises are imposed to be the same.The PIC stochastic model is obtained as follows: first, we perform ensemble simulations of theLU off-line and the LU on-line method over 1 day. Then, each ensemble realization is used as oneinitial random state for the
PIC off-line and the
PIC on-line simulations, respectively. For each20 igure 5. Snapshots of the vorticity field on the sphere for different models (with 20480 triangles) after 5 days.From left to right: reference, ensemble mean of LU online and deterministic LR. stochastic model, an ensemble run with 20 realizations is done. Besides, a deterministic coarse-grid simulation, denoted as LR , is also performed. For all these coarse models, the biharmonicviscosity coefficient is fixed to be the same as given in Table 1. Prediction of barotropic instability
In this section, we compare the predictions of the barotropic instability for different coarse mod-els to that provided by the reference simulation. The latter is obtained from the coarse-grainingprocedure through a bilinear interpolation of the high-resolution snapshots. In Figure 5, weillustrate snapshots of the vorticity fields on the sphere for the reference, LU and deterministicmodels after a simulation time of 5 days. We can clearly see that at that day the LU ensem-ble mean better captures the large-scale structure of the reference flow than the deterministicsimulation. To better distinguish the differences in the simulations, contour plots of the vortic-ity fields at day 4, 5 and 6, localized at the mid-latitude of the sphere, are given in Figure 6.From the evolution of the reference vorticity fields, we observe that the barotropic instabilityof the mid-latitude jet starts to develop at day 4. Subsequently, more and more small-scalefeatures emerge and the flow becomes turbulent. Furthermore, both LU on-line and LU off-linesimulations exhibit the stretched out wave at day 5 in the same way as the reference does, andthat some big vortices start to separate from the wave at day 6. On the other hand, thesecharacteristics are not correctly captured in both PIC off-line and LR simulations. We remarkthat the results of PIC on-line simulations are not include in Figure 6, since they behave quitesimilarly to the PIC off-line run.To physically interpret the above results, it is useful to analyze the energy spectra of differentmodels. From a basic knowledge of the two-dimensional turbulence theory (McWilliams, 2006),the potential enstrophy is transferred from the large scales to the small scales by the directcascade, whereas the kinetic energy is transferred from the small scales to the large scales bythe inverse cascade. However, introducing only a dissipation mechanism for coarse models oftenleads to an excessive decrease of the resolved kinetic energy (Arbic et al., 2013; Kjellsson andZanna, 2017). In our test case, this kind of issue is present in both PIC and the LR simulations,where the small-scale energy and enstrophy are over-dissipated, as illustrated in Figure 7. On theother hand, introducing the non-linear convection by the noise, the LU dynamical systems bringhigher turbulent energy and enstrophy to the small scales, which leads to better structuring21 ay 4 Day 5 Day 6 -3 -2 -1 0 1 -101 -4 -3 -2 -1 0 1 Figure 6. Comparison of the vorticity contour plots along the mid-latitude jet for different models (with 20480triangles) at day 4, 5 and 6 respectively. From top to bottom: reference, ensemble mean of LU on-line, ensemblemean of LU off-line, ensemble mean of PIC off-line and deterministic LR. The contour interval is fixed to 2 × − s − , the x-axis is longitude (in rad) and the y-axis is latitude (in rad).
22f the large-scale flow. For instance, the ensemble mean of the energy and enstrophy spectrafor both LU on-line and LU off-line simulations are much closer to that of the references atdifferent days. Note that these spectra on the sphere are calculated using the method proposedby Aechtner et al. (2015): first, the energy and enstrophy is interpolated onto a Gaussian grid,then the spherical harmonics basis are used to compute the power spectral density.
Kinetic energy Normalized enstrophy -4 -2 -4 -2 -4 -2 Figure 7. Comparison of the ensemble mean of the kinetic energy (left column) spectrums and the potentialenstrophy (right column) spectrums for different models (with 20480 triangles) at day 5 (1st row), 7 (2nd row)and 10 (3rd row) respectively. Note that the potential enstrophy is defined by the square of the potential vorticityand each potential enstrophy spectrum is normalized by its first value at the largest wavenumber. The dashedline is the k − (left column) and k − / (right column) power law. Evaluation of ensemble forecasts
Once the ensembles have been produced by the random models, we measure the reliability of theensemble forecast systems by some simple metrics. But before we do so, let us first demonstratequalitatively the time evolution of each ensemble spread and compare it with the observationtrajectory. To determine the latter, we evaluate the local vorticity field of the reference atdifferent grid points in the region of the mid-latitude jet. These points serve as observationpoints. The evolution of the spread of the ensemble forecast systems is then build by the 95%confident interval of its ensemble trajectories at each selected point. As shown in Figure 8,for the six local points chosen along the longitude Υ = − .
53 rad, the ensemble spreads ofthe LU off-line system are large enough to almost always include the observation trajectories,whereas the spreads of the PIC off-line system are quite small so that the observations are23ot always contained within the spread. For the latter, this will result in a wrong coupling ofthe measurement and the ensemble system, when performing data assimilation (Gottwald andHarlim, 2013; Franzke et al., 2015). -101 -4 -101 -4 -4 Figure 8. Comparison of the ensemble spread evolution over 20 days of the vorticity field for the LU-offline (redarea) runs and the PIC-offline (blue area) runs, at six different locations Θ = (0 . , . , . , . , . , .
2) radalong the longitude Υ = − .
53 rad. The observation trajectories are shown by the black lines.
To quantify whether the ensemble spread of the forecast system represents the true uncer-tainty of the observations, the rank histogram (Talagrand et al., 1997; Hamill, 2001) is widelyadopted as a diagnostic tool. This approach checks where the verifying observation usually fallsw.r.t. the ensemble forecast states which are arranged in an increasing order at each grid point.In an ensemble with perfect spread, each member represents an equally likely scenario, so theobservation is equally likely to fall between any two members. To construct the rank histogramin our test case, we proceed as follows:1. At every grid point x i , we rank the N e vorticity values { q ( j ) ( x i ) } j =1 ,...,Ne of the ensemblefrom lowest to highest. This results in N e + 1 possible bins which the observations can fallinto, including the two extremes;2. Identify which bin the observation vorticity q o ( x i ) falls into at each point x i ;3. Tally over all observations { q o ( x i ) } i =1 ,...,No to create a histogram of rank.As shown in Figure 9, the histograms of both random models exhibit a U-shape for a few days inthe beginning, while after a simulation time of about 10 days, the histograms of both LU on-lineand LU off-line systems become mostly flat. A U-shape indicates that the ensemble spread istoo small so that many observations are falling outside of the extremes of the ensemble while adome-shape indicates the contrary. A flat histogram, in contrast, indicates that the ensemblemembers and observations are sampled from a common distribution. We observe that the LUoff-line system performs slightly better than the LU on-line version. In contrast to these verygood ensemble spreads, the histograms of both PIC on-line and PIC off-line systems remainin a U-shape during the entire simulation period which indicates that these systems do notaccurately estimate the correct uncertainty around the observations.It is important to notice that a flat rank histogram does not necessarily imply good fore-casts, it only measures whether the observed probability distribution is well represented by the24 ay 5 Day 10 Day 15 Day 20
10 20
10 20
10 20
Figure 9. Comparison of the rank histograms for the LU on-line (1st row) runs, the LU off-line (2nd row) runs,the PIC on-line (3rd row) runs and PIC off-line (last row) runs, at day 5, 10, 15 and 20 respectively. i.e.
MSE ( t ) = 1 N o N o (cid:88) i =1 (cid:0) q o − (cid:98) E [ q ] (cid:1) ( t, x i ) ≈ (cid:16) N e + 1 N e (cid:17) N o N o (cid:88) i =1 (cid:100) Var[ q ]( t, x i ) = N e + 1 N e MEV ( t ) , (4.6)where (cid:98) E [ q ] = N e (cid:80) N e j =1 q ( j ) and (cid:100) Var[ q ] = N e − (cid:80) N e j =1 (cid:0) q ( j ) − (cid:98) E [ q ] (cid:1) denote the empirical mean andthe empirical variance, respectively.In Figure 10, we compare the differences in time between the MSE and the MEV, normalizedby the squared maximum of the initial vorticity, for the different random models from above.From these curves we can deduce that the LU off-line system exhibits the lowest errors duringthe entire simulation time of 20 days. In particular, during the first 10 days, these errors aresignificantly lower when compared to the other models, which can be explained by the fact thatthe LU off-line system incorporates data from the reference into the ensemble, which increasesthe reliability of the ensemble forecast. Although the errors between MSE and MEV of the LUon-line system is larger than the LU offline system from day 5 to day 10, they remain at low levelfrom day 10 onwards, implying that the reliability of the former increases for longer simulationtimes. In contrast, both PIC off-line and PIC on-line systems show higher error values at mostof the times and hence provide less reliable ensembles. We remark that other metrics, such asthe continuous ranked probability score (Resseguier et al., 2020; Weigel, 2012), can also be usedto measure a calibrated ensemble. days ( M SE - M EV ) / q Figure 10. Comparison of the differences between the mean square error (MSE) and the mean ensemble variance(MEV) of the ensemble vorticity fields for the LU on-line (red dashed line) runs, the LU off-line (red solid line)runs, the PIC on-line (blue dashed line) runs and the PIC off-line (blue solid line) runs. Note that these differencesare normalized by q = (cid:107) q (Υ , Θ , t = 0) (cid:107) ∞ . In this study, we introduced a stochastic version of the rotating shallow water equations underlocation uncertainty (RSW-LU). The derivation is based on a stochastic Reynolds transport26heorem, where the fluid flow is decomposed into a large-scale component and a noise termmodelling the unresolved small-scale flow. A benefit of this approach is that the total energyis conserved along time for any realization. In order to preserve this structure, we combinedan energy (in space) preserving discretization of the underlying deterministic equations of thisRSW–LU system with approximations of the stochastic terms that are based on standard finitevolume/difference operators.We could show for an f-plane test case that this approach leads for homogeneous noise to adiscretization of the RSW-LU system that preserves (spatially) the total energy. Moreover, usinginhomogeneous noise that well captures the impact of small scales to the large-scale flow, wedemonstrated that for a barotropically unstable jet on the sphere our proposed RSW–LU modelbetter predicts the development of the instabilities than a comparable deterministic model, whilethe ensemble spread of the RSW–LU system is more likely to contain the observations comparedto an ensemble of deterministic simulations with perturbed initial conditions (PIC). We alsoshowed that the RSW–LU forecast systems follows a common distribution of the observationsand is more reliable than the PIC system.Showing accurate ensemble spreads and reliable uncertainty quantification, we will next applyour developed RSW-LU system to data assimilation. We will also work towards discretizations ofstochastic flow models in the framework of LU that preserve total energy both in space and timeto which the present work provides a first step. Exploiting the modular approach of combiningdifferent discretizations for deterministic and stochastic terms, in future work we will explorethe possibility to consistently extend existing atmospheric and ocean models with stochasticparametrizations.
Acknowledgments
The authors acknowledge the support of the Mitacs Globalink Research Award and of theERC EU project 856408-STUOD. Besides, we would like to thank Alexander Bihlo and ScottMacLachlan for helpful discussions and thank Matthias Achtner for providing code to computethe energy spectrum on the sphere.
References
M. Aechtner, N. K.-R. Kevlahan, and T. Dubos. A conservative adaptive wavelet method for theshallow-water equations on the sphere.
Quarterly Journal of the Royal Meteorological Society ,141(690):1712–1726, 2015. doi: 10.1002/qj.2473.J. Anderson and S. Anderson. A Monte Carlo implementation of the nonlinear filtering problemto produce ensemble assimilations and forecasts.
Monthly Weather Review , 127(12):2741–2758,1999.B. K. Arbic, K. L. Polzin, R. B. Scott, J. G. Richman, and J. F. Shriver. On eddy viscosity,energy cascades, and the horizonal resolution of gridded stallite altimeter products.
Journalof Physical Oceanography , 43(2):283–300, 2013.W. Bauer and F. Gay-Balmaz. Towards a geometric variational discretization of compressiblefluids: the rotating shallow water equations.
Journal of Computational Dynamics , 6:1, 2019.W. Bauer, P. Chandramouli, B. Chapron, L. Li, and E. M´emin. Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach.
Journal ofPhysical Oceanography , 50(4):983–1003, 2020a.W. Bauer, P. Chandramouli, L. Li, and E. M´emin. Stochastic representation of mesoscale eddyeffects in coarse-resolution barotropic models.
Ocean Modelling , 151:101646, 2020b.27erner Bauer.
Toward goal-oriented R-adaptive models in geophysical fluid dynamics using ageneralized discretization approach . PhD thesis, Hamburg University Hamburg, 2013.R. Brecht, W. Bauer, A. Bihlo, F. Gay-Balmaz, and S. MacLachlan. Variational integrator for therotating shallow-water equations on the sphere.
Quarterly Journal of the Royal MeteorologicalSociety , 145(720):1070–1088, 2019.B. Chapron, P. D´erian, E. M´emin, and V. Resseguier. Large-scale flows under location un-certainty: a consistent stochastic framework.
Quarterly Journal of the Royal MeteorologicalSociety , 144(710):251–260, 2018.G. Da Prato and J. Zabczyk.
Stochastic equations in infinite dimensions . Encyclopedia ofMathematics and its Applications. Cambridge University Press, 2 edition, 2014.C. E. Franzke and A. J. Majda. Low-order stochastic mode reduction for a prototype atmosphericGCM.
Journal of the Atmospheric Sciences , 63(2):457–479, 2006.C. E. Franzke, T. J. O’Kane, J. Berner, P. D. Williams, and V. Lucarini. Stochastic climatetheory and modeling.
Wiley Interdisciplinary Reviews: Climate Change , 6(1):63–78, 2015.J. S. Frederiksen, T. J. O’Kane, and M. J. Zidikheri. Subgrid modelling for geophysical flows.
Philosophical Transactions of the Royal Society A: Mathematical, Physical and EngineeringSciences , 371(1982):20120166, 2013.J. Galewsky, R. K. Scott, and L. M. Polvani. An initial-value problem for testing numerical mod-els of the global shallow-water equations.
Tellus A: Dynamic Meteorology and Oceanography ,56(5):429–440, 2004.G. Gottwald and J. Harlim. The role of additive and multiplicative noise in filtering complex dy-namical systems.
Proceedings of the Royal Society A: Mathematical, Physical and EngineeringScience , 469(2155):20130096, 2013.G. Gottwald, D. T. Crommelin, and C. E. Franzke. Stochastic climate theory. In
Nonlinear andStochastic Climate Dynamics , pages 209–240. Cambridge University Press, 2017.I. Grooms and A. J. Majda. Stochastic superparameterization in quasigeostrophic turbulence.
Journal of Computational Physics , 271:78–98, 2014.Ernst Hairer, Christian Lubich, and Gerhard Wanner.
Geometric numerical integration:structure-preserving algorithms for ordinary differential equations , volume 31. Springer Sci-ence & Business Media, 2006.T. M. Hamill. Interpretation of rank histograms for verifying ensemble forecasts.
MonthlyWeather Review , 129:550–560, 2001.S. Kadri Harouna and E. M´emin. Stochastic representation of the Reynolds transport theorem:revisiting large-scale modeling.
Computers and Fluids , 156:456–469, 2017.J. Kjellsson and L. Zanna. The impact of horizontal resolution on energy transfers in globalocean models.
Fluids , 2(3):45, 2017.P. E. Kloeden and E. Platen.
Numerical Solution of Stochastic Differential Equations , volume 23.Springer-Verlag Berlin Heidelberg, 1992.H. Kunita.
Stochastic flows and stochastic differential equations , volume 24 of
Cambridge Studiesin Advanced Mathematics . Cambridge University Press, 1997.28. Majda, C. Franzke, and B. Khouider. An applied mathematics perspective on stochastic mod-elling for climate.
Philosophical Transactions of the Royal Society of London A: Mathematical,Physical and Engineering Sciences , 366(1875):2427–2453, 2008.Jerrold E Marsden and Matthew West. Discrete mechanics and variational integrators.
ActaNumerica , 10(1):357–514, 2001.J. C. McWilliams.
Fundamentals of Geophysical Fluid Dynamics . Cambridge University Press,2006.E. M´emin. Fluid flow dynamics under location uncertainty.
Geophysical & Astrophysical FluidDynamics , 108(2):119–146, 2014.J Blair Perot, Dragan Vidovic, and Pieter Wesseling. Mimetic reconstruction of vectors. In
Compatible Spatial Discretizations , pages 173–188. Springer, 2006.S. Pope.
Turbulent flows . Cambridge University Press, 2000.V. Resseguier, E. M´emin, and B. Chapron. Geophysical flows under location uncertainty, partI: Random transport and general models.
Geophysical & Astrophysical Fluid Dynamics , 111(3):149–176, 2017a.V. Resseguier, E. M´emin, and B. Chapron. Geophysical flows under location uncertainty, partII: Quasi-geostrophic models and efficient ensemble spreading.
Geophysical & AstrophysicalFluid Dynamics , 111(3):177–208, 2017b.V. Resseguier, E. M´emin, and B. Chapron. Geophysical flows under location uncertainty, partIII: SQG and frontal dynamics under strong turbulence.
Geophysical & Astrophysical FluidDynamics , 111(3):209–227, 2017c.V. Resseguier, L. Li, G. Jouan, P. Derian, E. M´emin, and B. Chapron. New trends in ensembleforecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics.
Archives of Computational Methods in Engineering , pages 1886–1784, 2020.L. Sirovich. Turbulence and the dynamics of coherent structures, part I: Coherent structures.
Quarterly of Applied Mathematics , 45(3):561–571, 1987.O. Talagrand, R. Vautard, and B. Strauss. Evaluation of probabilistic prediction systems.Workshop on Predictability, ECMWF, 1997.G. K. Vallis.
Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation .Cambridge University Press, 2 edition, 2017.A. P. Weigel. Ensemble forecasts. In