Reflection and transmission of gravity waves at non-uniform stratification layers
Christopher Pütz, Mark Schlutow, Rupert Klein, Vera Bense, Peter Spichtinger
RReflection and transmission of gravity waves at non-uniform stratification layers
Christopher P¨utz , Mark Schlutow , Rupert Klein , Vera Bense ,and Peter Spichtinger Department of Mathematics, Freie Universit¨at Berlin, Germany Institute for Atmospheric Physics, Johannes GutenbergUniversit¨at Mainz, Germany
Abstract
The present study focuses on the interaction of gravity waves in the atmosphere with the tropo-pause. As the vertical extent of the latter is small compared to the density scale height, wavepropagation is described by the Taylor-Goldstein equation as derived from the linearised Boussi-nesq approximation. Of particular interest in the construction of gravity wave parameterisationsfor the upper atmosphere are the transmission and reflection properties of the tropopause asthese determine the upward fluxes of energy and momentum carried by internal waves.A method is presented that decomposes internal waves explicitly into upward and downwardpropagating contributions, thus giving direct access to transmission and reflection coefficientsof finite regions of non-uniform stratification in a stationary atmosphere. The scheme utilizes apiecewise constant approximation for the background stratification and matches up- and down-ward propagating plane wave solutions in each layer through physically meaningful couplingconditions. As a result, transmission and reflection coefficients follow immediately.In the limit of an increasing number of layers the method leads to a reformulation of theTaylor-Goldstein equation in a particular set of variables. Numerical integration of this non-constant coefficient differential equation provides a representation of Taylor-Goldstein solutionsthat also distinguishes explicitly between the upward and downward travelling wave branchesof the dispersion relation and hence gives access to transmission coefficients also for smoothlystratified layers.The multi-layer solutions are not only shown to converge to the limit solution quadraticallywith the number of layers, but are also found to be surprisingly accurate – and hence efficient –for very small numbers of vertical layers. The results obtained for some test cases are in goodagreement with several existing results as well as with two-dimensional numerical solutions ofthe full non-linear pseudo-incompressible equations for a vertical slice. Yet, by revealing the up-and downward travelling wave components explicitly, the multi-layer solutions offer alternativeinsights into the interaction of gravity waves propagating through non-uniform stratification. a r X i v : . [ phy s i c s . a o - ph ] D ec he present paper focuses on internal wave eigenmodes of a stratified atmosphere. Yet, italso serves as the basis for the development of a new numerical method for the propagation ofnon-stationary wave packets described in a companion paper. Keywords: gravity waves non-uniform stratification tropopause Taylor-Goldstein equationnumerical simulation Introduction
Gravity waves arise from fluid displacements in a vertically stably stratified medium,in which the buoyancy acts as restoring force. This includes, to a large part, also theearth’s atmosphere. The propagation of atmospheric gravity waves has been the subjectof a number of earlier studies. One of the first investigations goes back to Scorer (1949),who focused on orographically generated waves and how they can be trapped in the leeof a mountain ridge.An important characteristic of gravity waves is the ability to transport energy hori-zontally as well as vertically. Eliassen and Palm (1961) analysed waves carrying energyupward and downward in the context of orographically excited gravity waves that reflectfrom vertically varying stratification and wind. They used a piecewise-constant approx-imation for both stratification and wind and found local solutions. These were matchedat the discontinuities of the approximation. The calculations were carried out for a two-and a three-layer atmosphere. The multi-layer method to be introduced below is alsobased on these concepts.By a similar approach Danielsen and Bleck (1970) examined mountain waves andapproximated key atmospheric parameters by piecewise exponential functions, which al-lowed them to solve the governing equations by combinations of Bessel functions. Suther-land and Yewchuck (2004) attended to the topic again and scrutinised the phenomenonof wave tunnelling, which describes the energy transport over a finite layer of decreased,or even vanishing, stratification. They undertook a mathematical analysis as well as lab-oratory experiments to support their findings. Brown and Sutherland (2007) expandedthe theory by allowing for shear flow over an unstratified layer. Both scenarios were laterexamined numerically by Nault and Sutherland (2007) who provided numerical solutionsfor plane wave transmission in arbitrary stratification and wind. As a concrete example,they performed simulations for an atmospheric stratification and wind profile that wasobserved over Jan Mayen island.Diffraction through a slit and back-reflection from a slope is covered by Mercier et al.(2008), mainly through experimental work. To keep track of the direction of the wavepropagation, they use a demodulation of the measured wave signal. B¨uhler (2009) uses amethod similar to that introduced by Eliassen and Palm (1961) for the two-dimensionalshallow water equations with water depth that changes only with one of the spatialcoordinates. It is approximated by a piecewise-constant depth which allows for explicitsolutions in each of those layers.The main focus of the present work is on the transmission and reflection of gravitywaves at confined regions of non-uniform stratification. The governing equations are thelinearised Boussinesq equations, which can be combined into a single equation for one ofthe dynamic variables. Since this “Taylor-Goldstein equation” does not allow for explicitsolutions in general, we use a piecewise-constant approximation of the stratificationwith finitely many layers to generate an analytically accessible approximate solutionfrom which we can compute the transmission and reflection of gravity waves. In eachlayer of constant stratification an explicit solution is constructed as a superposition ofupward- and downward-travelling plane waves. These are matched at the discontinuities3f the piecewise-constant approximation of the background rendering the perturbationsof vertical wind and pressure continuous across these interfaces. Conservation of energyis used to derive a transmission coefficient.By a reformulation, we find a new set of differential equations that describes the limitof an infinite number of discrete layers. These limit equations are solved numerically. Theresulting transmission coefficients are then compared with those calculated for a finitenumber of layers. For an increasing number of layers, the solution is found to converge tothe limit solutions with second order accuracy. Yet, even for a moderate number of layersthe results turn out to be rather accurate, and this allows for very efficient computationalestimates of the transmission and reflection properties of tropopause-like layers in theatmosphere. Another advantage of the presented method is its ability to keep trackof upward and downward propagating waves locally in arbitrary stratification profiles.Moreover, an extension of this method to a solver for non-stationary wave packets hasbeen achieved, and this is the topic of the companion paper (Puetz and Klein , 2018).The present paper is structured as follows. In section 2 we present the method de-veloped to analytically compute a transmission coefficient. Numerical computations arepursued in section 3 to provide evidence of the correctness of the method as well as anerror analysis. We show results for selected stratification profiles in section 4. Section 5summarizes results from numerical solutions of the full Boussinesq equations to furtherback up our theoretical findings. In section 6, we discuss further extensions as well aslimitations of the method. In particular, the adaptation to shear layers and wave packetsis of special interest for further studies.
We are interested in wave propagation through the tropopause. The latter is charac-terised by sharp changes in the stratification and also by strong jet winds. Above andbelow the tropopause, we have almost uniform stratification. The extent of the tropo-pause is small compared to the density scale height. Therefore, as we only focus on thepropagation through this rather shallow area, we can apply the Boussinesq approxima-tion.The starting points for our investigation are the two-dimensional inviscid Boussinesqequations with small amplitudes and an atmosphere at rest. If we linearise these equa-tions, they can be written as a single equation for one of the dynamic variables (see forexample Sutherland (2010)), here done for the vertical velocity w : (cid:18) ∂ ∂x + ∂ ∂z (cid:19) ∂ w∂t + N ∂ w∂x = 0 , (1)where N is the Brunt-V¨ais¨al¨a frequency, t denotes time and x and z are the horizontaland vertical coordinates. As long as N does not depend on x and t , the equation admitshorizontally and temporally periodic solutions, i.e., w ( x, z, t ) = ˆ w ( z ) exp( i ( kx − ωt )) (2)4y convention, we consider only ω > k >
0, i.e., the phase velocity pointsin the positive x -direction. This is no restriction, since the case k < k >
0. The partial differential equation (1) then can be transformed intothe ordinary differential equationd ˆ w d z + k (cid:18) N ω − (cid:19) ˆ w = 0 . (3)This actually corresponds to a Fourier transformation of equation (1) in the horizontaland time coordinates and the resulting equation is widely known as the Taylor-Goldsteinequation, which cannot be solved explicitly for non-constant Brunt-V¨ais¨al¨a frequency N .But since we are only interested in a confined region within which N is varying, we canapproximate N in this region by a piecewise constant function. To be more precise, weare given a function N ( z ) = N b , z < z b N c ( z ) , z b ≤ z ≤ z t N t , z > z t , (4)where z b , z t are the bottom and top of the tropopause, respectively (or any region ofinterest in general), N b , N t are constant values of N in the bottom and top layer, respec-tively, and N c is a (continuous) function of z with N c ( z b ) = N b and N c ( z t ) = N t . Nextwe partition the continuous part into a piecewise-constant function. This allows us tofind explicit solutions in each layer. This method is neither new nor ground-breaking,but the way we use it to compute wave transmission strikes a new path. Moreover, aswe will see later, we will gain structural insights into the solution and we are able to givenumerical evidence that justifies the ansatz, despite its simplicity.Let J be a positive integer. For now, J is fixed. Define an equidistant grid of J pointsfrom z b to z t : z j = z b + j − J − z t − z b ) for j = 1 , . . . , J (5)and set N := N ( z ) , N J +1 := N ( z J ) and N j := N ( z j ) + N ( z j − )2 for j = 2 , ..., J. (6)This can be understood as a piecewise function˜ N ( z ) = N j , z ∈ I j , (7)where I j = [ z j − , z j ) for j = 2 , . . . , J , I the troposphere region z < z b and I J +1 thestratosphere region z ≥ z t . In each single level, we are able to state the Taylor-Goldsteinequation, but now N takes a constant value. In particular, for the level I j , we have theequation d w j d z + k (cid:32) N j ω − (cid:33) w j = 0 . (8)5ach layer admits explicit plane wave solutions of the form w j ( z ) = A j exp( im j z ) + B j exp( − im j z ) , (9)where m j = − k (cid:115) N j ω − A j , B j are the amplitudes of the upward and downwardpropagating wave, respectively. Equation (10) is basically the transformed Boussinesqinternal gravity wave dispersion relation. The representation we use for the wave in equa-tion (9) corresponds to the Hilbert transform method for internal gravity wave analysisby Mercier et al. (2008). Although they were the first to apply this technique in theanalysis of observational data, it goes all the way back to the classical works of Eliassenand Palm (1961) and Booker and Bretherton (1967), who use this representation of planewaves in their theoretical setups. From a mathematical point of view, this demodulationis a rather natural approach to transform a real valued signal into a complex-valued one.The original real signal then corresponds to the real part of the complex one, in thiscase the right-hand side of equation (9). This representation has the advantage of track-ing upward- and downward-propagating waves. In Mercier et al. (2008), however, it isnecessary that the stratification is uniform or at most slowly changing. With the multi-layer method we are going to introduce, we will be able to keep track of upward- anddownward-propagating waves also in non-uniform background that can change stronglyover a fraction of the wavelength. To the best of our knowledge, this is the first tech-nique that is capable of providing this information by construction. In the upcomingcomputation steps, we will work with the Hilbert representation. To obtain the physicalsolution, we can always take the real part of what we computed. It is important to notethat the wave amplitude of the real solution is not computed by taking the real part ofthe complex amplitude, but also using the imaginary part. That said, a non-vanishingimaginary part of the amplitude basically acts like a phase shift of a cosine.Going back to the multi-layer approach, the solution indexed by 1 references thesolution below the tropopause while the solution indexed by J + 1 corresponds to thesolution above it. In particular, the amplitude A belongs to the incident wave, while A J +1 is the amplitude of the transmitted wave and B the one of the reflected wave.Since we are assuming a uniform stratification above the tropopause, there shall be noreflection from upper layers. Hence there is no wave hitting the tropopause from above,i.e., B J +1 = 0. This is a radiation condition which we are going to use later.We need to clarify the way we use the indices on the variables that depend on thelevels I j . All those variables depend implicitly on the (fixed) number J , i.e. for J (cid:54) = J ,for example N ( J ) j (cid:54) = N ( J ) j , where the superscript now reflects the dependence on thenumber of levels. To be precise and keep the variables comparable, one could indexthem by jJ or superscripting them with the number of steps. But since this is not onlycumbersome in writing and reading but also does not provide further benefit (most ofthe time, we are interested in the variables indexed with 1 and J + 1), we omit thisdependence but keep it in mind. 6o obtain a solution over the whole domain, we have to match the local solutionsat the interfaces in a proper way. Physically meaningful conditions require that thevertical wind speed and the pressure are continuous across the interface (see also Drazinand Reid (1981)). By using the polarisation relations in a Boussinesq fluid (see, e.g.,Achatz et al. (2010) for details), the conditions are equivalent to the requirement that w j and w (cid:48) j = d w j d z are continuous at the interfaces, i.e., w j ( z j ) = w j +1 ( z j ) (11a) w (cid:48) j ( z j ) = w (cid:48) j +1 ( z j ) (11b)A single pair of the form (11) gives us two equations for the four unknowns A j , B j , A j +1 , B j +1 . Hence we are able to derive a recurrence relation (cid:18) A j +1 B j +1 (cid:19) = M j (cid:18) A j B j (cid:19) , (12)where the matrix M j is of the form M j = (cid:18) c j d j d ∗ j c ∗ j (cid:19) . (13)The particular matrix entries are given by c j = 12 (cid:18) m j m j +1 + 1 (cid:19) exp( i ( m j − m j +1 ) z j ) (14a) d j = − (cid:18) m j m j +1 − (cid:19) exp( − i ( m j + m j +1 ) z j ) . (14b)For later reference, we introduce the ∗ -operation which changes the sign of the argumentof the exp-function, i.e., c ∗ j = 12 (cid:18) m j m j +1 + 1 (cid:19) exp( − i ( m j − m j +1 ) z j ) (15a) d ∗ j = − (cid:18) m j m j +1 − (cid:19) exp( i ( m j + m j +1 ) z j ) . (15b)As long as m j and m j +1 are real-valued, this corresponds to complex conjugation. Imag-inary values for m occur only when the waves are encountering a region of decreasedstratification, where N < ω . We will see later that these cases are harder to deal withwhen trying to find a limit for an increasing number of layers, hence they have to betreated very carefully.We can state a relation like equation (12) for all j = 1 , . . . , J and combine them toobtain a chain of equations: (cid:18) A J +1 B J +1 (cid:19) = M J (cid:18) A J B J (cid:19) = M J M J − (cid:18) A J − B J − (cid:19) = . . . = (cid:89) k = J M k (cid:124) (cid:123)(cid:122) (cid:125) =: M (cid:18) A B (cid:19) (16)7e have to be careful about the order of the matrix multiplication, since it is in generalnot commutative.Now we need to recall what the different amplitudes with their respective indicesrepresent. A J +1 is the transmitted wave, A is the incident wave and B J +1 = 0, i.e., thereis no downward propagating wave in the uppermost layer. To compute a transmissioncoefficient, we have to relate A and A J +1 . We have A J +1 = M , A + M , B (17)0 = M , A + M , B , (18)where M k,l are the entries of M . Solving the equation system for A and A J +1 showsthat A J +1 A = (cid:18) M , − M , M , M , (cid:19) = det( M ) M , . (19)To compute a meaningful transmission coefficient, we need to find a quantity that isconserved over the whole domain. Since we did not allow for dissipation or backgroundhorizontal wind in the equations, wave energy (sometimes called perturbation energy) isconserved. It consists of kinetic and potential energy. Moreover, we are in a horizontallyperiodic domain, so it is convenient to have a look at the horizontally averaged energy (cid:104) E (cid:105) = (cid:104) E kin (cid:105) + (cid:104) E pot (cid:105) = 12 ρ b N ω | A w | . (20)Here ρ b is the background density and A w is the amplitude of the vertical wind. The unitof (cid:104) E (cid:105) is energy per unit volume, therefore the correct term would be energy density.Energy can be derived from this expression by integrating over a fixed control volume.But since there is no danger of confusion, we stick to the term ”energy” for (cid:104) E (cid:105) . Theconservation equations for kinetic and potential energy can be derived directly from theBoussinesq equations. Horizontal averaging and adding the equations yield ∂ (cid:104) E (cid:105) ∂t + ∂ (cid:104)F z (cid:105) ∂z = 0 , (21)where (cid:104)F z (cid:105) = c g z (cid:104) E (cid:105) is the vertical wave energy flux and c g z = ∂ω∂m denotes the verticalcomponent of the group velocity. A full derivation of the above equations can be foundin chapter 3 . ∂ (cid:104)F z (cid:105) ∂z = 0 (22)basically says that the vertical mean wave energy flux is constant. Again, we use thework of Sutherland (2010) to find out that this flux for a Boussinesq wave in uniformstratification is given by (cid:104)F z (cid:105) = 12 ρ b N ω k sin( α ) cos ( α ) | A w | (23)8here α = arctan( mk ) is the angle between the wave vector and the horizontal. Usingthe identities sin(arctan( x )) = x √ x (24)cos(arctan( x )) = 1 √ x (25)and the internal gravity wave dispersion relation ω = N k √ k + m (26)we obtain that (cid:104)F z (cid:105) = ρ b mω k | A w | . (27)Since we have a region with uniform stratification below and above the tropopause, wecan compare the upward energy fluxes in both of those regions. Recall that ω and k arechosen to be constant. Moreover, we made the assumption that the density does notvary too much over the tropopause, so that we take a reference value (cid:37) for both regions.The transmission coefficient is then defined as the ratio of the upward energy flux aboveand below the tropopause T C := (cid:104)F z (cid:105) above,up (cid:104)F z (cid:105) below,up = m J +1 m (cid:12)(cid:12)(cid:12)(cid:12) A J +1 A (cid:12)(cid:12)(cid:12)(cid:12) = m J +1 m (cid:12)(cid:12)(cid:12)(cid:12) det( M ) M (2 , (cid:12)(cid:12)(cid:12)(cid:12) . (28)In a similar fashion we can define a reflection coefficient, which compares the upwardflux with the downward flux below the tropopause: RC := (cid:12)(cid:12)(cid:12)(cid:12) B A (cid:12)(cid:12)(cid:12)(cid:12) (29)By conservation of vertical energy flux, given by equation (22), we have that T C + RC = 1 . (30) This section is devoted to the investigation of an increasing number of layers, eventuallytending to infinity. First, we have a view on how the transmission coefficients for agiven parameter set change when the number of layers is altered. Then, we will be ableto derive an expression for the number of layers tending to infinity, which results in areformulation of equation (3) in a set of variables that allows the distinction betweenup- and downward propagating wave modes. This immediately gives rise to transmissionand reflection coefficients for an infinite number of layers. We also investigate how quickthe multi-layer method converges to this limit.9t should be mentioned that the multi-layer method itself can be used to construct anapproximate solution to equation (3). When reformulating equation (3) to a first-ordersystem, where one of the variables corresponds to the first derivative of the solution, themulti-layer method with the corresponding matching conditions (11) can be written as aone-step method, similar what is done in Lara (2004), where also a proof of convergenceis given. This puts the multi-layer method on mathematically solid ground.
Our first step towards showing that the method is convergent with increasing numberof layers is to have a look at how the transmission coefficient is influenced by this verynumber, again denoted by J . As a reference setup, we take a stratification profile thatincreases linearly from a value N b to a value N t > N b over a confined region [ z b , z t ] withlength ∆ z := z b − z t : N ( z ) = N b , z < z b N b + z − z b z t − z b ( N t − N b ) , z b ≤ z ≤ z t N t , z t < z. (31)A visualisation of the profile can be seen in the right panel of figure 4. Our computa-tion domain in this investigation covers a frequency range from 0 to nearly N b and awavelength range from 100 ∆ z to about ∆ z . The results for various values of J can beseen in figure 1. The images show a colour plot of the transmission coefficient, where thex-axis corresponds to frequency and y-axis to horizontal wavelength. Here, we are notinterested in the interpretation of the images in themselves, but rather in the transitionfrom small to large J . As one can see, except for the case J = 4 (upper left panel),one can barely spot any difference in the pictures. This leads to the hypothesis that themethod converges to a certain limit for increasing J , and from what we see qualitatively,this convergence seems to be pretty fast. The next step will be to find the limit. Unfortunately, finding the limit has turned out to be much harder than anticipated.Since the matrices M j have complex entries, there seems to be no closed formula forlim J →∞ (cid:89) k = J M ( J ) k , (32)because most known formulas blow up because of the non-vanishing imaginary part(under certain conditions, we managed to find a formula, but in general, these conditionscannot be fulfilled). Hence, a new approach was needed.The idea is to reformulate the limit process as a differential equation for a vectorconsisting of the amplitudes for the upward and downward propagating wave. We knowthat the depth of each layer is h = ∆ z J , so the limit process J → ∞ can also be seen as10igure 1: The panels show the transmission coefficient for profile (31) for given horizontalwavelength λ x and frequency ω for different number of layers J : J = 4 ontop-left, J = 16 on top-right, J = 64 on bottom-left and J = 256 on bottom-right. Except for J = 4, the pictures are nearly indistinguishable. As we willsee in figure 3, we have a convergence rate of about 2.11 z J → h →
0. Moreover, the j -indexed variables A j , B j , m j are approximations oftheir continuous counterparts at z j . In fact, they are approximations at z j + h , whichgoes to z j in the limit h →
0. Also, a reformulation such that e.g. m j = m ( z j ) gives thesame result. By using the recurrence relation (12) we can write (cid:18) A j +1 B j +1 (cid:19) − (cid:18) A j B j (cid:19) = (cid:16) M ( J ) j − I (cid:17) (cid:18) A j B j (cid:19) , (33)where I is the 2-by-2-identity matrix. Dividing now by h and taking the limit h → z -derivative of the amplitudes. Using the short-handnotation A = A ( z ) = (cid:18) A ( z ) B ( z ) (cid:19) for the vector of amplitudes, we haved A d z = lim h → (cid:16) M ( J ) j − I (cid:17) h A . (34)This is now a differential equation for the amplitudes in A . If we want to have anychance of solving it (either analytically or numerically), we have to execute the limitprocess lim h → (cid:16) M ( J ) j − I (cid:17) h . (35)This is done component-wise. The upper-left entry of the matrix inside the limit in (35)is c j − h . If we let h →
0, we get f ( z ) := lim h → c j − h = − m (cid:48) ( z )2 m ( z ) − im (cid:48) ( z ) z. (36)With a similar computation for the off-diagonal entry d j h , we obtain g ( z ) := lim h → d j h = m (cid:48) ( z )2 m ( z ) exp( − im ( z ) z ) . (37)The full derivation for the limit in equations (36) and (37) can be found in the appendix.The respective limits for the starred entries yield the same except for a replacement of i by − i . Hence the differential equation for the amplitudes isd A d z = (cid:18) f ( z ) g ( z ) g ∗ ( z ) f ∗ ( z ) (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) =: M ( z ) A . (38)This equation can only be solved analytically if M ( z ) M ( z ) = M ( z ) M ( z ) holds forall z , z in the integration domain. Unfortunately, this is in general not true for ar-bitrary stratification profiles. What we do know, however, is that if the functions areanalytic over the interval [ z b , z t ], equation (38) has a unique analytic solution for arbi-trary initial data A ( z ) = A , z ∈ [ z b , z t ]. See, for example, Teschl (2012) for the theory12n complex ODEs. It is easy to check that this is the case for stratification profiles andwave parameters such that there is no reflection layer, i.e. a point z r where N ( z r ) = ω .Although securing the existence of solutions, finding analytic or even explicit expres-sions for them will be a nearly hopeless undertaking. Another approach to finding atleast approximate solutions are power series methods. Since analytic functions on anopen subset coincide locally with a convergent power series (see , e.g., Stalker (1998) fordetails), we can make a power series ansatz for the solution of equation (38). In orderto do so, we extend [ z b , z t ] to an open subset of the complex numbers, in which M isstill analytic. But although the matrix has no singularities for real values, it has somefor certain complex numbers, which drastically restricts the radius of convergence of thepower series solution to a value that is not guaranteed to be large enough to cover thewhole region of interest. However, it would be possible to partition the interval intosmaller segments and finding the power series solution in each segment, but this proce-dure is very tedious and still only yields a solution up to a certain precision. We will seein the upcoming error analyses that the multi-layer method yields very accurate resultsnotwithstanding that it is a much easier-to-apply technique. Hence, the evaluation ofequation (38) will be done numerically.This leaves us with another challenge. If the stratification profile contains a reflectionlayer z r such that N ( z r ) = ω , the entries of the coefficient matrix M tend to infinity,because m ( z r ) = 0. By regarding equation (38) as a system of complex differentialequations, the point z r is an isolated singularity. At first sight, it seemed to us that thesingularity is a first order pole, since functions of the form f (cid:48) f , where f has a zero of anyorder at some point z do have a first order pole at z . Unfortunately, the extra termsthat are prevalent in the coefficients, i.e. im (cid:48) ( z ) z and exp( − im ( z ) z ) involve squareroots, which results in f and g not being holomorphic in a punctured disk around z r ,since the complex square root has two branches and therefore it is only holomorphic in adisk with one half-axis removed. Therefore, all theorems about existence and structureof solutions can not be applied. Existing research in this case, such as Sutherlandand Yewchuck (2004), who investigate propagation of gravity waves through a layer ofsudden reduced or vanishing stratification, suggest some sort of ”wave tunnelling” (aterm coined by the comparison to quantum tunnelling of electrons in quantum physics)through this region, dependent on the wavelength of the incident wave. The findings fromour multi-layer method confirm those results also for continuous transitions to a lowervalue of the Brunt-V¨ais¨al¨a frequency (see section 4.2). Hence, an intensive investigationof the behaviour of gravity waves near reflection layers would require a scale analysis fordifferent regimes of vertical wavelengths. This task is taken on by the authors, but liesoutside the scope of this paper. This subsection is dedicated to the numerical integration of equation (38). To this end,proper boundary conditions are needed. Since we are interested in waves that are initiallytravelling upwards and encountering a non-uniform stratification over a confined regionand eventually reaching a region of uniform stratification again, we require that there is13igure 2: The left panel shows the multi-layer method with J = 512 steps and theright panel shows the values of the transmission coefficient computed from thenumerical evaluation of the limit approach. With bare eye, they are indis-tinguishable. At every single point in the domain, the error is smaller than10 − only a wave incident on the non-uniform region from below. Hence, above this region,there is no wave that is travelling downwards. Moreover, we are only interested in ratiosbetween the incident and the transmitted wave. Therefore we are free to choose thevalue of the incident or the transmitted wave, since the equation is linear. To be moreprecise, if z t is the top of the non-uniform region, then we take the boundary conditions A ( z t ) = (cid:18) A ( z t ) B ( z t ) (cid:19) = (cid:18) (cid:19) . (39)We want to analyse the error between the limit solution and the multi-layer solution,which will give strong evidence of the correctness of our method. As model setup, weagain choose the linearly increasing profile (31). Our analysis consists of two parts. Thefirst one is an error computation over a large domain of wave numbers and frequencieswhile keeping the number of discretisation levels constant at J = 512 to show that theerror is small over the whole wave number-frequency-domain. The second one choosesseveral specific points in this domain and tracks the error for an increasing number oflevels J , up to J = 10 .The results of the first part can be seen in figure 2. The pictures show the transmissioncoefficient for a certain domain of wavelengths and frequencies, similar to figure 1. But we14 a) (b) (c) Figure 3: Relative error for profile (31) with N t = 2 N b , ω = N b √ and different horizontalwavelengths λ x : λ x = ∆ z in panel (a), λ x = 2 ∆ z in (b) and λ x = 10 ∆ z in(c). The mean slope of the different plots is nearly the same, namely µ ≈ − J = 512 levels, the right panel shows the transmission coefficientcomputed from solving equation (38) numerically. It is impossible to spot any differencebetween the two frames. Computing the relative error yields the estimatemax ω,λ x | T C d ( ω, λ x ) − T C l ( ω, λ x ) || T C l ( ω, λ x ) | < · − (40)For the second analysis, we fix specific wave parameters, i.e, a pair ( λ x, , ω ) of wave-length and frequency, and analyse how the relative error develops for increasing J . Inparticular, we perform the calculation for three different wavelength-frequency-pairs. Wechoose ω = N b √ for all three cases and have a look at the wavelengths ∆ z , ∆ z and 10 ∆ z .The results can be seen in figure 3.We computed the relative error of the limit solution and the discrete solution for severalnumbers of layers J , that were logarithmically spaced between 10 and 10 . For any twoadjacent points, we computed the slope in the log-log diagram and for every wavelength,the mean and the standard deviation of all computed slopes. We found the mean slopes15o be µ a = − . ± . , µ b = − . ± . µ c = − . ± . J = 100 layers, the relative error is about 10 − .Based on this, we will use J = 128 layers for the forthcoming computations. Thisguarantees fast run times as well as results that are sufficiently accurate. In this section, we present some results that are obtained for several stratification profiles.First, we focus on the linearly increasing profile we already used for some model com-putations. Afterwards, we show that the model also supports wave tunnelling that wasalready described by Sutherland and Yewchuck (2004). At the end, we present resultsfor a profile that has characteristics of a real tropopause. All profiles share a commonstructure, namely that we have a region of non-constant stratification of depth ∆ z thathas a region of constant stratification with value N b below and with value N t above it.This is about what we can observe in the atmosphere: The stratification changes rapidlyin the tropopause and is nearly constant in the free troposphere and the stratosphere.We will focus on a frequency range from 0 to N b , since waves with frequencies largerthan N b are evanescent. The wavelength spectrum differs in the various examples butranges in general from ∆ z to 100 ∆ z . We use a parameter space of 300 horizontal wavenumbers and 300 frequencies, resulting in a total of 90000 grid points. The number oflayers equals 128. The computations are performed on the first author’s office computerwith a standard Intel ® Core ™ i7-3770 CPU and 8 GB RAM. The software we used is Matlab . The computation times lie within a range of 70 to 80 seconds when computedon a single core. Compared to the numerical method of Nault and Sutherland (2007),who need about 1 day to simulate 300 ×
300 parameters on a ”typical desktop computer”at that time, this is a decrease in computation time by a factor of about 1000. Evenwhen considering the slightly higher clock rate and RAM, the multi-layer method ismuch more efficient.
The first case for which we will present results is the case of the linearly increasingstratification, as we defined earlier in equation (31). We focus on three regimes ofhorizontal wavelengths: comparable to ∆ z , longer than ∆ z and shorter than ∆ z . Recallthat we are using the Boussinesq approximation, hence ∆ z , as it is also the scale ofvariation of N , is small compared to the density scale height H ρ . So wavelengths thatare small and comparable to ∆ z are also small compared to H ρ . This is the regime wherethe classical WKB theory is applicable. Ray theory, that is based on WKB assumptions,predicts perfect transmission for those waves in a linearly increasing profile and this is16igure 4: Here we see the results for a stratification profile, that increases linearly over afinite region of depth ∆ z . The right panel is an enlargement of the rightmostpart of the left one. The 3 thick lines represent constant vertical wavelengthin the bottom layer (it changes while the wave propagates upwards).exactly what we are able to find. Even for large horizontal wavelengths, the transmissionis high, at least up to a certain point. As we can see in figure 4, there is stronger reflectionwhen we are moving to the right and to the bottom in the figure, that means that thewave frequency gets closer to the lower value N b and the horizontal wavelength (andeventually the vertical wavelength) is growing. Figure 4 grants us a closer view intothe area of interest, together with some additional information. The solid lines are linesof constant vertical wavelength, which are determined by the gravity wave dispersionrelation (10) (they are hyperbolas in the frequency-wave number-space). We see thatthe vertical wavelength increases by moving to the right and to the bottom in our domain.If the vertical wavelength exceeds ∆ z by about an order of magnitude, waves start totransmit worse and worse.Waves whose wavelengths are small, or at least comparable to the scale of variation inthe Brunt-V¨ais¨al¨a-frequency, adapt smoothly to those changes and have a high transmis-sion, while for waves with large wavelengths, the change still happens abruptly. Wavesthat have frequencies close to the Brunt-V¨ais¨al¨a frequency are almost purely horizontaland in the limit of ω → N b for fixed λ x , the transmission eventually gets zero. Thisseems reasonable, since there is no vertical wave structure and hence no vertical energytransport. In the limit λ x → ∞ (or equivalently k →
0) for fixed ω < N b , we can finda closed formula for the matrix product in (16) and hence a formula for the amplituderatio, which coincides with the classical transmission coefficient of a two-layer model: T C L = A t A i = 2 m i m i + m t , (41)17igure 5: The figure shows the transmission coefficients for a profile that has a region ofdecreased stratification. We see that waves, whose frequency is larger than0 . N b can transmit, if their wavelength is long compared to the region ofdecreased stratification.where the indices describe incident and transmitted wave properties respectively. Weadapt the slightly different notation of a transmission coefficient as the ratio of the waveamplitudes due to the classical work that was done in this area, e.g., by Eliassen andPalm (1961). An analytic execution of these limits can be found in the appendix. We consider now a case where the stratification drops from some N b to a value N d < N b and eventually increases again back to N b : N ( z ) = N b , z < z b N b + z − z b z d − z b ( N d − N b ) , z b < z ≤ z d N d , z d ≤ z ≤ z d N d + z − z d z t − z d ( N b − N d ) , z d < z ≤ z t N b , z t < z. (42)In the example we present, z d − z b = 0 . ∆ z = z t − z d and N d = 0 . N b . The resultsfound for this case are very different from what the classical theory tells us. Ray theorypredicts that waves reflect totally from a layer, where ω ≥ N . However, when there18s only a finite small region where ω ≥ N holds, wave propagation through this regionis possible under certain conditions. Sutherland and Yewchuck (2004) described thisphenomenon for a sharp drop to a weak or even vanishing stratification, and our resultsshow that tunnelling also exists in the case of a ”smooth” transition. In figure 5, one cansee the transmission coefficient for profile (42). In this particular example, ray theorypredicts that every wave with frequency ω ≥ . N b would fully reflect from this layer,but we can observe that if the wavelength is large compared to the extent of the regionwith weak stratification, it is possible to obtain high wave transmission.As already mentioned in section 3, this case is hard to deal with analytically and nu-merically. As we will see later in section 5, numerical computations of the full Boussinesqequations show the existence of wave tunnelling (although there are some other difficul-ties with this case). Moreover, there are lab experiments, for example by Sutherland andYewchuck (2004), to further support these findings. However, the limit ODE found insection 3 cannot be solved numerically due to a singularity in the coefficients. Althoughthere is no evidence that the results are incorrect, the results should mathematically stillbe taken with a grain of salt. Since we are ultimately interested in the behaviour of atmospheric gravity waves and theirinteraction with the tropopause, we now want to consider a ”realistic” tropopause profile.The stratification is constant with a value N b below the tropopause. At the temperatureinversion layer, the Brunt-V¨ais¨al¨a frequency has a very sharp increase to a peak value N p , almost like a jump, followed by a relaxation to a value N t with N b < N t < N p that is the constant value of the stratification in the stratosphere. We realize this by apiecewise-defined continuous function: N ( z ) = N b , z < z b N b + z − z b z p − z b ( N p − N b ) , z b < z ≤ z p az + bz + c, z p < z ≤ z t N t , z t < z, (43)where a = N p − N t ( z p − z t ) , b = − z t N p − N t ( z p − z t ) and c = N t + z t N p − N t ( z p − z t ) . The values were chosensuch that the profile is continuous at z p and z t and differentiable at z t . In the examplewe show here, the values were set to be z p − z t = 0 . ∆ z , N p = 3 N b and N t = 2 N b . Theprofile as well as results for the transmission coefficient can be seen in figure 6.For the limits ω → N b for fixed λ x and λ x → ∞ for fixed ω , we have the exact samebehaviour as in the linearly increasing case, which is no surprise, since in the first limit,we still have no vertical energy flux and in the second limit we again approach the two-layer model. The rest of the picture however gives some interesting insights. By makingthe sharp increase asymptotically thin, i.e., making it a (discontinuous) jump, we obtain,for wavelengths comparable to ∆ z , a composition of the transmission coefficient for atwo-layer model (that describes the jump) and the one for the smooth profile that followsafter the jump. 19igure 6: Values for the transmission coefficient for a tropopause profile that can beeseen in the right panel. In the limit for long waves, the transmission coefficientapproaches again the two-layer solution. For moderately long waves, we cansee a combination of 2 effects: The sharp increase, which is almost like a jumpand blocks a part of the waves and the smooth relaxation afterwards that hashigh transmission. We also observe that in the classical WKB regime, thetransmission is still very high In this section, results for transmission coefficients retrieved from numerical simulationsof atmospheric motion are presented, discussed and compared to the results from themulti-layer method introduced above.We use two different models to tighten the theoretical findings. The first one is EU-LAG, a Eulerian/semi-Lagrangian fluid solver described in Prusa et al. (2008), the sec-ond one is PincFloit (Pseudo-incompressible flow solver with implicit turbulence model),which was developed by Rieper et al. (2013).
The calculations are done on a two-dimensional x - z -domain with periodic boundariesin x -direction. Both models, EULAG and PincFloit, are run in Boussinesq mode andare set up in the following way in order to resemble the scenarios for wave transmissiondiscussed above the most. An absorption layer is located at the bottom of the domain,where the flow is relaxed towards a state that is here chosen to be a plane wave field20hich oscillates in space and time. In order to do this properly, the Brunt-V¨ais¨al¨afrequency needs to be constant in this sponge layer. This results in the excitation ofa plane wave that can propagate freely above the sponge layer that covers the lowest40% of the domain, until it eventually reaches a region of non-uniform stratification, thetropopause, which is defined as in the three cases discussed in section 4, where it getsreflected and partially transmitted. The reflected part gets soaked up by the bottomsponge while the transmitted part travels in a stratospheric region, i.e. a region where N is again constant before getting damped away by the top sponge. This sponge coversthe last 20% of the domain.The initial wave amplitude has to be chosen very small, since the theoretical workuses the linearised Boussinesq equations. Moreover, the waves should not break duringthe simulations. Hence we initialise the amplitude as b = a N b m b (44)where a < b s = N b m b , found in, e.g., Achatz et al. (2010) is reached. The quantities of the wave fieldon the bottom layer can be derived via the polarisation relations for Boussinesq waves,which can be found in, e.g., Achatz et al. (2010) or Sutherland (2010). u = b m b k ωN b cos( kx + m b z − π − ωt ) (45) w = b ωN b cos( kx + m b z + π − ωt ) (46)The top sponge is just an absorption layer, which damps the wave and relaxes to asteady, hydrostatic background.The region of non-uniform stratification starts at 60% of the domain and is in generalresolved with 100 grid points. For this value, the theoretical results are very acurate(absolute error is smaller than 10 − , as seen in section 3.3). In the reflection layer case,we changed the size of the tropopause rather than the wavelength since we wanted tosimulate large ratios between wavelength and tropopause depth, so it was more efficientto shrink the extent of the tropopause. The fluid solver EULAG, described by Prusa et al. (2008), is used in the Boussinesqframework (anelastic setup with special choice for the hydrostatic basic state in density ρ and in potential temperature θ , namely constant). The calculations are done ona two-dimensional model domain that extends up to 8500 m in the vertical and has ahorizontal width ranging from 1000 m to 3000 m, depending on the wavelength that isprescribed. The resolution is chosen to be around 10 m in x and z direction. The verticalextent of the domain is limited by the Boussinesq assumption to a value below whichthe basic state profiles remain physically meaningful.21he transmission coefficient for each time step t is calculated as T C ( t ) = m | A | m | A | = m (max | w ss ( t ) | ) m (max | w ts ( t ) | ) . (47)The transmitted amplitude A is determined by finding the maximum absolute verticalwind speed in a layer above the region of interest (stratosphere: w ss ) while the amplitude A is retrieved from a reference simulation with constant stratification and taken as themaximum absolute vertical wind speed w ts just below the tropopause altitude. Thedepth of these regions is roughly 1 km for the cases shown here. This is sufficient sincehorizontal and vertical wavelengths are small enough to ensure a maximum or minimumin this box. Since we expect a static situation for the transmission of a plane wavethrough a layer of non-uniform stratification, the simulations are done up to a simulationtime of 15 h when the simulated flow has stabilized and shows little temporal variationin the quantities of interest. We then take the mean value of T C for the
T C calculatedat simulation times 8 h , ,
10 h , ...,
15 h.
PincFloit has a built-in switch for a Boussinesq atmosphere, i.e. constant backgrounddensity, but due to extensions implemented by B¨ol¨oni et al. (2016), it is possible to havea non-constant profile for the Brunt-V¨ais¨al¨a frequency. Hence it was an easy task toimplement the test cases of this manuscript into the solver. Moreover, as can be seen forexample in B¨ol¨oni et al. (2016) and Schlutow et al. (2017), PincFloit gives very robustresults in simulations of atmospheric gravity waves.The model domain covers 10000 m in the vertical and one wavelength in the horizontaldirection, which lies between 1000 m and 3000 m. The vertical resolution is 10 m, andthe horizontal resolution is 25 m.In all simulated runs, we could observe a stable steady state for several hours of simu-lation time, which allows us a very good computation of the transmission coefficient. Theformula for the computation is, of course, the same as for the EULAG simulations, thechoice of incident and transmitted amplitude however is a little bit different. Nonethe-less, no method is superior over the other and both yield equally good results. As initialamplitude, we take the mean amplitude of the excited wave in the middle of the bot-tom sponge layer, as transmitted amplitude we take the mean amplitude of the wavein the stratospheric region, i.e. between 7000 and 8000 m. Moreover, the transmittedamplitude is computed every 10 min of simulated time and the transmission coefficientis taken as mean over all computed values after reaching the steady state.
In this section we show exemplary simulations that cover a portion of the informationobtained in the transmission coefficient figures (see figures 4, 5 and 6) that have beenpresented before, as well as snapshots of the wave field after reaching a steady state.This gives a better intuition for wave transmission and reflection.22 z λ x In table 1, we see the comparison of the transmission coefficients obtained from the multi-layer method on one hand and from the numerical simulations of PincFloit and EULAGon the other hand. In general, the values from both methods are in good agreement withthe theoretical results.Generally speaking, the numerical transmission coefficients are a bit lower than thetheoretical ones. This has two reasons. The first one is numerical dissipation whichcauses a slight decrease of the amplitude over the course of the domain. This is becausethe wave is initialised in the bottom part of the domain and then freely propagatesupward. The second reason is concerned with the method we compute incident andreflected amplitude. Due to the nature of the sponge, there are some fluctuations in theinitialised wave amplitude and hence also in the transmitted amplitude. We take careof this by averaging over a larger area, but this produces minor errors in the value forthe transmission coefficient. Plots of the simulated wind field however show that there islittle to no reflection at the tropopause in the cases from table 1, hence the transmissionseems to be even closer to the theoretically predicted values.
The simulation setup in this case uses fixed values for horizontal and vertical wave-length, but changes the depth of the tropopause, as this is numerically more convenient.Technically, we have the same wave frequency, but change the ratio between horizon-tal/vertical wavelength and tropopause depth, so this corresponds to a vertical slice inthe left panel of figure 5. By choosing λ x = λ z , as we did here, this results in a fre-quency ω = √ N b ≈ . N b . For tropopause depths that are comparable to λ x , weexpect strong reflection, while for a very short tropopause, we should obtain a wavetunnelling effect. The results can be seen in table 2. The PincFloit simulations used λ x = λ z = 1000m, while EULAG used λ x = λ z = 2000m.Both models match the theoretical prediction very accurately. While the transmission23igure 7: Vertical wind field of a PincFloit simulation for the linear case with a horizontalwavelength of 2500 m and an incident vertical wavelength of 1000 m. We seeno chequerboard pattern, but only a refraction of the wave due to the changein vertical wavelength. This implies that almost no wave reflection occurs.The smaller wave amplitude in the stratosphere (i.e. between 7000 m and8000 m) is a result of the stronger stratification in this altitude. The axes arenormalised by 1000 m. ∆ z λ z λ z = λ x and varying tropopause depth ∆ z .This corresponds to a vertical cut along ω ≈ . Table 3 shows a comparison of values for the transmission coefficient obtained from themulti-layer method and values computed from PincFloit and EULAG simulations. Asin the two other cases, we can see a very good agreement between theory and numerics.These cases were particularly interesting since there is only a partial reflection of thewave. Both codes were able to reproduce this phenomenon. A snapshot of how a wavefield looks like in this case can be seen in figure 9.25 z λ x The simulations deliver very accurate results. The qualitative as well as the quantitativebehaviour could be reproduced in every single of the atmospheric test cases with differentparameters. Even wave tunnelling could be observed. This leads to the conclusion thatthe multi-layer method is capable of predicting the transmission and reflection of gravitywaves on a satisfying level. Since the computation is very fast and efficient, this could finduse as a black box in numerical weather models, whose resolution is in general too coarseto resolve all gravity wave structures. The prediction of how much of the wave energytransmits through the tropopause can lead to a better parametrisation of gravity wavesin the middle and upper atmosphere, and could also explain why some of the waves thatare supposed to have broken at a certain height are still present and stable. Moreover,stronger temperature inversions in the tropopause lead to larger wave reflection with acorresponding downdraught of energy which could affect the tropospheric balance.
This section is dedicated to limitations as well as the possible extensions of the model.Since we are using the Boussinesq approximation, there are clearly some bounds onthe applicability in case the fluid density is changing significantly. But there is alsoroom for improvement and generalisation. Background wind is an important factor inthe atmosphere as well as wave packets, i.e. locally confined wave movements whoseenvelope is moving with the group velocity.
Waves with a frequency higher than the Brunt-V¨ais¨al¨a frequency are evanescent. How-ever, we have seen in Section 4.2 that if the region over which the waves are evanescentis small compared to the wavelengths, the waves can survive and eventually propagateagain. Hence one would expect a similar behaviour if waves are excited with a frequency26igher than the Brunt-V¨ais¨al¨a frequency, but after a short extent reach a region whereit gets and stays larger than their own frequency. Although there is nothing that speaksagainst this hypothesis, our approach as it is presented here, is not capable of reproduc-ing this effect. We are assuming an infinitely extended domain of constant stratificationabove and below the non-uniform region and that incident waves come from ”far away”.Hence a wave that is initially evanescent will completely vanish (i.e. amplitude A = 0)when reaching the non-uniform region.A workaround for this would be to assume a stratification profile that has an infinitebottom layer that allows for wave propagation with a certain frequency ω , followedby a jump to a confined region where the Brunt-V¨ais¨al¨a frequency is smaller than ω ,which in turn is succeeded by the profile we are interested in. The relevant amplitudefor the transmission coefficient computation would then be the upward amplitude in theconfined region of reduced stratification. This is possible since the multi-layer methodkeeps track of upward and downward propagating wave amplitudes. Including background wind
All previous computations are done for an atmosphere at rest. The case with backgroundwind can be more difficult, since not wave energy, but rather wave action is conserved.However, since we are interested in the wave action fluxes, the formula for the transmis-sion coefficient does not change. We have to distinguish two different cases: constantand non-constant background wind.In the case of constant background wind U = U ≡ const., the results are pretty muchthe same as with no background wind, except that one has to use the intrinsic frequencyˆ ω = ω − kU instead of the extrinsic frequency ω (which are in fact equal if there is nobackground wind). Waves can then propagate for ˆ ω ∈ (0 , N b ). Apart from that smallchange, the multi-layer method, as it is presented here, can be applied without furthermodification.If the background wind is changing with height, it will have an impact on the verticalpropagation of the gravity waves just as much as the non-uniform stratification has.In this case, we have to take into account the change in absolute frequency over theshear region as well as the curvature of the mean background wind, which is the rateof change of the shear. Moreover, we can have reflection layers, similar to the case ofweakening stratification, and also a new phenomenon, called critical layers, arises. Theyoccur at locations where the absolute frequency vanishes, i.e., when the horizontal phasespeed matches the background wind speed. Around critical layers, non-linear effectsbecome more and more important. For example, wave-mean flow interaction, where thewave deposits energy to or draws it from the mean flow, is an important factor that ourapproach does not account for. But apart from critical layers, which are also an issuein direct numerical simulations, the implementation of non-constant background wind ispossible, but will be postponed to a companion paper Puetz and Klein (2018).27 ave packets Wave packets are an important part of atmospheric gravity wave analysis, since a long-lasting source of wave generation, such as steady flow over a mountain ridge, is a rareevent. A wave packet can be seen as an amplitude-modulated plane wave, with almostcompact support, i.e., the envelope vanishes or approaches zero outside a specific region.Mathematically, this can be described as the superposition of infinitely many plane waveswith different wave numbers. They destructively interfere almost everywhere except fora confined region in which the interference is constructive. Fourier transformation canbe used to break the wave packet down into the different wave numbers with respectiveamplitudes and phase shifts.It is possible to apply the method introduced in this paper to wave packets by ap-proximating them as a superposition of finitely many plane waves with correspondingamplitudes. Since it is not only about applying the multi-layer method for some selectedwave parameters, but also requires some additional set-up, we will not include the de-tails here. Alongside the inclusion of background wind, this will be the main issue of acompanion paper Puetz and Klein (2018).
We developed a method with which we can compute the transmission of gravity wavesthrough a finite region of non-uniform stratification by modelling this region as a multi-layer fluid which has a uniform stratification in each layer. The solutions we found foreach layer were matched across the fluid interfaces and we were able to relate incidentand transmitted wave. The method is applicable to any stratification profile. Moreover,it is able to keep track of the upward and downward propagating wave amplitude at anypoint in the domain.We were also able to find the limit for the number of layers tending to infinity, lead-ing us to a system of ordinary differential equations for the amplitudes of upward anddownward propagating wave. For stratification profiles not including a reflection layer,this ODE system could be solved numerically in a very efficient way by using a basicintegration scheme. An error analysis showed that the multi-layer method converges tothe limit solution.Numerical simulations of the full Boussinesq equations were carried out to support ourtheoretical findings. We found the results of the simulations to be in good agreement tothe predictions from the multi-layer method.
Acknowledgements
The data for this paper are available upon request from the authors. The authorsthank the German Research Foundation (DFG) for support through the research unitMulti-Scale Dynamics of Gravity Waves (MS-GWaves) and through grants KL611/24-1,KL611/25-1 and SP1163/5-1. They also thank the group of Prof. Ulrich Achatz for the28rovision of PincFloit. 29 eferences
Ulrich Achatz, Rupert Klein, and Fabian Senf. Gravity waves, scale asymptotics andthe pseudo-incompressible equations.
J. Fluid Mech. , 2010.John R. Booker and Francis P. Bretherton. The critical layer for internal gravity wavesin a shear flow.
J. Fluid Mech. , 1967.G. L. Brown and Bruce R. Sutherland. Internal wave tunnelling through non-uniformlystratified shear flow.
Atmosphere-Ocean , , 2007.Gergely B¨ol¨oni, Bruno RIbstein, Jewgenija Muraschko, Christine Sgoff, Junhong Wei,and Ulrich Achatz. The interaction between atmospheric gravity waves and large scaleflows: an efficient description beyond the nonacceleration paradigm. J. Atmos. Sci. ,2016.Oliver B¨uhler.
Waves and mean flows . Cambridge University Press, 2009.Edwin F. Danielsen and Rainer Bleck. Tropospheric and stratospheric ducting of sta-tionary mountain lee waves.
J. Atmos. Sci. , , 1970.P. G. Drazin and W. H. Reid. Hydrodynamic Stability . Cambridge University Press,1981.Arndt Eliassen and Enok Palm. On the transfer of energy in stationary mountain waves.
Geofys. Publ. , :1–23, 1961.Luis Lara. A numerical method for solving a system of nonautonomous linear ordinarydifferential equations. Appl. Math. Comput. , 2004.Matthieu J. Mercier, Nicolas B. Garnier, and Thierry Dauxois. Reflection and diffractionof internal waves analyzed by the hilbert transform.
Phys. of fluids , 20, 2008.J. T. Nault and Bruce R. Sutherland. Internal wave transmission in nonuniform flows.
Phys. Fluids , (016601), 2007.J. M. Prusa, P. K. Smolarkiewicz, and A. A. Wyszogrodzki. Eulag, a compuationalmodel for multiscale flows. J. Comput. Fluids , 2008.Felix Rieper, Stefan Hickel, and Ulrich Achatz. A conservative integration of the pseudo-incompressible equations with implicit turbulence parametrization.
Mon. WeatherRev. , 2013.Mark Schlutow, Rupert Klein, and Ulrich Achatz. Finite-amplitude gravity waves in theatmosphere: travelling wave solutions.
J. Fluid Mech. , 826, 2017.Richard S. Scorer. The theory of waves in the lee of mountains.
Quart. J. R. Met. Soc. , , 1949. 30ohn Stalker. Complex Analysis . Birkh¨auser, 1998.Bruce R. Sutherland.
Internal Gravity Waves . Cambridge University Press, 2010.Bruce R. Sutherland and Kerianne Yewchuck. Internal wave tunnelling.
J. Fluid Mech. ,2004.Gerald Teschl.
Ordinary Differential Equations and Dynamical Systems , volume 140 of
Graduate studies in mathematics . American Mathematical Society, 2012.Christopher P¨utz and Rupert Klein. Initiation of ray tracing models: Evolution of small-amplitude gravity wave packets in non-uniform background.
Theor. Comp. Fluid Dyn. ,submitted. 31 ppendix A
Here is the proper execution of the limit process mentioned in equation (36): c j − (cid:18) m ( z j ) m ( z j +1 ) + 1 (cid:19) exp ( i ( m ( z j ) − m ( z j +1 )) z j ) −
1= 12 (cid:18) m ( z j ) m ( z j + h ) + 1 (cid:19) exp ( i ( m ( z j ) − m ( z j + h )) z j ) −
1= 12 (cid:18) m ( z j ) m ( z j ) + hm (cid:48) ( z j ) + o ( h ) + 1 (cid:19) · exp( ih ( m (cid:48) ( z j ) + o ( h )) z j ) −
1= 12 (cid:18) m ( z j ) m ( z j ) + hm (cid:48) ( z j ) + o ( h ) + 1 (cid:19) · (cid:0) ihm (cid:48) ( z j ) z j + o ( h ) (cid:1) −
1= 12 (cid:18) m ( z j ) + hm (cid:48) ( z j ) m ( z j ) + hm (cid:48) ( z j ) + o ( h ) (cid:19) (cid:0) − ihm (cid:48) ( z j ) z j + o ( h ) (cid:1) − (cid:18) − hm (cid:48) ( z j )2( m ( z j ) + hm (cid:48) ( z j ) + o ( h )) (cid:19) − h (cid:18) m ( z j ) + hm (cid:48) ( z j )2( m ( z j ) + hm (cid:48) ( z j ) + o ( h )) (cid:19) im (cid:48) ( z j ) z j + o ( h )= h (cid:18) − m (cid:48) ( z j )2( m ( z j ) + hm (cid:48) ( z j ) + o ( h )) (cid:19) − h (cid:18) m ( z j ) + hm (cid:48) ( z j )2( m ( z j ) + hm (cid:48) ( z j ) + o ( h )) (cid:19) im (cid:48) ( z j ) z j + o ( h ) . (48)We are able to compute the limit of c j − h f ( z ) := lim h → c j − h = − m (cid:48) ( z )2 m ( z ) − im (cid:48) ( z ) z. (49)A similar derivation can be made for d j h . Appendix B
We investigate the limit ω → N b for fixed k and the limit k → ω . For thefirst limit, we have a look at the following:lim ω → N b m j = lim ω → N b − k (cid:115) N j ω − − k (cid:115) N j N b − m j , (50)If j (cid:54) = 1, then N j (cid:54) = N b (remember: linear increasing profile). Hence, ˜ m j (cid:54) = 0. Sincedet M j = m j m j +1 , M j is regular for J (cid:54) = 1. For j = 1, N = N b , hence ˜ m = 0 and32herefore det M = 0. Since the determinant is multiplicative, the determinant of thematrix product (16) is 0, and by the definition of the transmission coefficient (28), thetransmission is also 0.For the second limit, we write m j = k (cid:115) N j ω − k ˆ m j (51)and consider the matrix entries (14) in the limit k → c j := lim k → c j = lim k → (cid:18) m j m j +1 + 1 (cid:19) exp ( i ( m j − m j +1 z j ))= lim k → (cid:18) k ˆ m j k ˆ m j +1 + 1 (cid:19) exp ( ik ( ˆ m j − ˆ m j +1 z j ))= lim k → (cid:18) ˆ m j ˆ m j +1 + 1 (cid:19) (1 + ik ( ˆ m j − ˆ m j +1 z j ) + o ( k ))= 12 (cid:18) ˆ m j ˆ m j +1 + 1 (cid:19) (52)In a similar fashion, we follow that˜ d j := lim k → d j = 12 (cid:18) ˆ m j ˆ m j +1 − (cid:19) . (53)If we assume now, that there is no reflection layer, i.e. m j ∈ R for all j , then the matrices (cid:102) M j := lim k → M j are real and symmetric. Hence the product of any two matrices ofthis kind is again symmetric. We have a closer look at (cid:102) M j (cid:102) M j − . The diagonal entriesare ˜ c j ˜ c j − + ˜ d j ˜ d j − = 14 (cid:18) ˆ m j ˆ m j +1 + 1 (cid:19) (cid:18) ˆ m j − ˆ m j + 1 (cid:19) + 14 (cid:18) ˆ m j ˆ m j +1 − (cid:19) (cid:18) ˆ m j − ˆ m j − (cid:19) = 14 (cid:18) ˆ m j − ˆ m j +1 + ˆ m j ˆ m j +1 + ˆ m j − ˆ m j + 1+ ˆ m j − ˆ m j +1 − ˆ m j ˆ m j +1 − ˆ m j − ˆ m j + 1 (cid:19) = 12 (cid:18) ˆ m j − ˆ m j +1 + 1 (cid:19) . (54)33or the off-diagonal entries, we have˜ c j ˜ d j − + ˜ d j ˜ c j − = 14 (cid:18) ˆ m j ˆ m j +1 + 1 (cid:19) (cid:18) ˆ m j − ˆ m j − (cid:19) + 14 (cid:18) ˆ m j ˆ m j +1 − (cid:19) (cid:18) ˆ m j − ˆ m j + 1 (cid:19) = 14 (cid:18) ˆ m j − ˆ m j +1 − ˆ m j ˆ m j +1 + ˆ m j − ˆ m j −
1+ ˆ m j − ˆ m j +1 + ˆ m j ˆ m j +1 − ˆ m j − ˆ m j − (cid:19) = 12 (cid:18) ˆ m j − ˆ m j +1 − (cid:19) . (55)A simple induction argument shows that the entries of (cid:102) M = (cid:18) ˜ c ˜ d ˜ d ˜ c (cid:19) := (cid:89) j = J (cid:102) M j (56)are ˜ c = 12 (cid:18) ˆ m ˆ m J +1 + 1 (cid:19) (57)˜ d = 12 (cid:18) ˆ m ˆ m J +1 − (cid:19) . (58)The relation between the amplitudes is the same as in (19), i.e. A J +1 A = det( (cid:102) M )˜ d = m m J +1 (cid:18) (cid:18) ˆ m + ˆ m J +1 ˆ m J +1 (cid:19)(cid:19) − = 2 ˆ m ˆ m + ˆ m J +1 ,,