{'fraction_non_alphanumeric': 0.0580625, 'fraction_numerical': 0.02609375, 'mean_word_length': 3.9529484599907136, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical weather prediction in two dimensions with topography, using a finite volume method

We aim to study a finite volume scheme to solve the two dimensional inviscid primitive equations of the atmosphere with humidity and saturation, in presence of topography and subject to physically plausible boundary conditions to the system of equations. In that respect, a version of a projection method is introduced to enforce the compatibility condition on the horizontal velocity field, which comes from the boundary conditions. The resulting scheme allows for a significant reduction of the errors near the topography when compared to more standard finite volume schemes. In the numerical simulations, we first present the associated good convergence results that are satisfied by the solutions simulated by our scheme when compared to particular analytic solutions. We then report on numerical experiments using realistic parameters. Finally, the effects of a random small-scale forcing on the velocity equation is numerically investigated. The numerical results show that such a forcing is responsible for recurrent large-scale patterns to emerge in the temperature and velocity fields.


Introduction
We consider the inviscid primitive equations of the atmosphere with humidity. Besides its mathematical importance, the inviscid primitive equations are used for numerical weather predictions in geophysics. Considerable effort and attention have been devoted to this subject over the last few decades. The theory of the inviscid primitive equations usually does not resemble the theory of the Euler equations. In addition, it is well-known that the inviscid primitive equations are ill-posed for any set of boundary conditions of local type; see e.g. [27] and [37]. Hence the theoretical and numerical understanding of this topic is very scarce and remains as an important open problem.
In the presence of the topography and the divergence free term in the primitive equations (conservation of mass), the numerical methods require a careful adaptation. We propose a suitable finite volume scheme to overcome these difficulties. In addition, we consider a compatibility condition which comes from the boundary conditions. To enforce the compatibility condition, we introduce a version of a projection method. We perform various numerical simulations which include deterministic and stochastic cases. The aim of this article is to see the influence of the mountain; starting with an unsaturated humidity, we observe that the rain appears near the mountain. In addition, by incorporating an additive noise to the model formulation, it is numerically illustrated that recurrent large-scale patterns can emerge from the combined effect of a random small-scale forcing with those of the topography.
Since Lions, Temam, and Wang proposed new formulations of the primitive equations in [25], the mathematical studies of primitive equations have been developed in many different directions. Considering the viscosity, one can find many mathematical results for the primitive equations in e.g. [3], [22], [23], and [30]. In the absence of viscosity, two of authors of this article have studied these equations with a set of nonlocal boundary conditions in both theoretical and computational sides; see e.g. [6], [7], [8], [9], [33] and [34]. The primitive equations with humidity have been investigated in the classical references [19], [20], [21], and [32]. The authors from [11] and [12] have proposed the problem of water vapor in presence of saturation for the simplified model. We focus, in this article, on a numerical method for the realistic and complex model of the inviscid primitive equations.
This article is organized as follows. In Section 2 we present our model equations and physically plausible boundary conditions. We also discuss our projection method for the velocity, and how the pressure relates to the topography. Then in Section 3, we present a finite volume method to solve numerically the model equations. Due to the topography, the classical finite volume schemes produce errors near the topography. To resolve this problem, we propose a new scheme which is a modified Godunov type method that exploits the discrete finite-volume derivatives by using the so-called Taylor Series Expansion Scheme (TSES) introduced in [2] and [18].
Finally, in Section 4, we report on numerical results based on the scheme introduced in Section 3, in the deterministic as well as stochastic context. In the deterministic setting, it is shown for physically plausible parameter values, how the projection method proposed in Section 2.3 allows for the compatibility condition (2.12)-(2.13) -that the vertical integration of the horizontal component of the model's velocity field must satisfy -to be satisfied to a better numerical accuracy compared to when the projection method is not used. The effects of a small-scale additive noise on the (horizontal component) of the velocity equation is then numerically investigated. As a main result, it is shown that such a small-scale random forcing can significantly impact the model's dynamics at the large scales, leading to the appearance of waves that although evolving irregularly in time, manifest characteristic frequencies across a low-frequency band which is more pronounced in the temperature field than in the velocity field.

The primitive equations with humidity and saturation
Our goal in this Section is to introduce the primitive equations of the atmosphere with humidity and saturation, then describe the boundary conditions of the problems.
The two dimensional inviscid model, in coordinates (x, p), under consideration accounts for the conservation of horizontal momentum (in one direction), conservation of mass and energy. The hydrostatic equation is introduced as well as the equation of conservation of water vapor with saturation; see e.g. [1], [11], [20], [21], [28], [32], [25], and [30]. In the simple humidity model that we consider, following [20], [21], and [32], the water vapor leaves the system when it condenses.
2.1. The model equations. We consider the two-dimensional inviscid primitive equations which depends on two spatial like variables, x ∈ [0, L] and p ∈ [p A , p B ], the pressure. We denote by M the (pseudo) spatial domain M = (0, L) × (p A , p B ); the function p B = p B (x, t) refers to the pressure at the bottom of the atmosphere (if we do not have topography, we simply set p B (x, t) = p 0 = 1000), and p A refers to the pressure at the top of the atmosphere. We choose the value p A = 200. The two-dimensional inviscid primitive equations then read: where H is the Heaviside function H(x) = 1 2 (1 + sign(x)), based on equation (9.14) in [20]; see also [12], [11]. Note that in the first equation ( [19]). • R = 287 J K −1 kg −1 is the gas constant for dry air (p. 597 in [19]).
• q s (T, p) is the saturation specific humidity. From equation (9.6) in [20] we have: where e s (T ) is the saturation vapor pressure. We approximate its value with equation ( We aim to find: • T = T (x, p, t): local temperature.
• u = u(x, p, t): the velocity along the x axis.
The variables ω is a diagnostic variables which will be computed using the prognostic variables u. We treat the geopotential φ = φ(x, p, t) separately using (2.1) 5 .
The boundary conditions. We supplement these equations with the physically relevant boundary conditions. At the top (p = p A ) and bottom (p = p B ) we consider the impermeable boundary conditions: and we do not assign any boundary condition for T and q, at At x = 0 we consider Dirichlet boundary conditions for T , q and u. For ω we impose homogeneous Neumann boundary condition. At x = L we consider the homogeneous Neumann boundary conditions for T , q, u and ω, that is, where g T , g q , and g u are sufficiently smooth functions defined on [p A , p B (0)]. We also assume the following conditions on z B and p B to make our numerical computations easier: This means that the topography is flat near x = 0 and L (see e.g. Figure 2.1). The boundary conditions (2.5) and (2.6) plays an important role in the computation of ω from u. Indeed, with the boundary condition ω(p A ) = 0 and (2.1) 4 , we obtain Remark 2.1. We assume that ω = 0 at the top (p = p A ), but, in general, ω(p B ) = 0 in the presence of topography. We have two boundary conditions for ω: at top (2.5) and at the bottom (2.6). We directly enforce the boundary condition at the top by computing ω in (2.10). We also need to verify that the equation (2.10) at p = p B will satisfy (2.6), for this reason we propose a compatibility condition for u. We first rewrite equation (2.6) as In view of (2.11) and (2.10) we deduce that ∂u ∂x dp + ω(p B ) = 0 (by (2.10)). (2.12) Thus, we see that u satisfies the following compatibility condition: We recall our system of equations and we interpret (2.1) 6 as an equation for z and (2.14) φ(x, p, t) = z(x, p, t)g.
If we know φ and one of the characteristics of the topography, either z B (x) or p B (x), we can find the other one using (2.14). For the computation of φ x , we differentiate in x equation (2.1) 5 and we obtain then, from the boundary condition (2.7), we deduce

2.2.
Reformulation of the equations. In this subsection we rewrite and simplify (2.1) in view of the numerical simulations in Sections 3 and 4. Thanks to the divergence free condition (2.1) 4 , the equations of T , q and u in (2.1) 1 -(2.1) 3 read: We assume that T and φ are perturbations of a stratified configuration, (T (p),φ(p)) satisfying the hydrostatic equation and this implies that p B does not depend on t as announced, that is: For the sake of simplicity, we set (2.21) In view of of the notations above, (2.10), (2.16), and (2.17) we arrive at the following boundary value problem: and the boundary conditions are the same as in (2.8): (2.23) 2.3. Projection of u. In this section, we develop a projection method to ensure the compatibility condition (2.13) in Remark 2.1; as we will see, this projection method is similar -but simpler -than the projection method in incompressible fluid mechanics Temam [36] and Chorin [10]. Letũ ∈ L 2 (M) be the solution in (2.22) 3 that does not satisfy (2.13). We now construct its projection u = Pũ on the appropriate set of functions that satisfy the compatibility condition (2.13): We denote by (I − P )L 2 the orthogonal complement of P L 2 in L 2 (M).
Proposition 2.1. The orthogonal complement of P L 2 in L 2 (M) can be characterized as follows Proof. We first show that (I −P )L 2 ⊃ R.H.S of (2.25). Let α = α(x) ∈ L 2 ([0, L]) be independent of p and L 0 α(x)dx = 0. Letᾱ be the primitive function of α vanishing at x = 0, that is, so that α =ᾱ x andᾱ ∈ H 1 (0, L). Note also thatᾱ(L) = 0 since L 0 α(x)dx = 0. Then for v ∈ P L 2 , we deduce that Conversely, we prove that (I − P )L 2 ⊂ R.H.S of (2.25). Letũ ∈ L 2 (M) and let u be its orthogonal projection on P L 2 , thenũ − u = (I − P )ũ and so ∂φ/∂p ∈ P L 2 , which gives This implies that ∂ ∂p (u −ũ) = 0 in the sense of distributions on M. Therefore u −ũ does not depend on p, that is, u −ũ = α ∈ L 2 (0, L) (see [35] and Section 4.4 in [30] for more details about distributions independent of one variable). It remains to prove that For the last implication we have chosen an arbitrary v ∈ P L 2 such that p B p A vdpwhich is constant in x and independent of p -is not zero. Thus, we have and (2.25) is proved. Now we assume that φ and ω are known and denote the third component of the solution of equation (2.22) 1 byũ ∈ L 2 (M); so far ω is not required yet to satisfy (2.22) 2 . Let u ∈ P L 2 be the orthogonal projection ofũ onto P L 2 and let λ x (x) be its orthogonal complement so that λ, λ x ∈ L 2 (0, L). Hence, we obtain We integrate (2.27) in p from p A to p B and differentiate in x using (2.13): p Aũ dp.
From (2.28) we deduce the equation for λ x : p Aũ dp.
2.4. Treatment of p B (x). As explained in Section 2.1, the topography is determined by either p B (x) or z B (x). To find the relation between z B and p B , we use (2.1) 6 and (2.18): then by integrating in p At z = 0, we have p = p 0 (virtual pressure on the whole segment) and this gives C: Moreover, from the fact that z = z B at p = p B , we deduce that Since the topography is known, that is the function z B = z B (x) is given, and we can then compute p B (x) from (2.32) using the Newton method. However, in our simplified calculations, we choose p B (x) and deduced the topography from (2.32) to avoid the repeated use of the Newton method.
3. Numerical scheme: the finite volume method In this article, we use the Godunov's method as described in Chapter 23 of [24] combined with a spatial discretization by finite volumes which has to be performed with a special care here due to the topography. In the presence of topography, the spatial domain in x and p is not rectangular. Such a geometry gives rise to computational difficulties since the classical methods for the directional derivatives are not accurate. To avoid this problem, we first propose a specific discretization in a given spatial domain in Section 3.1. Then, we introduce the first order finite volume scheme to compute T , q andũ in Section 3.2. We then study the discrete projection method, in Section 3.3, to obtain the solution u. Following the computation of the diagnostic variables, we consider the computation of ω, in Section 3.4, and the computation of φ in Section 3.5, the two prognostic variables. Finally, in Section 3.6 we present our discretization scheme in time which is the classical Runge-Kutta 4th-order method.

Space discretization. Let us set the spatial domain
which contains a topography. On the interval [0, L] in the x-direction, we define , and x i+ 1 ). We now consider the barycenter of the inside cells (see Figure 3.2). Using the diagonals, we split the quadrilateral cell C i,j , for 1 ≤ i ≤ N x and 1 ≤ j ≤ N p , into four different triangles and find the barycenter of each of them:  The barycenter of the quadrilateral cell C i,j is the point of intersection between two lines which pass throughx 1 andx 3 , andx 2 andx 4 , respectively. Thus, the barycenter (x m , y m ) is We also define the centers of the East and West edges to compute the fluxes in Section 3.2 below. For a given cell C ij , let us call the centers of the West edge as (x i− 1 2 , p i− 1 2 ,j ) and the centers of the East edge as (x i+ 1 2 , p i+ 1 2 ,j ) , respectively. Then, the centers of the edges are defined below: (3.6) We introduce the flat control volumes along the boundary of M to impose the boundary conditions: We set the centers of the flat control volumes in (3.7) as follows: We also define the segments Γ i,j+ 1 2 and Γ i+ 1 We now introduce the finite volume space V h : We then write where χ C i,j is the characteristic function on C i,j . For the computation of ω in Section 3.4, and φ in Section 3.5, we construct the quadrilateral cells C i,j+ 1 2 to employ the finite volume derivatives as in [2], [18], and [17]: , and x i, j+1 , Figure 3.3. We aim to compute the gradients satisfying the equations below, for 0 ≤ i ≤ N x and 1 ≤ j ≤ N p : We then obtain a l i+ 1 2 ,j+ 1 2 for l = 1, 2, 3, 4 by fixing one of the four variables; see [2] for more details. We then compute u i+ 1 Now, we can define the non singular matrices M i,j+ 1 2 whose row vectors represent the diagonal of C i+j+ 1 2 : Then, we obtain the gradient ∇ h u h : 3.2. Finite volume scheme: Godunov's scheme. In our simulations, we derive our finite volume method from an upwind finite volume method, see e.g. [24], to implement our schemes. From (2.21) -(2.23) we define the finite volume space for u The finite volume space for ω is For our last unknown φ x , the finite volume space is By integrating (2.22) on each cell to project (2.22) onto the finite volume spaces, we obtain and u, S, and Φ x are as in (2.21). In this subsection, we focus on the fluxes and find u = (T, q, u) using upwind schemes. We then look for ω and φ x separately in Sections 3.4 and 3.5.
Using the divergence theorem, we obtain that where n is the outer normal vector of the cell C i,j , 1 ≤ i ≤ N x and 1 ≤ j ≤ N p . We then rewrite the second term of (3.19) as where the vertical fluxes G i,j+ 1 2 and G i,j− 1 2 are respectively the up and down fluxes, and the horizontal fluxes F i+ 1 2 ,j and F i− 1 2 ,j are respectively the West and East fluxes. We aim to find Let n i+ 1 2 ,j be the vector for the East and West boundaries such that (3.25) n i+ 1 2 ,j = (1, 0). The vertical fluxes are defined as follows: for 1 ≤ i ≤ N x and 0 ≤ j ≤ N p , We note that the barycenter of the trapezoidal cells are well-aligned in p-direction, then we reconstruct u i, , ω i,j+ 1 2 using the interpolation method. For instance, we approximate u i,j− 1 2 and ω i,j− 1 2 , for 1 ≤ i ≤ N x , 1 ≤ j ≤ N p + 1, by (3.28) see e.g. Figure 3.4.
) on the edges where (u i,j−1 , ω i,j−1 ) and (u i,j , ω i,j ) are at the center of each trapezoid cell.
For the horizontal fluxes F i+ 1 2 ,j and F i− 1 2 ,j the normal vectors are (1, 0) due to the proposed spatial discretization. (see Figure 3.1). Hence, for 0 ≤ i ≤ N x and 1 ≤ j ≤ N p , the horizontal fluxes are Figure 3.5 shows how we interpolate u i,j and u i+1,j to obtain u i+1/2,j . The expression of u i+ 1 2 3.3. Computation of the projection methods. In general, the initial condition of u does not follow the compatibility condition (2.11). Hence, the projection method in Section 2.3 plays an important role in our problem. We first set α h λ x such that where x i are the x-coordinates of the barycenter of the cells. We adopt a forward difference scheme for the derivative in x, then (2.29) becomes where andũ is the solution u in Section 3.2. We utilize the mean zero condition in (2.25) to impose a boundary condition for α Nx : From equations (3.33) and (3.35) we obtain the value of α h using an LU decomposition.
Remark 3.1. We consider the Euler method in time, as an example, and write our projection method as follows: We define the intermediate steps where φ n x = (φ x ) h as in (3.42) below, then we find the solutionũ n+1 . We note that u n+1 does not necessarily satisfy the compatibility condition. We therefore apply the projection method to obtain u n+1 ; Then we easily obtain (2.22) 3 by adding (3.36) to (3.37).

Computation of ω.
We look for ω considering the incompressibility in (2.22) 2 and the given value u. We first write a discretized form of (2.22) 2 such that h and ∇ p h denote the discrete directional derivatives in x and p, respectively. Thanks to the proposed spatial discretization, we choose the standard finite difference methods (FDM) for ∇ p h . However, for ∇ x h , it is not accurate to use standard FDM. Instead, we utilize finite volume derivatives on C i,j+ 1 2 which is defined in (3.15). Then, we rewrite (3.38) with ω i,0 = 0 because ω h ∈ W h . We rewrite (3.39) in a matrix form where F h (u h ) is the right-hand side and A h ω h is the left-hand side of (3.39). Since we have ω i,0 = 0, for 1 ≤ i ≤ N x , equation (3.40) has a unique solution ω ∈ W h for a given u h ∈ V h .

3.5.
Computation of φ x . We recall (2.16) to compute (2.22) 3 . We then project φ x onto the space K h in (3.18), and write . We utilize the finite volume derivative on C i,j+ 1 2 to compute T x as in Section 3.4. Then, (2.16) becomes Considering the boundary condition in (3.18), we complete the computation in (3.42).

Numerical simulations
In this section we carry out numerical experiments. We first modify (2.22) by removing S and adding a source terms so that we see the effectiveness of the proposed scheme from Section 3.6. In Section 4.1, we test our scheme and estimate its rate of convergence, numerically. In Section 4.2 and 4.3, we perform physically plausible computations by solving the full equations (2.22) supplemented with the proper boundary conditions in (2.23).

Analytic case.
In this section we use the following system of equations:

(4.4)
Using these analytic functions, we find the rate of convergences for the proposed scheme. In the simulations, we set ∆t = 10 −2 , and [N  Figure 4.1 show the relative L 2 errors at t f = ∆t × k, where k = 100, for different spatial discretizations. Here, we define the relative L 2 error for e.g. T by (4.5) where T EX is the exact solution in (4.2) and T N U M is the numerical solution of (3.23) in (4.1). Also, we denote |C i,j | is the area of the (i, j)-th cell. The relative L 2 errors for u and ω are obtained in the same way. In Table 4.1, we observe the rate of convergence of the numerical solutions T , u and ω. Figure 4.1 shows a first order convergence of the numerical solutions for our scheme, with t f = 100∆t. We observe the same rate of convergence for greater t f .
In our simulation we average T for every 18 time step, i.e. m = 18n. For u, ω, and φ, we average the cells in the same way but we take m = n.

The initial conditions.
We recall that the temperature can be written as T (x, p, t) =T (p) + T (x, p, t) and we take T (x, p, t = 0) = 0; therefore where T 0 = 300K and ∆T = 50. The initial value of the humidity q is where q s is the saturation defined in (2.3). Figure 4.2 shows the shapes of the initial condition for q and the saturation q s at a certain height (around 200m away from the earth). Note that we choose a slightly under saturated initial condition for q to see how the mountains produce saturation and rain, that is, To describe the initial condition for the velocity u, we first introduce the intermediate  valueũ which the velocity before it is projected. We then writẽ u(x, p, t) =ū + u (x, p, t), where (4.9)ū = 7.5m/s, u (t = 0) = 2 cos pπ p 0 cos 2nπx L .
Note thatũ does not follow the compatibility condition in (2.13). Hence, after applying the projection method described in (3.36) and (3.37), we obtain the initial condition for u, that is, u = F h (ũ). Figure 4.3 reports on the striking numerical advantage of using such a method for a simulation of u in good agreement with the natural constraints associated with the problem at hand such as the compatibility condition defined in (2.13) emphasizing the constant horizontal profile that a vertical integration of u must satisfy. As one can observe on Fig. 4.3, when the projection method is applied, the deviations from such a constant horizontal profile are reduced by a factor 10 4 , reducing in other words, the error of the simulated u in satisfying (2.13) by the same factor.
For the boundary condition at x = 0, we use g T (p) =T (p), g q (p) = q s (T , p), g u (p) = F h ū + 2 cos pπ p 0 . (4.11) In the following simulation, we choose the parameters:  p A u dp, before and after the projection. Since u is positive, the flow moves from West to East. On the upstream side of the mountain, the temperature is lower and the humidity is higher, whereas on the downstream side of the mountain, the temperature is higher and the humidity is lower for sufficiently large t. These results are coherent with the physical context. Note that in Figure 4.4 we magnify the value of T and q near the ground to see in detail the behavior of T and q. Figure 4.5 shows behaviors of u (horizontal velocity) and ω (vertical velocity). In Figure 4.6, we present the 1D curve for T and q at t = 15000 along the mountain, i.e. along the dotted line in Figure 4.7. We observe that the figures show asymmetries: it is more humid on the left, and warmer on the right. There are consistent with the fact that it rains on the left side of the mountain, where the wind comes from. Figure 4.8 shows the time-evolution of the (spatial) L 2 -norm of T , q, u, and ω. We observe that the numerical solutions reach a steady state for sufficiently large t (e.g. t > 15000), while the time-evolution the L 2 -norm of T , q, and u exhibit a transient growth that suggests the presence of nonnormal modes which in the present context can be explained as resulting from the topography which breaks the symmetry of the (physical) domain leading typically to non-orthogonal modes associated with the linearized operator. Under these circumstances, disturbances can develop in the system that is favorably configured to undergo rapid transient growth, even in the absence of any growing modes. Because the modes are non-orthogonal, they have a non-zero projection on one another, so it is possible to superpose them to produce disturbances that initially grow rapidly. Such a growth can be further amplified by nonlinear effects or by noise such as documented in the literature; see e.g. [29] and the supporting information of [5] for an analogous situation in a simple model. We turn now to the investigation of such a phenomena in the next section by including some random small-scale disturbances in the model formulation that as we will see cause the appearance of propagating (modulated) waves in response to such a noise-forcing. We refer to [14], [15], [26], [31], [39], and [38] for manifestations of nonnormal modes in various settings borrowed from hydrodynamic or geophysical fluid dynamics.       u σ (x, p, t) for σ = 5, for x ∈ [43.5 km, 75 km] (mainly over which the bombardment occurs) and at a "fixed" value of the topography level. Bottom panel: Same for the temperature profile T σ (x, p, t). Both fields exhibit modulated waves traveling eastward. Here 480 cells are randomly "bombarded" according to (4.14) as time flows, for a spatial resolution of the model corresponding to N x = N p = 100. The figure shows that large-scale recurrent patterns arise from such a small-scale random forcing.

4.3.
Recurrent large-scale patterns from random small-scale forcing. As motivated at the end of the previous subsection, we analyze here from a numerical viewpoint, the effects of a stochastic perturbation to the model formulation. In that respect, we adopted to stochastically perturb only the u-equation in the inviscid primitive equations considered here, which turned out to be enough to illustrate our purpose. More precisely, where η D (x, t) is a random "bombardment" over the region D [45km, 75km] × [0, 2000] (at the east side of the mountain) which takes the following form: (4.14) η(x, t) =  u(·, p, t) L 2 (0,L) (blue curve) and of t → T (·, p, t) L 2 (0,L) (red curve). Bottom panel: Anomalies of t → u(·, p, t) L 2 (0,L) and t → T (·, p, t) L 2 (0,L) , normalized by their respective (empirical) standard deviation.
Here, W t denotes a one-dimensional Brownian motion and the N centers x j (t) of the balls B(x j (t), r j ) are drawn uniformly in D as t flows. In practice the radii r j are chosen to take a fixed value r so that B(x j (t), r) ⊂ D for all t, with r to be characteristic of some small spatial scales for the problem at hand (in what follows r ≈ 1 km). Such a noise term can be argued to model physical processes that are not accounted for by the given equations such as for instance, the vortices that would arise on the east side of the mountain from a sufficiently large initial horizontal component of an eastward wind.
It has been observed that on a spatial resolution of the model corresponding to N x = N p = 100, such a noise term can help trigger interesting dynamics such as illustrated in Figure 4.9. Indeed for σ = 0, the system is a stationary regime (steady state) whereas as σ starts to increase, a complex spatio-temporal dynamics takes place in the u σ -fields (horizontal velocity with (4.13)) as well as the T σ -fields (temperature with (4.13)); see Figure 4.9 below for σ = 5. 1 A closer look at the horizontal velocity and the temperature fields 2 reveals that recurrent large-scale patterns, while evolving irregularly in time 3 , are now the dominant ingredients of the time evolution of these scalar fields. These patterns are mainly expressed here, as waves traveling eastward whose amplitude is irregularly modulated as time flows; the dominant "quasi-period" for the T σ -field being noticeably larger than for the u σ -field as can be observed on Figure 4.9. Without entering in a detailed analysis of the space-time variability of such patterns that could be performed for instance by some multivariate data-adaptive spectral methods [16], a simple unidimensional spectral analysis of the time-evolution of the L 2 -energy (in the x-direction) contained in the respective fields gives already good information about the recurrence characteristics of such fields. 4 In what follows, we will denote by t → u(·, p, t) L 2 (0,L) and t → T (·, p, t) L 2 (0,L) these respective timedependent energies. Figure 4.10 below reports on a standard numerical estimation of the power spectrum 5 associated with the time-variability of these energies, such as obtained from their corresponding autocorrelation functions [13,16]. The latter are estimated from the respective anomalies reported on the bottom panel of Figure  4.10, after normalization by the standard deviation to plot the curves on a similar order of magnitude.
The numerical results indicate that the signal t → T (·, p, t) L 2 (0,L) contains a broadband peak that stands above an exponentially decaying background at lowfrequencies within the band [0, 2.5 × 10 −3 ]; see red curve on top panel of Figure 4.10. Associated with the signal t → u(·, p, t) L 2 (0,L) , a broadband peak stands also above an exponentially decaying background, but with much less energy contained in it and spread over a broader range of frequencies (almost over [0, 5 × 10 −3 ]); see blue curve on top panel of Figure 4.10. 4 In other words, the choice of the L 2 -energy as observable allows here to capture key features of the variability of the given spatio-temporal fields. We mention that such comments have to be understood within the language of the spectral theory of dissipative dynamical systems; see [4] for a brief introduction on the topic. 5 also known as the power spectral density. Interestingly, these patterns develop much less spatial irregularities as time flows compared to those of Figure 4.9. The statistics of t → T (·, p, t) L 2 (0,L) are still non-Gaussian (not shown). The latter property is here again a signature of some nonlinear effects triggered by the noise term (4.14); nonlinear effects which however give rise to a less complex spatial structure of the patterns that develop as time flows, compared to the case with topography.
The fact that the broadband peak associated with the time evolution of u(·, p, t) L 2 (0,L) spreads over the range [2.5×10 −3 , 5×10 −3 ] in the frequency domain, is consistent with the higher-frequency time evolution exhibited by the field u compared to the field T as can be observed on Figure 4.9 in the space-time domain, as well as the higher-frequency oscillations exhibited by the evolution of u(·, p, t) L 2 (0,L) compared to the one of T (·, p, t) L 2 (0,L) , in the time domain alone; see Figure 4.10 bottom panel.
The spatio-temporal evolution of u and T exhibit thus recurrent patterns that are modulated in time and occur on a multiplicity of scales, whose the dominant ones are significantly larger than the spatio-temporal scales on which the random forcing act upon. Such a complex dynamical behavior can be argued to be consistent with the idea that the noise has triggered some nonlinear effects that were not expressed when σ = 0; idea further supported by the non-Gaussian character of the model's dynamics such as observed on the probability density function associated with t → T (·, p, t) L 2 (0,L) ; see Figure 4.11.
Interestingly the characteristics of the large-scale patterns are noticeably different in the case without topography compared to the case with topography, while still resulting from nonlinear effects triggered by the noise; see Figure 4.12. To summarize, the random forcing (4.14) combined with the numerical scheme developed in this study, allow for a nice illustration of the fact that the topography is a determining factor for the generation of complex large-scale patterns in geophysical fluid models.