Long's Equation in Terrain Following Coordinates

Long's equation describes two dimensional stratified atmospheric flow over terrain which is represented by the geometry of the domain. The solutions of this equation over simple topography were investigated analytically and numerically by many authors. In this paper we derive a new terrain following formulation of this equation which incorporates the terrain as part of the differential equation rather than the geometry of the domain. This leads to new analytic insights about the solutions of this equation and enable us to compute steady state gravity wave patterns over complex topography.


Introduction
Long's equation [1,2,3,4] models the flow of stratified incompressible fluid in two dimensions over terrain. When the base state of the flow (that is the unperturbed flow field far upstream) is without shear the numerical solutions (in the form of steady lee waves) of this equation over simple topography (i.e. one hill) were studied by many authors [5][6][7][8][9][10][11][12][13]. The most common approximation in these studies was to set Brunt-Väisälä frequency to a constant or a step function over the computational domain. Moreover the values of two physical parameters which appear in this equation were set to zero. (These parameters control the stratification and dispersive effects of the atmosphere -see Sec 2.) In this (singular) limit the nonlinear terms and one of the leading second order derivatives in the equation drop out and the equation reduces to that of a linear harmonic oscillator over two dimensional domain. Careful studies [8] showed that these approximations set strong limitations on the validity of the derived solutions [9].
Long's equation also provides the theoretical framework for the analysis of experimental data [14,15,29] under the assumption of shearless base flow. (An assumption which, in general, is not supported by the data). An extensive list of references appears in [16,17,18].
An analytic approach to the study of the solutions of this nonlinear equation was initiated recently by the current author [19,20,21]. We showed that for a base flow without shear and under rather mild restrictions the nonlinear terms in the equation can be simplified. We also identified the "slow variable" that controls the nonlinear oscillations in this equation. Using phase averaging approximation we derived for self similar solutions of this equation a formula for the attenuation of the stream function perturbation with height. This result is generically related to the presence of the nonlinear terms in Long's equation. The impact that shear has on the generation and amplitude of gravity waves was investigated by us in [20]. A new representation of this equation in terms of the atmospheric density was derived in [21].
One of the weak aspects of Long's equation is related to the fact that the terrain is represented by the shape of the domain and the boundary conditions. As a result the impact of different terrains on the solution of this equation can only be studied numerically. Furthermore discretization errors which occur in the representation of the terrain render it impractical to consider complex terrain.
In part these errors are due to the scale of the terrain relative to the computational domain.
Accordingly only simple topographies which were represented by one hill were considered in the literature. Furthermore even for these simple topographies only approximate boundary conditions were applied at the terrain. (See discussion in Sec 2).
With this motivation it is our objective in this paper to derive a terrain following formulation of Long's equation in which the terrain is incorporated as part of the coefficients of the differential equation, and the computational domain is always a rectangle. This new representation makes it possible to derive new analytic insights about the solution of this equation in some limiting cases.
It will make it easier also to study how the solution varies as a function of the terrain and other parameters that appear in the equation.
The plan of the paper is as follows: Sec. 2 presents a short review of Long's equation and some aspects of its solutions. In Sec. 3 we derive the new formulation of this equation. Sec 4 considers some analytic and geophysical aspects of this new formulation while Sec 5 compares its numerical solution over three different terrains. These simulations are motivated by recent experiments in the Alps region to educe properties of gravity waves from experimental data [14]. We end up in Sec 6. with a summary and conclusions.

Long's Equation -A Short Overview
In two dimensions (x, z) the flow of a steady inviscid and incompressible stratified fluid is modeled by the following equations: where subscripts indicate differentiation with respect to the indicated variable, u = (u, w) is the fluid velocity, ρ is its density p is the pressure and g is the acceleration of gravity.
We can non-dimensionalize these equations by introducinḡ where L represents a characteristic horizontal length, and U 0 ,ρ 0 represent respectively the free stream velocity and averaged base density (i.e. hereρ 0 is a constant). N 2 0 is an averaged value of the Brunt-Väisälä frequency where ρ 0 = ρ 0 (z) is the base density.
In these new variables eqs (2.1)-(2.4) take the following form (for brevity we drop the bars) β is the Boussinesq parameter [13] (this name has nothing to do with the "Boussinesq approximation") which controls stratification effects (assuming U 0 = 0) and µ is the long wave parameter which controls dispersive effects (or the deviation from the hydrostatic approximation). In the limit µ = 0 the hydrostatic approximation is fully satisfied, [10,11].
In view of eq. (2.7) we can introduce a stream function Ψ so that Using this stream function we can rewrite eq. (2.8) as J{ρ, Ψ} = 0 (2.14) where for any two (smooth) functions f, g Eq. (2.14) implies that the functions ρ, Ψ are dependent on each other and we can express each of them in terms of the other. Thus we can write Ψ as Ψ(ρ) (or ρ as ρ(Ψ) [21]).
After a long algebra one can derive the following equation for Ψ [22,1,13] is the nondimensional Brunt-Väisälä frequency. We observe that in this definition N 2 is a function For example if we let i.e consider a shearless base flow with lim and eq. (2.16) becomes: It is evident from this derivation that different profiles for the base flow as x → −∞ will lead to different forms of S(Ψ) [20].
For a general base flow in an unbounded domain over topography with shape f (x) and maximum height H 0 the following boundary conditions are imposed on Ψ where the constant in eq.(2.22) is (usually) set to zero. As to the boundary condition at x → ∞ it is appropriate to set lim x→∞ Ψ(x, z) = Ψ 0 (z) (in spite of the fact that Long's equation contains no dissipation terms). However over finite computational domain only radiation boundary conditions can be imposed in this limit. Similarly as z → ∞ it is customary to impose (following [7]) radiation boundary conditions. (The imposition of these boundary conditions is discussed in detail in Sec 4.1).
For the perturbation from the shearless base flow We observe that when |τ | ≪ 1 the boundary condition (2.22) can be approximated by When N is constant eq. (2.24) is invariant with respect to translations in x, z and hence admits self-similar solutions of the form η = f (kx + mz) [19]. These solutions are interpreted as gravity waves that are generated by the flow over the topography.
From a numerical point of view it is a common practice [7,8,13] to solve eq. (2.24) in the limit β = 0 and µ = 0 with constant N over the domain. However observe that the definition of N in Long equation is given by (2.17) and it depends on Ψ. In some other numerical simulations the computational domain is divided into subdomains where N is constant in each subdomain but this led to numerical instabilities at the interface between these subdomains.
In these limits Eq. (2.24) reduces then to a linear equation We observe that the limit β = 0 can be obtained either by letting U 0 → 0 or N 0 → 0. In the following we assume that this limit is obtained as U 0 → 0 (so that stratification persists in this limit and the leading term in N 0 is not zero).
Eq. (2.27) is a singular limit of Long's equation as one of the leading second order derivatives drops when µ = 0 and the nonlinear terms drops out when β = 0 and N is constant. This approximation and its limitations were considered numerically and analytically [6,7,19,20] and was found to be justified only under strong restrictions even under the assumption that the base flow is shearless. Nevertheless it is used routinely in the actual analysis of atmospheric data [14,15,16].
The general solution of eq. (2.27) is where the functions p(x), q(x) have to be determined so that the the boundary conditions derived from eq. (2.22),(2.26) and the radiation boundary conditions are satisfied. These lead in general to an integral equation for p(x) and q(x) and it easy to show [13] that is the Hilbert transform of q(x). The boundary condition on the terrain becomes; This integral equation has to be solved numerically [6,7,13].
To derive a terrain following formulation of Long's equation which incorporates the terrain in the coefficients of the differential equation (rather than the shape of the domain) we introduce Gal-Chen transformation. If the height of the (bottom) terrain is described by a sufficiently smooth where H is a constant, then this transformation is given bȳ Under this transformation we have Furthermore the expression of the Laplace operator becomes Under this transformation the continuity equation (2.7) becomes then it is a simple algebra to show that eq. (3.5) can be rewritten as From this equation we see that we can introduce a "terrain following stream function" ψ so that Multiplying eq. (2.8) by √ G we can rewrite this equation in the following form: Using eq. (3.8) this can be rewritten asJ whereJ is defined as in eq. (2.15) but with differentiations with respect to (x,z). Equation (3.10) implies that ρ(x,z) = ρ(ψ(x,z)) (and vice versa).
To eliminate the pressure term from eqs. (2.9) and (2.10) we differentiate (2.9) by z and apply the operator µ 2 ∂ ∂x to (2.10) and subtract. We obtain Using eq. (2.8) the first two terms in this equation can be written as Using eqs (3.6) and (3.8) to re-express u 2 + µ 2 w 2 we have The third and the fourth terms in eq. (3.11) can be rewritten using (2.7) as where χ = µ 2 w x − u z is the vorticity. Expressing χ in terms of ψ we have is the "terrain following Laplace operator".
Finally for the right hand side of eq. (3.11) we have Combining all the results contained in eqs. (3.12)-(3.17) we can re-express eq. (3.11) in the following form: where N 2 (ψ) is defined as in eq. (2.17). Hence it follows that, This is the terrain following form of Long's equation. At this juncture it might be asked why one can not "save" this derivation and apply the terrain following transformation (3.1) directly to (2.16). Doing so will yield an extremely complicated equation. This has been avoided in our derivation by the use of the "terrain following stream function" in (3.8).
To determine the function S(ψ) in eq. It follows then that and Long's equation becomes: where n is the normal to the topography which is described by the curve h(x). Hence this normal is given by n = (−h ′ (x), 1). Using (3.6),(3.8) this leads to the boundary condition ψ(x, 0) = constant (3.22) and this constant can be chosen to be zero. The other boundary condition that has to be imposed on ψ is a radiation boundary condition atz = H (which implies that the outgoing wave is not reflected by the boundary).
To obtain an equation for the perturbation from the base state we set ψ(x,z) =z + η(x,z).

Analytic solutions of Long's equation
In the traditional representation of Long's equation the topography determines the shape of the flow domain and as a result it is not feasible to obtain analytic solutions to this equation even in some limits of the parameters β and µ . We now show that this problem can be overcome in some limiting cases when the terrain following formulation of this equation is used.
For brevity we drop in the following the bars over x, z.

The Limiting
Case β = 0, µ = 0 In this case eq. (3.21) simplifies to whose general solution is Here ν = N √ G and A(x), B(x) are functions which have to be determined from the boundary conditions .

The boundary condition (3.22) implies A(x) = −h(x).
To determine B(x) we must apply the radiation boundary condition as z → ∞ on the solution. To this end we must insure that the vertical group velocity of the wave is positive. Using the dispersion relation for hydrostatic flow given in [16,p.181] this group velocity is: where k is the horizontal wave number. We deduce then that the vertical group velocity is positive when kν ≥ 0.
To impose this condition on the solution (4.2) we express A(x), B(x) in Fourier integral form where k is the horizontal wave number. We deduce then that the solution (4.2) can be written as To satisfy the radiation boundary condition for z → ∞ the first and second integral must vanish for k < 0 and k > 0 respectively. Therefore a(k) and b(k) must satisfy We compare now these analytic results with the solution methodology that has been used previously in the literature as was discussed Sec 2. First we note that this analytic solution requires only the direct (and simple) computation of the Hilbert transform of the terrain function h(x) . This is a straightforward procedure even if it has to be done numerically. On the other hand to compute q(x) using (2.29) requires in general the solution of an integral equation. To do so one must use an iterative algorithm which might turn out to be unstable or non-convergent over complex terrain. Furthermore there is the issue of applying the boundary conditions on ψ at the terrain. To this end the procedure discussed in Sec 2 requires the use of the approximations that lead to (2.26). As a result the equation that is used to compute q(x) (eq. (2.29)) is also an approximate equation which will yield at best approximate solution for this function. On the other hand the application of the boundary conditions using the procedure discussed in this section is exact and does not place constraints on the height of the terrain.
From an experimental geophysical point of view it has been a common practice to assume that the gravity wave generated by a flow over terrain is of the form sin(kx + mz) (or similar) [14,15,28].
This has led to difficulties in the eduction of this wave from experimental data. Our results show that this form of the wave is incorrect (at least in principle). Furthermore the numerical simulations that we carry in the next section demonstrate that complex terrain can alter drastically the shape and amplitude of this wave due to interference effects.

The Limiting
Case β = 0, µ = 0 In this limiting case eq. (3.21) becomes Since this is a nonlinear equation we can find an approximate analytical solution using first order perturbation expansion under the assumption that α 2 = N 2 β 2 ≪ 1 (which is satisfied in most practical situations). Expressing ψ approximately as ψ = ψ 0 + α 2 ψ 1 and substituting this expression in (4.7) we obtain to order zero and one in the parameter α 2 the following equations The boundary conditions on ψ 0 , ψ 1 are given by eq. (3.22) at z = 0 and radiation boundary conditions as z → ∞.
Solving these (linear) equations for ψ 0 and ψ 1 we obtain the following expressions for their (4.10) The determination of the functions A(x), B(x), C 1 (x), C 2 (x) from the boundary conditions can be done using the same procedure outlined in the previous subsection. However as it is algebraically cumbersome we omit the details.

Numerical Simulations over Complex Terrain.
In the previous section we presented analytical solutions to Long's equation in some limiting cases. In general one has to resort to numerical simulations of Long's equation. However as it was said in the introduction the terrain following form of Long's equation obviate the need to discretize the terrain in these numerical procedures and therefore lead to better representations of complex terrains. Also the application of the boundary condition at the terrain is being simplified considerably and is treated exactly contrary to the approach discussed in Sec 2 (see eqs. (2.26) and (2.29)).
These simulations are carried in order to validate our equations and explore the different flow patterns (and possible special effects) that can be predicted over complex topography. The geophysical motivation for these simulations is related to recent experiments to measure and educe gravity waves and their properties over the Alps from balloon data. This endeavor faced several difficulties. In part these difficulties can be traced to the fact that the terrain over which the measurements were made contained (at least) two summits rather than one [15].
In this section we carry simulations over three terrains. These consist of one, two and three hills with the following shape functions In all cases we solved (3.24) for the perturbation with the following parameters; Here N represents the nondimensional Brunt-Väisälä frequency which was defined in (2.17).
To solve for the perturbation η over a finite two dimensional domain which after the trans- To solve (3.24) under these settings we used Matlab [28]. Central finite differences approximations were used to discretize Long's equation on the domain grid and a fixed point iterative algorithm was implemented to solve the resulting equations. The convergence criteria for the iterations was that the step error |η m+1 − η m | was less than 1.10 −9 where m is the iteration number.
Convergence was achieved in less than 100 iterations. --------------------It should be noted however that η in these simulations represents the perturbation from the "terrain following stream function" which was defined in (3.8). Using (3.6), (3.8) we find that the perturbation to the base flow field u p , w p is given by wherez is defined in (3.1). From this vector field we can compute the "regular" stream function where ν(x, z) satisfies A zoom-plot of φ(x, z) (on part of the computational domain) on the region around the three topographies discussed above is presented in F igs 1, 2, 3. We observe that the phase lines of the stream function in all these figures tilts with height. In all three figures the gravity waves are in the lee of the topography.
These figures demonstrate the impact that complex terrain can have on the structure of gravity waves.
F ig 1 displays the results of the simulation for one hill (described by (5.1)). The plot corresponds to the classical results that appeared in the literature for this case. The results for two hills in F ig 2 show some wave activity between the two hills and less organized waves with somewhat smaller amplitude in the lee of the two hills. This might be due to interference effects between the waves generated by the two hills. For three hills the simulation shows the same wave activity between the two hills (on the left) as in F ig 2 but then a "quiet" zone between the second and third hills in which φ is almost constant. The waves in the lee of the three hill appear regular as in F ig 1.
6 Summary and Conclusions.
We derived in this paper a terrain following formulation of Long's equation in which the topography is "absorbed" in the coefficients of the differential equation representing the flow rather than being part of the boundary conditions. We used this representation to solve Long's equation analytically in some limiting cases and numerically over complex topography. The new formulation also opens the possibility to develop analytical estimates which compare the solutions of this equation over different topographies .
From a geophysical point of view it well known that some present models for the generation of gravity waves over estimate this effect [25,27]. Partially, this is due to the fact that they use oversimplified representation of the terrain. Furthermore they do not take into account the effects that are due to complex terrain (as demonstrated by our simulations). We believe that the new form of Long's equation will make it easier to consider more realistic representations of the terrain and its effect on the generation and propagation of gravity waves.