A perfectly matched layer formulation for the nonlinear ... - CiteSeerX

A perfectly matched layer formulation for the nonlinear shallow water equations models: The split equation approach I.M. Navon Department of Mathematics and Computational Science and Information Technology Florida State University Tallahassee, FL 32306-4130 B. Neta Naval Postgraduate School Department of Mathematics Monterey, CA 93943 M.Y. Hussaini Computational Science and Information Technology Florida State University Tallahassee, FL 32306-4120 July 22, 2001

1

Abstract Perfectly matched layer (PML) equations for the treatment of boundary conditions are constructed for the two dimensional linearized shallow-water equations. The method uses the splitting technique, i.e. the absorbing layer equations are obtained by splitting the governing equations in the coordinate directions and absorbing coefficients are introduced in each split equation. The shallow water equations are discretized using the explicit Miller Pierce finite difference approach. The method was tested and the transparency of the boundaries was demonstrated on the nonlinear shallow-water equations including the Coriolis factor on a limited-area rectangular domain for a Gaussian pulse with convective mean flow. Different scenarios involving various values of the mean flow speeds, depth of the PML layer and different PML absorbtion coefficients are tested. The numerical results indicate that outgoing waves are leaving the domain without perturbing the flow in the physical domain. No filters were used in the numerical experiments which included both a stationary as well as a convected Gaussian. The exterior domain was ended using well posed boundary conditions for the shallow water equations. The L2 error in the height field between the PML layer domain and a reference solution along a line located inside the interior domain was computed also for the case using only characteristic boundary conditions as well as the second order Engquist and Majda boundary conditions for comparison purposes. The efficacy of the PML scheme over the other two aforementioned schemes for the nonlinear shallow water equations was confirmed . PML layers of increasing thickness yielded better results than either the characteristic treatment or the second order Majda Engquist absorbing boundary conditions.

2

1

Introduction

Absorbing boundary conditions and the treatment of lateral boundaries have been a perennial problem in the formulation of lateral boundary conditions since the early days of numerical weather prediction. In a limited-area model the lateral edges are not physical boundaries of the flow but constitute artificial constraints imposed by computational considerations. Hence they do not have a physical counterpart. We must impose conditions at these artificial boundaries in order to solve the problem in an efficient and accurate manner. Several good review papers were written on the topic of both actual and artificial boundary conditions, see alsoa special issue (Givoli and Harari, eds.) of Computer Methods in Applied Mechanics and Engineering dedicated to exterior problems of wave propagation. A general discussion was presented by Turkel [48] while a review of non-reflecting boundary conditions presenting the issue in a uniform manner is that of Givoli [15]. Nonreflecting boundary conditions aim to achieve some of the following goals for a problem defined in a computational domain Ω whose artificial boundary is Γ: • The problem in Ω together with the boundary condition on Γ should be mathematically well-posed and be a good approximation of the original problem in the infinite domain. • The boundary condition on Γ should be highly compatible with the numerical discretization scheme used in the interior computational domain Ω. It is also important that the numerical method used together with the boundary condition should result in a numerically stable scheme, while the amount of spurious reflection generated by the boundary condition on Γ should be small. Another consideration is that the implementation of the boundary condition should not lead to a large computational effort. A recent comprehensive and extensive review of artificial boundary conditions was presented by Tsynkov [47]. He shows that any algorithm for getting artificial boundary conditions can be viewed as a compromise between locality and practical efficacy but implying insufficient accuracy. On the other hand highly accurate nonlocal boundary techniques result in impractical algorithms. Gustafson and Sundstrom [17], treated the linearized shallow-water equations using a Fourier transform in all but one Cartesian direction and used energy estimates to treat the resulting one dimensional systems to generate well-posed boundary conditions for the half-space problem. Tsynkov [47] also shows that the method of introducing sponge layers near the computational boundaries occupies an intermediate position between local and nonlocal approaches This is since on one hand, there are no global integral relations along the boundary and when numerical simulations are conducted, the model equations inside the sponge layer are solved using a numerical method close or identical to the one employed inside the computational domain. On the other hand a certain amount of nonlocality persists since the computational domain is artificially magnified. Two basic approaches have been used for dealing with limited-area boundary conditions. Simple boundary conditions for hyperbolic equations consist in having a specification requiring that the characteristics leaving from the interior domain be zero. In addition 3

one must combine the absorbing boundary conditions with the given inflow conditions. A strictly mathematical treatment of absorbing boundary conditions (ABC) for hyperbolic equations was presented by Engquist and Majda [13, 14] based on pseudo-differential operators subsequently expanded to non-local operators to get local well posed ABC’s. These methods can be viewed as a generalization of the Sommerfeld radiation condition and the characteristic approach. A second alternative approach is the use of sponge or damping layers to damp out disturbances prior to their reaching the artificial boundary. A variation of this process is to construct a layer where the outgoing waves will slow down rather than decay. Hence, the waves will not reflect back into the limited-area forecast domain of interest except at very late times (see Perkey and Kreitzberg [43], Davies [8], Israeli and Orszag [26]). In meteorological practice the boundary relaxation scheme introduced by Davies [8, 10] is the most frequently used for limited area forecasting using mesoscale modeling. See also the book by Kalnay [29] and the excellent review of McDonald [35]. Typically the forecast equations at the boundary are modified by the addition of a Newtonian relaxation term of the form ∂u ∂u +c = −K(u − u¯) (1) ∂t ∂x where u¯ is prescribed from the host model. It damps differences between the mesoscale and host model at inflow boundaries while mitigating effects of overspecification at the outflow boundaries. Recent research in the atmospheric science community considers transparent boundary conditions for the shallow water equations. McDonald [36],[37] has experimented with transparent boundary conditions for the linearized shallow water equations constructed by incorporating the boundary conditions into equations which describe uni-directional waves. Holstad and Lie [24] and Lie [34] have tested transparent boundary conditions based on characteristics using a mixed finite element formulation of the shallow water equations with encouraging results radiating out the waves for the 2-D nonlinear shallow water equations except for problems with the corners. The quest for an ABC that produces negligible reflections has been, and continues to be, one of the most active areas of research. Most of the popular ABC’s can be grouped into those that are derived from differential equations or those that employ a material absorber. A more recent development is the use of perfectly matched layers (PML) method introduced by Berenger [4]. This approach can be viewed as an improvement on the original idea of sponge layers, since the PML approach allows the solution in the layer to decay for all angles and frequencies. In more detail, one surrounds the computational domain with a finite-thickness layer of specially designed model medium which would either slow down or attenuate all the waves that propagate from inside the domain. The parameters of the layer are chosen such that the wave either never reaches the external boundary, or, if it reaches it, it does reflect back and by the time it reaches the interface between the absorbing (sponge) layer and interior computational domain its amplitude is so small that it will not contaminate the solution. The interface between computational domain and the layer should cause minimal or zero reflection, the latter case being called the PML (see also Clement[5], Karni [31], Kosloff et al. [32], Hu [25], Collino [6], Hayder 4

et al [21], Hayder and Atkins [20]). The work of Hayder et al. [21] is unique in the sense that it applies the PML approach in a nonlinear setting. Dr. R´ emy Baraille drew our attention to a preliminary work [7] to implement PML to the linearized shallow water equations model in oceanography. Abarbanel and Gottlieb [1] provided the general mathematical analysis of the PML method while Abarbanel et al [3] provided a well posed version of PML for advective acoustics. Abarbanel and Gottlieb [2] provided the mathematical framework for use of PML in computational acoustics. The well posedness of PML for linearized Euler equation and for the Cauchy problem is discussed in Rahmouni [44] and Metral [39], respectively. In many applications the PML approach has been shown to provide significantly better accuracy than most other ABC’s. Hayder et al [21] have carried out experiments with the nonlinear version of Euler’s equations. Thus, the PML technique has become the standard to which other techniques must be compared. Because of the quality afforded by the PML scheme, much of the previous literature on ABC’s was made obsolete. This paper is organized as follows. In section 2 we consider the approach of Engquist and Majda [13, 14] for absorbing boundary conditions based on the theory of pseudodifferential operators and we briefly review its application to the linearized shallow-water equations. We also describe the sponge-layer method applied to the linearized 1-D shallowwater equations as a precursor of the perfectly matched layer approach. In section 3 we introduce the PML split-equations approach for the linearized 2-D shallow-water equations. We also explore the implications of the split shallow-water model and weak hyperbolicity of the numerical solution of the PML system thus obtained. An extensive analysis displays the stability problems of this approach. The weak hyperbolicity of the split PML is analyzed for a general case where the shallow water equations include the Coriolis factor. Reflection and transmission of waves at the interface between two perfectly matched layers is also discussed. Numerical implementation of the PML approach for a limited area nonlinear shallow water equations model is presented along with the numerical results using a realistic test. The results obtained point to split PML approach performing satisfactorily for the nonlinear shallow water equations model. To illustrate the performance of the PML method for the nonlinear shallow water equations we compared the accuracy of the computed solution with that arrived at by using only characteristic variables as well as to that arrived at by using the second order aborbing boundary condition of Engquist and Majda. Several numerical experiments were carried out for different values of the mean flow as well as for different parameters of the PML layer. Summary and conclusions are presented in section 4.

2

Engquist-Majda boundary conditions

The approach of Engquist and Majda [13, 14] is based on the theory of pseudodifferential operators (see Taylor [46]). A sequence of local approximate nonreflecting boundary conditions (NRBC) of increasing order is obtained. Consider first the 2-D wave equation in Cartesian coordinates : 1 utt = uxx + uyy c2 5

(2)

Substitute the exponential solution u = Dei(−ωt+k1 x+k2 y) in (2) to obtain the dispersion relation ω2 k 2 = k12 + k22 = 2 (3) c Consider a straight segment of the artificial boundary B with outward normal in the k2 (with |s| ≤ 1) we have from (3) positive x direction. Denoting s = k √ (4) k1 = ±k 1 − s2 on B where ± represent outgoing and incoming waves, respectively. If we wish to obtain an equation on B admitting only outgoing waves, the branch corresponding to + sign is chosen. If we consider (4) as 1-D dispersion relation of an equation P u = 0 on B

(5)

obtained by applying an inverse Fourier transform to (4) which is an exact relation on B. Since k1 = k1 (s) is an irrational function of s, P in (5) is not a differential operator but rather a pseudodifferential operator, which is nonlocal in both time and space. Engquist and Majda [13] approximate the √ nonlocal operator P by a local differential operator E. This is done by approximating 1 − s2 by a rational function (such as Pade approximants). By using rational approximations of increasing accuracy they obtain local boundary conditions Em u = 0 on B of increasing order m. For the case above: 1 ∂ ∂ − )u =0 (6) E1 u = ( ∂x c ∂t x=0 1 ∂2 1 ∂2 1 ∂ 2 E2 u = ( − + )u =0 (7) c ∂x∂t c2 ∂t2 2 ∂y 2 x=0 q

i( ξ 2 − ω 2x + ξt + ωy) , we obtain the In terms of the modified expression (2) u = e symbol of the boundary condition (see Engquist and Majda) as: s

ω2 d − iξ 1 − 2 dx ξ q

2

2

which can be approximated as 1 − ωξ2 = 1 + O( ωξ2 ) at normal incidence ( i.e. ω = 0) 1∂ where iξ corresponds to gives (6) - (7) . c ∂t The next approximation (either Taylor or Pade) to the square root is : s

1−

ω2 1 ω2 ω4 = 1 − ( ) + O( ) ξ2 2 ξ2 ξ4

and multiplying by iξ yields the symbol ∂ 1 iξ + ξ2 − ω2 ∂x 2 or 1 ∂2 1 ∂ 2 1 ∂2 − + )u = 0. ( c ∂x∂t c2 ∂t2 2 ∂y 2 x=0 6

2.1

Application to linear shallow-water equations

The linearized shallow water equations assume the form of a system given by: 











a 0 c b 0 0 0 f 0 ∂ wˆ  ˆ  ˆ   ∂w  ∂w  = 0 a 0  + 0 b c  +  −f 0 0  wˆ ∂t ∂x ∂y c 0 a 0 c b 0 0 0 Physical restrictions on the constants are : c > 0 and 0 < a2 + b2 < c2 . We assume the ∂ wˆ matrices are constant. We diagonalize the normal matrix (multiplying ) by multiplying ∂x by U−1 where the unitary map U is  

U= 

√1 2

0 − √12

0 1 0

√1 2

 

0  

√1 2

and let w = U−1 wˆ to obtain b − √c2 a−c 0 0 ∂w  ∂w  b − √c2 a 0  = 0 +   ∂t ∂x c √ 0 0 a+c 0 2 









0  ∂w  c √  − √f2 +  2  ∂y b 0 0

√f 2

0 √f 2



0  − √f2  w 0

We want to write this equation in the form ∂w ∂w ∂w =A +E + Bw ∂x ∂t ∂y We get:

∂w = ∂x

        

1 a−c

0

0

1 a

0

0

0

1 a+c

0

        

b −a − c  

∂w + ∂t

          



+

           

c √ (a − c) 2

0

− ab

c − √ a 2

c √ a 2 0 0

f √ a 2 0

c √ − b (a + c) (a + c) 2  f √ − 0  (a − c) 2  −

0

−

f √ (a + c) 2

f √ a 2 0

            

∂w ∂y

    w     

Boundary conditions will depend upon whether we have linearized about an inflowing state, i.e. a < 0 or about an outflowing state (with a > 0). For the inflow case Engquist 7

and Majda [13] obtained a sequence of increasing order local nonreflecting boundary conditions 1 0 > Λ1 = a−c   1 0  0 < Λ2 =  a 1 0 a+ c The first absorbing approximation yields w1 w2 and the second approximation

!

=0 x=0

a ∂w2 a ∂w1 +√ − √ f w2 =0 ∂t 2 ∂y c 2 x=0

2.2

Sponge-layer approach

Relatively few studies were carried out for lateral sponge layers (LSL). Davies [8, 10] has formulated LSLs that relax interior flows to the external flow at the boundaries making use of spatially varying sponge coefficients. These coefficients can be either specified empirically or determined optimally in some sense (see Davies [10], and Lehmann [33]) by trying to minimize reflections of outgoing gravity-inertia waves. Following Kar and Turco [30], one can write the one-dimensional shallow-water equations linearized around a basic state as ∂h c ∂u +g = −αu + β h ∂t ∂x H

(8)

∂h ∂u H +H = β u − αh (9) ∂t ∂x c where H is the height of free surface from basic state, h is a deviation of height from H for a perturbed state, u is velocity √ perturbation in x-direction, c is the phase speed of surface gravity waves, i.e. c = gH, and α > 0, β either positive or negative are the sponge layer coefficients. Introducing characteristic variables C + and C − defined by C+ = u + −

C =u−

r

g h H

r

g h H

yields characteristic equations ∂C + ∂C + +c = −(α − β)C + ∂t ∂x ∂C − ∂C − −c = −(α + β)C − ∂t ∂x 8

Using normal mode solution of the form (u, h) = Re

ˆ uˆ(x), h(x) e−iνt

h



i

(10)

ˆ ν being the angular frequency, and substituting (10) in (8) - (9) we obtain uˆ(x), h(x) which implies that we get uˆ(x) = Deik1 x + Eeik2 x r gˆ h(x) = Deik1 x − Eeik2 x H To create a sponge layer (for the non-discretized case) one expands the domain 0 ≤ x ≤ L by an extent d, i.e. for a right boundary L ≤ x ≤ L + d using equations (8) - (9) with |β| ≤ L. We define wave reflectivity as the ratio of complex-valued amplitudes of incoming part of the wave to the outgoing part of wave solution. Using normal mode solutions one obtains at x = L R = e−2αd/c In the discrete case we have ∂ui g c hi+ 21 + hi− 21 + (hi+ 1 − hi− 1 ) = −αi ui + βi 2 2 ∂t ∆ H 2 ∂hi+ 1

2

∂t

+

H H ui + ui+1 (ui+1 − ui ) = β˜i+ 1 −α ˜ i+ 1 hi+ 1 2 c 2 2 ∆ 2

where α, ˜ β˜ are sponge-layer coefficients at n points. Using same techniques, where the width of the sponge layer region L ≤ x ≤ L + d has n grid intervals, one obtains the reflectivity R in functional form as ˜ k) ˜ R = F (n, α ˜ , β, ˜ β∆ with α ˜ = α∆ c , and β = c . Here ∆ is the mesh interval in the x-direction. Numerically, Kar and Turco [30] 0.5c 3. and chose empirically a value α∆ ≥n obtained a satisfactory LSL for n = 6, α ≈ c ∆ The LSL with selective damping works for gravity waves with horizontal phase speed s greater than basic-state zonal velocity. In atmospheric sciences the limited-area model, called also the guest model, uses values supplied to the lateral boundaries by the host model, which usually has a coarser meshwhile the guest is the finer mesh limited-area model. As sharp differences may arise between imposed boundary fields and adjacent mesh points in the limited-area domain a relaxation towards the host model fields is imposed in a zone close to the boundary, i.e. the interior fields are relaxed towards host model fields, in a sponge-layer zone close to the boundary. In the case of the 1-D advection equation the scheme assumes the form: ∂φu ∂u (x, t) + U0 (x, t) = −K(x)(u(x, t) − uE (x, t)) ∂t ∂t 9

(11)

where uE (x, t) is the externally prescribed field by the host and K(x) is a coefficient which is nonzero only in the sponge-layer close to the boundary varying from a large number on the boundary to zero on the interior beyond the sponge layer. For U0 = 0 and K, uE being constants we obtain u(0)e−Kt + uE [1 − e−Kt ]

(12)

i.e. for large K the solution approaches the external field while for K = 0, at the start of the sponge layer, the solution is unchanged. To minimize spurious reflections at the boundary judicious choices of K and the width of the sponge layer have to be made. As shown in McDonald [35] time discretization of (11) leads to ui = (1 − αi )uIi + αi uE i with 0 ≤ αi ≤ 1 , α1 = 0 and αi = 0 for lines beyond the relaxation width or sponge layer zone of n lines in width. Kallberg [28] arrived at a value of αi = 1 − tanh( i−1 ) for 2 a sponge layer with n = 8 lines. Other profiles for α were put forward by Jones, Murphy and Nogeur [27] using a linear i−1 profile of αi = 1 − , with a value of n = 4. McDonald and Haugen [38] proposed a n cosine profile of the form (i − 1)π 1 + cos n αi = (13) 2 In what follows we will see that the sponge layer approach was indeed the precursor of the perfectly matched layer approach (see also Davies [9], and Lehmann[33]).

3

The perfectly matched layer methodology

The perfectly matched layer (PML) technique was first introduced by Berenger [4] for absorbing electromagnetic waves using a finite difference method for solving the Maxwell equations. A sponge-layer or medium is introduced in a region adjacent to the artificial boundary of a computational limited area domain. The novelty of the PML approach is that it is designed such that outgoing waves are absorbed in the layer with theoretically no reflection, i.e. the waves will decay in all directions of propagation in the layers. Work by Hu [25], Abarbanel and Gottlieb [1, 2] Hesthaven [22], Tam et al [45], Hayder et al [21] has extended the technique of PML to the linearized Euler equations and advective acoustics. It is understood that the sponge or buffer-zone solutions themselves need not be physical and they only serve to prevent contamination of the solution in the physical domain of interest by the reflection from the computational boundaries. Usually the absorbing layers are terminated using characteristic boundary conditions. We thus observe an obvious analogy between the PML approach and the sponge boundary conditions used previously in meteorology.

10

In his implementation of the PML, Hu [25] suggested letting the absorption coefficients in the layer vary spatially in the form: σ = σm



d D

λ

where D is the thickness of PML, d is the distance from the interface with the interior domain and λ is a constant. In what follows we will describe a first implementation of PML to the shallow-water (S-W) equations.

3.1

Application of the split PML to the linearized S-W equations

Consider a write-up of the S-W equations as (following Oliger and Sundstrom), [42] d uH + ∇H φ + F = 0 dt d φ + φ∇H · uH = 0 dt where

(14) (15)

∂ d = + uH ∇H dt ∂t φ > u21 + u22

(16) (17)

or in vector from where q = (u1 , u2 , φ)T X ∂ ∂ q+F =0 q+ 2Aj (q) ∂t ∂xj j=1 



(18) 



u2 0 0   A2 =  0 u2 1  0 φ u2

u1 0 1   A1 =  0 u1 0  φ 0 u1

where the eigenvalues of Aj are

uj , uj + c, uj − c, where

1

c = φ2

(19)

(20)

One can simultaneously symmetrize A1 and A2 by the symmetric positive definite transformation matrix 



A linearized version is

φ 0 0  ∗ −1 −1 T R= 0 φ 0  =T 0 0 1 ∂ ′ ∂ ′ X 2Aj (q) q + q + Cq ′ + F ′ = 0 ∂t ∂xj j=1 11

(21)

(22)

The two dimensional linearized S-W equation assumes the form ∂ ∂ ∂ u+U u+ φ=0 ∂t ∂x ∂x

(23)

∂ ∂ ∂ v+U v+ φ=0 ∂t ∂x ∂y

(24)

∂ ∂ ∂u ∂v φ + U φ + c2 ( + )=0 ∂t ∂x ∂x ∂y

(25)

Let us define the following splitting of the fields u, v, φ into, (u1, u2 ), (v1 , v2 ), (φ1 , φ2 ) for the perfectly matched layer: ∂u1 ∂(φ1 + φ2 ) + σx u1 = − ∂t ∂x ∂(u1 + u2 ) ∂u2 + σx u2 = −U ∂t ∂x ∂v1 ∂(φ1 + φ2 ) + σy v1 = − ∂t ∂y ∂(v1 + v2 ) ∂v2 + σx v2 = −U ∂t ∂x ∂φ1 ∂(φ1 + φ2 ) ∂(u1 + u2 ) + σx φ1 = −c2 −U ∂t ∂x ∂x ∂(v1 + v2 ) ∂φ2 + σy φ2 − c2 ∂t ∂y

(26) (27) (28) (29) (30) (31)

In the above the coefficients σx and σy have been introduced for the absorption of waves in the PML layer. We will refer to them as absorption coefficients in this work and they will be assumed to be non negative. We notice that when σx = σy = 0

(32)

we are reduced to the original linearized 2-D shallow-water equations with u = u1 + u2

(33)

v = v1 + v2

(34)

φ = φ1 + φ2

(35)

As such the linearized 2-D S-W equations may be viewed as a special case of the PML equations. The spatial derivatives involve only the total fields of u, v and φ which are assumed to be continuous at the interface between the interior domain and the PML layer. Two types of interface are being created, namely, the interfaces between the interior domain and the PML domain and those between two PML domains.

12

We will proceed to calculate the wave propagation and absorption properties within a perfectly matched layer and attempt to calculate the reflection and transmission coefficients at the interface between two PML domains for the linearized shallow water equations. Let a plane wave in the PML domain be expressed as [u1 , u2 , v1 , v2 , φ1 , φ2 ] = [u10 , u20 , v10 , v20 , φ10 , φ20 ]ei(kx x+ky y−ωt)

(36)

the subscripts being used to denote the amplitude of the components. Substituting (36) into (26)-(31) we obtain (ω + iσx )u10 = kx (φ10 + φ20 ) (37) (ω + iσx )u20 = kx U(u10 + u20 )

(38)

(ω + iσy )v10 = ky (φ10 + φ20 )

(39)

(ω + iσx )v20 = kx U(v10 + v20 )

(40)

(ω + iσx )φ10 = c2 kx (u10 + u20 ) + kx U(φ10 + φ20 )

(41)

(ω + iσy )φ20 = c2 ky (v10 + v20 )

(42)

Let us consider the case when ω − kx U + iσx 6= 0

(43)

Then the components in (36) can be expressed as

φ10

kx U u10 ω − kx U + iσx kx U = v10 ω − kx U + iσx

u20 =

(44)

v20

(45)

ky U c2 kx (ω + iσx ) v10 u10 + = 2 (ω − kx U + iσx ) ω + iσy φ20 =

We also obtain

!

c2 k y ω + iσx v10 ω − kx U + iσx ω + iσy (

kx ω + iσx u10 ) = ω + iσy v10 ky

(46) (47)

(48)

By substituting expressions for φ10 and φ20 into (37) and (39) we obtain u10 = (ω + iσy )v10 =

kx2 c2 (ω + iσy )u10 + c2 kx ky (ω + iσx )v10 (ω − kx U + iσx )2 (ω + iσy )

(49)

c2 kx ky (ω + iσy )(ω + iσx )u10 + c2 ky2 (ω + iσx )2 v10 (ω − kx U + iσx )2 (ω + iσy )

(50)

We should compare this wave for 1-D shallow-water equations with a mean flow U parallel to the x-axis. The S-W equations dispersion relation (neglecting the Coriolis force) is: 13

(ω − Ukx )2 = c2 (kx2 + ky2 )

(51)

φ(x, y, t) = φ0 exp(i(kx x + ky y − ωt))

(52)

where an assumption of wave solutions of the form

is made, i.e., for 1-D shallow water equations, U ± c[1 + (ky2 /ω 2)(U 2 − c2 )]1/2 kx = ω U 2 − c2 (

)

(53)

(See Durran et al., [11]). A first-order Taylor expansion of the square root yields: kx =

cky2 ω ± U ∓c 2ω

(54)

Solving (51) for ky yields, after linear approximation of the square root, ky = ±[

ω − Ukx ckx2 − ] c 2(ω − Ukx )

(55)

For PML we obtain the shallow-water dispersion relationship in the form: (ω − kx U + iσx )2 (ω + iσy )2 − kx2 c2 (ω + iσy )2 − ky2 c2 (ω + iσx )2 = 0

(56)

This reduces, inside the interior domain where σx = σy = 0, to (ω − kx U)2 ω 2 = c2 ω 2(kx2 + ky2 )

(57)

i.e., yielding an identical expression as (51). Collecting u10 terms in (49), dividing by ω + iσy and making use of (48) we get, u10 kx = ±(ω − kx U + iσx ) q 2 c u210 + v10

(58)

Substituting this in (48) we get an expression for ky as ky = ±(ω − kx U + iσy )

ω + iσy v10 q ω + iσx c u2 + v 2 10 10

(59)

If u10 /v10 is real for solution of the S-W equations then one could express u10 and v10 as u10 = A cos φ

(60)

v10 = A sin φ

(61)

where A is a complex number and φ is a real number. If we substitute (60) and (61) into (58) and (59) and we solve for kx and ky , we obtain (by taking only the + sign) kx =

(ω + iσx ) cos φ c + U cos φ 14

(62)

ky =

(ω + iσy ) sin φ c + U cos φ

(63)

Using these expressions we can simplify expressions (44)-(47) so as to get the plane wave solutions to (37)-(42) to be expressed as          

u1 u2 v1 v2 φ1 φ2



         

=

    A    

cos φ U/c cos2 φ sinφ U/c cos φ sin φ c cos2 φ + U cos φ c sin2 φ

         

!

cos φ sin φ iω x+ y−t c + U cos φ c + U cos φ e ×

(64)

σx cos φ σy sin φ (− x− y) c + U cos φ e c + U cos φ

This expression represents a wave propagating with gravity wave speed relative to the mean flow and making an angle φ with the x axis. When both σx and σy 6= 0, the magnitude of the wave decreases exponentially as it propagates in either the x or y directions, respectively.

3.2

Weak hyperbolicity of the split PML

3.2.1

Analysis of the linearized S-W equations

The 2-D linearized shallow water model assumes the form ∂ ∂ ∂ u+U u+ φ − fv = 0 ∂t ∂x ∂x ∂ ∂ ∂ v+U v+ φ + fu = 0 ∂t ∂x ∂y ∂ ∂ ∂u ∂v φ + U φ + c2 ( + )=0 ∂t ∂x ∂x ∂y We can write the system as ∂W ∂W ∂W +A +B + CW = 0 ∂t ∂x ∂y where the vector W is

and the matrices A, B and C are





u    v  1 φ c 



U 0 c   A= 0 U 0  c 0 U 15

(65) (66) (67)





0 0 0   B= 0 0 c  0 c 0





0 −f 0  0  C= f 0  0 0 0

The lower order term CW does not affect the well-posedness and therefore, we study this issue without considering it. Since A and B are symmetric and the equation is hyperbolic, thus the system is strongly well-posed [17]. 3.2.2

Analysis of the split-PML linearized shallow water equations including the Coriolis factor

The 2-D linearized shallow water model assumes the form ∂ ∂ ∂ u+U u+ φ − fv = 0 ∂t ∂x ∂x ∂ ∂ ∂ v+U v+ φ + fu = 0 ∂t ∂x ∂y ∂ ∂u ∂v ∂ φ + U φ + c2 ( + )=0 ∂t ∂x ∂x ∂y

(68) (69) (70)

If we include the absorption coefficients for the PML prior to the splitting we have: ∂ ∂ u+U u+ ∂t ∂x ∂ ∂ v+U v+ ∂t ∂x

∂ φ − f v = −σu u ∂x ∂ φ + f u = −σv v ∂y

∂ ∂ ∂u ∂v φ + U φ + c2 ( + ) = −σφ φ ∂t ∂x ∂x ∂y

16

(71) (72) (73)

The inclusion of the Coriolis factor in the nonlinear shallow-water equations requires the following modification of the PML split form: ∂u1 ∂φ + = −σx u1 ∂t ∂x ∂u ∂u2 +U = −σx u2 ∂t ∂x ∂u3 − fv = 0 ∂t ∂v1 ∂φ + = −σy v1 ∂t ∂y (74) ∂v2 ∂v +U = −σx v2 ∂t ∂x ∂v3 + fu = 0 ∂t ∂φ1 ∂φ ∂u + c2 +U = −σx φ1 ∂t ∂x ∂x ∂v ∂φ2 + c2 = −σy φ2 ∂t ∂y with φ = φ1 + φ2 , u = u1 + u2 + u3 and v = v1 + v2 + v3 . A similar approach was used for the linearized shallow-water equations in oceanography by Darblade et al (1998, [7]). A similar analysis can be applied as in the previously presented case. These equations can be written as ∂W s ∂W s ∂W s + As + Bs + C sW s = 0 ∂t ∂x ∂y where the vector W s is now 

1 1 φ1 φ2 c c

T

0 0 0 0 0 0 c c U U U 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 U U U 0 0 0 0 0 0 0 0 0 0 c c c 0 0 0 U U 0 0 0 0 0 0 0 0



u1 u2 u3 v1 v2 v3

and the matrices As , B s and C s are 

As =

             

17

             



Bs =



             

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

σx 0   0 σx   0 0   0 0  Cs =   0 0   f f   0  0 0 0 The matrix S that diagonalizes As is 

−1 −1 

S=

The inverse is U2   c2 − U 2  

S −1 =

                       

                        

1 U 2c−U

1 c 2c+U

0 0 0 0 0 0 0 c

0 0 0 0 0 0 0 c

0 0 0 0 0 0 0 c

0

0

0

0

1

0

0

0

0

0

0 0 0 0

1 0 0 0

0 −1 1 0

0 0 0 0

0 −1 0 1

0 0 1 0

0

0

0

−1

0

0

0

0

1

0

0 0 0 0 0

1 0 0 0 0

0 −1 0 0 1

1 c 2c+U

              

−

1 U 2c−U

1 c 2c+U 18

             

1

1

U c

0

0



c U

−

0

0

1 U 2c−U

0 0 0 c 0 0 0 0

0 0 0 0 c−U 0 U

U2 c2 − U 2

−

0 0 0 c 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 −f −f −f 0 0 0 σy 0 0 c c 0 0 σx 0 0 0 f 0 0 0 0 0 0 0 0 0 σx 0 0 0 0 0 0 σy

c2 c2 − U 2

0 0 0 0 0

−

0 0 0 0 0 0 0 0

0

0 0 0 −1 0 0 0 1 1 1

0

0

0

0

0

0

0 0 0 0 c+U c 0

                          

Uc Uc − 2 − 2  2 c −U c − U2   

0 0 0 0 0 1 U 2c−U

0 0 1 0 0 1 U 2c−U

1 c 2c+U

1 c 2c+U

                       



S −1 As S =

             

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 U 0 0 0 U −c 0 0 0 U +c

              

Note that S −1 As S is a diagonal matrix with 5 zero eigenvalues being introduced as a result of the splitting. It is straightforward to show that these five additional eigenvalues imply that S and S −1 cannot transform B s into a matrix that can be made symmetric by multiplication with a positive definite diagonal matrix (see Hesthaven [22]). This observation is not conclusive in terms of lack of strong well-posedness. The most general diagonalizer of As is S = T R [3] where the columns of T are the eigenvectors of As and R is a matrix such that the columns of S are the most general representation of the eigenvectors of As .

3.3

Reflection and transmission at an interface between two perfectly matched layers

Let us consider the wave reflection and transmission at an interface between 2 PML domains. This includes the interface between the interior limited area domain and a sponge-like PML domain. The absorbing coefficients σx and σy will be chosen such that σy is the same across the interface normal to y. The linearized S-W equations can be viewed as PML equations with both absorption coefficients being zero across an interface normal to x and y between an interior domain and a PML domain. To illustrate this we present in the following diagram the absorption coefficients on a corner of the computational domain. P ML

(0, σy ) (0, 0) interior

(σx , σy ) (σx , 0) P ML

We want to show that reflection coefficient at an interface downstream normal to x is zero for incident gravity waves, a result that carries over for the interfaces. Let the interface be located downstream at x = 0 and let absorption coefficients be denoted by σx1 and σy on one side of it and σx2 and σy on the other. We want to consider incident, reflected and transmitted waves. a. Incident wave

19

         

u1 u2 v1 v2 φ1 φ2



         

=

    Ai     

cos φi (U/c) cos φi sin φi (U/c) cos φi sin φi c cos2 φi + U cos φi c sin2 φi

         

!

cos φi sin φi iω x+ y−t c + U cos φ c + U cos φ i i e × σy sin φi σx1 cos φi x− y) (− c + U cos φi e c + U cos φi (75)

b. Reflected wave          

u1 u2 v1 v2 φ1 φ2



         

=

    Ai     

− cos φr (U/c) cos φr sin φr −(U/c) cos φr sin φr c cos2 φr − U cos φr c sin2 φr

         

!

− cos φr sin φr iω x+ y−t c − U cos φ c − U cos φ r r e × σy sin φr σx1 cos φr x− y) ( c − U cos φr e c − U cos φr

(76)

c. Transmitted wave          

u1 u2 v1 v2 φ1 φ2

         



=

    Ai     

cos φt (U/c) cos φt sin φt (U/c) cos φt sin φt c cos2 φt + U cos φt c sin2 φt

         

!

sin φt cos φt x+ y−t iω c + U cos φ c + U cos φ t t × e σx2 cos φt σy sin φt (− x− y) c + U cos φt e c + U cos φt (77)

The transmitted wave can have another component of ω − kx U + iσx = 0 in which case It also follows u10 = v10

ω + iσx U = 0 and from (37), (39)

(78)

kx =

(79)

φ10 + φ20 = 0

(80)

u20 = −B sin ψ

(81)

One can express u20 and v20 as

20

v20 = B cos ψ

(82)

For ky one obtains upon dividing (41) by (42) and using (79)-(82) ω + iσy tan ψ U

ky =

(83)

So in this case the plane wave solutions will take form          

u1 u2 v1 v2 φ1 φ2



         

=

        

0 −B sin ψ 0 B cos ψ −Bc2 sin ψ/U Bc2 sin ψ/U



!

(84)

1 tan ψt x+ y − t (− σx2 x − σy tan ψt y) U U U U e

(85)

1 tan ψ x+ − t (− σx x − σy tan ψ y) U U U e U

   iω  e    

Thus the transmitted wave may have a component          

0 −Bt sin ψt 0 Bt cos ψt −Bt c2 sin ψt /U Bt c2 sin ψt /c



   iω  e    

!

At the interface we impose continuity conditions, i.e. that u1 + u2 , v1 + v2 , and φ1 + φ2 be continuous. Since continuity is valid for all values of y at the interface it follows that the coefficients of y in the exponents of (a), (b), (c) (incident, reflected and transmitted values), must be the same. This yields sin φi sin φr = (86) c − U cos φr c + U cos φi

We assume

sin φt sin φi = c + U cos φt c + U cos φi

(87)

tan ψt = tan ψi

(88)

sin φi tan ψi . = U c + U cos ψi

(89)

We then find from the above that −1

φr = 2 tan

1−U φi tan 1+U 2

!

(90)

φt = φi

(91)

ψt = ψi

(92)

21

Using the continuity of u1 + u2 , v1 + v2 and φ1 + φ2 we have, Ai (c + U cos φi ) cos φi − Ar (c − U cos φr ) cos φr − Bi sin ψi = At (c + U cos φt ) cos φt − Bt sin ψt

(93)

Ai (c + U cos φi ) sin φi + Ar (c − U cos φr ) sin φr + Bi cos ψi = At (c + U cos φt ) sin φt + Bt cos ψt

(94)

Ai (c + U cos φi) + Ar (c − U cos φr ) = At (c + U cos φt )

(95)

Using (90)-(92), the above 3 equations may be rewritten as a linear system of 3 homogeneous equations for Ai − At , Ar , and Bi − Bt , provided the coefficient determinant is different from zero. A straightforward but messy calculation shows that the coefficient determinant is not zero for any angle of incidence. Thus the only solution is the trivial, i.e. Ar = 0, At = Ai , and Bt = Bi . Thus for the simple linearized S-W equations if we are at the interface between two PML domains downstream normal to x-axis with PML absorption coefficients (σx1 , σy ), (σx2 , σy ) we get null reflection and transmitted waves will maintain same direction and amplitude as incident waves. This subject to the restrictions imposed by assumptions of linearity. Since it has been shown by Hu [25], Tam [45], and Abarbanel et al [3] that the PML approach is only weakly stable, a low pass filter of the form ! 10 k∆x 1 − sin (96) 2 where k is the wave number and ∆x the grid spacing can be used. If we take gradually varying absorption coefficients of the form σ = σm

d D

!β

(97)

where D is the thickness of the PML domain and d is the distance from interface with interior domain. We can assume that if we apply at the end of the PML domain a radiation boundary condition, then the wave is reflected and upon reflection in the domain its total reflection factor is expressed as 2σm D cos φ − 2 2 2 e β + 1 c − U cos φ

(98)

for gravity waves and for Rossby waves it can be estimated as σm D − e U(β + 1)

(99)

In experiments with linearized Euler equations a gradual σ variation with a value of σm D ≈8 β+1 Usually β ≈ 3, i.e. σm D = 32. 22

(100)

3.4

Numerical implementation considerations

A 2-D nonlinear shallow-water equations solver based on all explicit Miller-Pearce scheme was used [40]. The scheme has a CFL stability condition ∆x c where c is the speed of external gravity waves. Spatial differencing of the nonlinear shallow water equations was carried out on a rectangular domain of 61 × 61 grid points, with a spatial horizontal grid length of ∆x = ∆y = 20km. To demonstrate the implementation of the PML equations we used a number of variable thickness stencils of the PML region - thus allowing a gradual variation in the PML domain. At the edge of the PML domain we applied boundary conditions consisting of a well posed characteristic boundary condition. We compared the results with a control simulation computed on a much larger domain i.e. a domain of [ −240, 240 ] × [ −240, 240 ] which is not affected by the boundary conditions for the integration time span. The interior domain where the unaltered nonlinear shallow-water equations are applied is [ −20, 20 ] × [ −20, 20 ]. To assess the accuracy of the proposed scheme, we have computed either the maximum error in the height field along the line x = 18 inside the inner domain as a function of time. This is done both in the infinity norm as well as in the L2 norm. This allows us to verify the efficiency of the PML absorbing layer as a function of the width of the layer. Characteristic shallow water boundary conditions are used to terminate the PML layer domain. We observe that even for a PML layer of thickness of only 6∆x, the PML scheme outlined above outperforms the characteristic boundary conditions in terms of accuracy, while increasing the width of the PML layer a significant increase in accuracy is obtained. The present results were arrived at without the use of filtering, despite the fact that linear analysis of the scheme points to weak well-posedness ([22]). As expected the height field wave propagates undisturbed from the computational domain with no visible reflections. The important factors to be mentioned are that we solved the PML split version of the nonlinear shallow water equations and the results obtained show that the split scheme including the Coriolis factor is stable and confirms the decaying properties of fields inside the PML layer. The L2 norm difference between the computed and reference solutions along a line inside the inner domain two grid points away from the PML boundary namely at x = 18 (in absolute distance) is presented in Figure 1 as a function of time. This measures the magnitude of the reflected wave at the outflow boundary. PML domains whose thickness are 6, 10, 12, 14, 16, 20 and 30 grid points have been used. The initial conditions are as follows: ∆t ≤

h = hav

(x − x0 )2 + (y − y0 )2 L2 + Ae −

23

(101)

g y − y0 − u = 2 A e f0 L2

(x − x0 )2 + (y − y0 )2 L2

(x − x0 )2 + (y − y0 )2 g x − x0 L2 v = −2 A e f0 L2 −

(102) (103)

where A is the amplitude, (x0 , y0 ) is the initial center of Gaussian vortex, hav is a reference height, f0 being the Coriolis parameter at a reference latitude.

3.5

Results

Here we report on the results of our numerical experiments with various values of mean flow speeds and absorption coefficients for the PML layer. In all experiments we let the center of the Gaussian be at the origin, i.e. (x0 , y0) = (0, 0), the Gaussian amplitude is A = 6 and L = 4. In the first experiment the PML parameters are β = 2 and σm = 0.0012. The mean flow speed varied from Um = 0m/sec to Um = 100m/sec. In figures 1-4 we take Um = 30m/sec. Figure 1 shows the L2 error in the height field between the PML layer domain and the reference solution along a line located inside the interior domain, two grid points away from the PML layer boundary. For comparison we present the results obtained with models having different values of the thickness of the PML layer (6∆x-30∆x) as well as results arrived at using only characteristic variables or results obtained using the second order Engquist boundary conditions. ¿From Figure 1 we note that, during a first stage, the L2 error differences for all the various layer thicknesses start by growing up quickly, while at the second stage, they display a tendency to flatten out, after about 7000 seconds, while at the third stage, they grow up again attaining a size of the order of 10−2 , and become constant thereafter (this is related to the exit of the Gaussian vortex from the model domain). The errors computed using only the characteristic boundary conditions and the second order Engquist boundary conditions are given for comparison. The depth of PML layer obviously impacts on the results, namely the thicker the PML layer is, the better are the results achieved. In the same figure, it can also be seen that all the results for the PML layers of various thicknesses are better than those obtained by either the characteristic boundary condition (highest dashed line) or second order Engquist boundary conditions (second highest line). The time history of the Gaussian pulse transported by the mean flow and propagating towards the right boundary is presented in Figures 2 - 4 in terms of height contours. Figure 2 shows the contour lines initially, while Figure 3 shows the contours when the outer edge of the Gaussian reaches the PML boundary. Figure 4 shows the vortex passing through the PML layer. One notes that the Gaussian vortex passes without reflection through the PML boundary. In the PML region, the contours are deformed, since their amplitudes are reduced towards the outer boundary. In the next experiment we have increased the mean flow speed to Um = 40m/sec. Figure 5 shows the evolution of the L2 error along the same vertical line inside the interior domain, two grid points away from the PML layer boundary. The results are similar to those obtained with Um = 30m/sec, however the times at which the different stages are attained differ, due to faster movement of the Gaussian vortex towards the right boundary. 24

Again when using the characteristic boundary conditions or the second order Engquist boundary conditions we get larger L2 errors. The time history of the Gaussian pulse transported by the mean flow and propagating towards the right boundary is presented in Figures 6 - 8 in terms of height contours. Figure 6 shows the contour lines at the initial time, while Figure 7 shows the contours when the outer edge of the Gaussian reaches the PML boundary. Figure 8 shows the vortex passing through the PML layer. Again one notes that the Gaussian vortex passes without reflection through the PML boundary. In the PML region, the contours are deformed, since their amplitudes are reduced towards the outer boundary. The next experiment takes a high mean flow speed of Um = 100m/sec. Figure 9 show a comparative L2 error norm in the height field. For this high mean flow speed strong oscillations are noticed for the PML layer model after 15000 seconds, when the vortex has moved out of PML region. Again the errors are largest when using either the characteristic boundary conditions or the Engquist boundary conditions. The time history of the Gaussian pulse transported by the high mean flow speed of Um = 100m/sec and propagating towards the right boundary is presented in Figures 10 - 12 in terms of height contours. We also tested the PML method for a non-convective case where Um = 0. Figure 13 represents the evolution of the L2 error. Again the error computed using only characteristic boundary conditions and second order Engquist boundary conditions is given for comparison. In this case the three stages come very early and by 2000 seconds the curves flatten out. As before, the PML performs better than the characteristic and Engquist boundary conditions. See Figures 14 - 16 for the height contours of the Gaussian at different times. Note that since this is a non-convective case the center stays at the origin at all times. We carried out the same numerical experiments with another set of values of the absorption coefficients for the PML layer, i.e. σm and β. Figures 17-28 show the results of an experiment with β = 3 and σm = 0.006 and various values of mean flow speed. The results obtained for Um = 30m/sec display a better separation between the different layer thicknesses of the PML domain. As in all previous cases the PML errors are smaller than those using the characteristic or Engquist boundary conditions. Figures 18 - 20 show the height contours of the Gaussian at different times for a mean flow velocity of Um = 30m/sec. There is no difference between the first two figures 18-19 and the corresponding ones, figures 2-3. The effect of the absorption coefficient manifest itself in Figure 20 as a wider vortex whose right boundary closer to the right boundary of the domain. Figure 21 represents the evolution of the L2 error in the height field for Um = 40m/sec. The error computed using only characteristic boundary conditions or second order Engquist boundary conditions is given for comparison. Again these two didn’t perform as well as the PML schemes. Figures 22 - 24 show the height contours of the Gaussian at different times. As in the previous experiment, we see the same effect of the absorption coefficients on the height contours at the three times, i.e. only at the last time plotted. Figure 25 represents the evolution of the L2 error for Um = 100m/sec. The error computed using only characteristic boundary conditions or second order Engquist boundary conditions are, as in other experiments, the largest. The case with the smallest (6∆x) 25

thickness of PML layer is giving as large errors as the second order Engquist boundary condition when the time is between 5000 to 10000 seconds. This is the only case where the PML layer has a difficulty due to the lower value of σm tested combined with the high value of the mean flow field. See Figures 26 - 28 for the height contours of the Gaussian at different times. A numerical analysis of the propagation of fast gravity waves is given in Appendix A. All the other results obtained show that the PML layer experiments are successful in terms of reduced reflectivity when compared to either charcteristic well-posed boundary conditions or to the second order Engquist boundary conditions for the nonlinear shallowwater equations model. To measure the efficiency of the PML scheme in the inflow layer, we plotted (Figure 29) the L2 error in the height field along the line x = −18 (in the left layer) for the same case described in Figure 1. The maximum error is smaller for the characteristic and Engquist second order boundary conditions. As for the PML, the only difference is at the onset of the third stage when curves flatten out.

4

Summary and conclusions

In this paper we have described and implemented the PML split equations approach for the nonlinear shallow water equations based on an explicit Miller-Pearce scheme finite difference discretization. We have derived the numerical method and developed the code for a rectangular grid domain. The results obtained without using any filters show the robustness of the method. Contrary to results obtained by Hu [25], using a Shapiro 8-th order filter has negligible effects in our case. Extensive numerical testing with an advected Gaussian pulse has shown that the PML boundary conditions radiates out the waves efficiently. Comparison with a set of characteristics well posed boundary conditions as well as the second order Engquist-Majda absorbing boundary conditions show that the PML scheme outperforms the characteristic boundary condition as well as the second order Engquist-Majda absorbing boundary conditions in terms of accuracy, while increasing the width of the PML layer leads to significant increase in accuracy. The research carried out here has a natural extension to the formulation of boundary conditions for advanced mesoscale models, such as the MM5 and the new MRF models, and promises to improve upon the combination of nudging and sponge layer presently used in such models. Issues of computational efficiency combined with numerical experience show that sizable benefits are already obtained with the present split-version of the nonlinear shallow water equations with a PML layer of a width not exceeding 16 cells, a fact which matches practical experience in mesoscale meteorology. Work with PML in the framework of mesoscale models will mean that gravity waves can not only leave the domain but also enter it without hindrance. The fields imposed by the PML must be well balanced - since while the host fields are well-balanced on their own grid - subsequent interpolation to the guest model may upset the balance [37].

26

Our results are encouraging and constitute a step towards using the PML transparent boundary conditions for full 3D atmospheric and ocean models. One avenue to achieve this goal is to implement the PML boundary conditions to a 3D multi-layer shallow water equations model as a way to proceed towards full 3D models. This can be done for the linearized hydrostatic equations by carrying out a normal mode decomposition yielding a shallow water equation for each vertical mode. Development of a non-split version of both the linearized and the nonlinear version of the shallow water equations based on ideas put forward by Abarbanel and Gottlieb [2], Abarbanel et al [3] and Hesthaven [22] is presently also being investigated. Acknowledgements The first author would like to express his gratitude for the support extended to him by the Center of Excellence (COE) at Florida State University during part of the write-up of this paper. This is a contribution from the climate institute, a center of excellence funded by the FSU research foundation.He also acknowledges support from NSF grant ATM97B1472 managed by Dr. Pamela Stephens. He would also like to express his gratitude for the support and encouragement provided to him during his short stay at NCAR, fall 1999 made possible by DR Ying-Hwa Kuo, Head Mesoscale Prediction Group and by Dr Joseph Klemp, Head Mesoscale Dynamics Group both at NCAR/MMM. The secretarial support and the expert librarian support is also acknowledged.

References [1] Abarbanel, S. and D. Gottlieb, 1997: A mathematical analysis of the PML method, J. Comput. Phys., 134, 357-363. [2] Abarbanel, S. and D. Gottlieb, 1998: On the construction and analysis of absorbing layers in CEM, Appl. Numer. Math., 27, 331-340. [3] Abarbanel, S., Gottlieb, D. and J. S. Hesthaven, 1999 : Well-posed perfectly matched layers for advective acoustics, J. Comput. Phys., 154, 266-283. [4] Berenger,J. -P., 1994:A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys., 127, 185-200. [5] Clement, A., 1996: Coupling of two absorbing boundary condition for 2D timedomain simulation of free surface gravity waves, J. Comput. Phys., 126, 139-159. [6] Collino, F., 1996: Perfectly matched absorbing layers for paralaxial equations, J. Comput. Phys., 131, 164-180. [7] Darblade, G., Baraille, R., le Roux, A.-Y., Carton, X. and D. Pinchon, 1997: Conditions limites non r´ efl´ echissantes pour un mod` ele de Saint-Venant bidimensionnel barotrope lin´ earis´ e, C. R. Acad. Sci. Paris, 324, 485-490. [8] Davies, H. C., 1976 : A lateral boundary formulation for multilevel prediction models, Quart. J. Roy. Meteor. Soc., 102, 405-418. 27

[9] Davies, H. C., 1983 : Techniques for limited area modeling. Seminar 1983: numerical methods for weather prediction. ECMWF, 5–9 September 1983,213-236. [10] Davies, H. C., 1985 : Limitations of some common lateral boundary schemes used in regional NWP models, Mon. Wea. Rev., 111, 1002-1012. [11] Durran, D. R., Yang, M. J., Sinn, D. N. and R. G. Brown, 1993: Toward more accurate wave permeable boundary conditions, Mon. Wea. Rev., 121, 604-621. [12] Elvius, T., and A. Sundstrom, 1973: Computationally efficient schemes and boundary conditions for a fine mesh barotropic model based on the shallow-water equations, Tellus, 25, 132-156. [13] Engquist, B. and A. Majda, 1977 : Absorbing boundary conditions for the numerical simulation of waves, Math. Comp., 31, no 139, 629-651. [14] Engquist, B. and A. Majda, 1979 : Radiation boundary conditions for acoustic and elastic wave calculations, Comm. Pure. Appl. Math., 32, 313-357. [15] Givoli, D., 1991: Nonreflecting boundary conditions, J. Comput. Phys., 94, 1-29. [16] Goodrich, J.W. and T. Hagstrom, 1997 : A comparison of two accurate boundary treatments for computational aeroacoustics. AIAA paper 97-1585.11pp. [17] Gustafson, B., and A. Sundstrom, 1978: Incompletely parabolic problems in fluid dynamics, SIAM J. Appl. Math., 35, 343-357. [18] Haltiner, G. J. and R. T. Williams, 1980: Numerical Prediction and Dynamic Meteorology, John Wiley & Sons, New York. [19] Harari, I., Slavutin, M. and E. Turkel, 2000: Analytical and numerical studies of a finite element PML for the Helmholtz equation, J. Comp. Acoustics, 8, 121-137. [20] Hayder, M. E. and H. L. Atkins, 1997: Experience with PML boundary conditions in fluid flow computation in T.L. Gears Ed.. Collection of Abstract of Symposium on computational methods for unbounded domains. Univ. of Colorado at Boulder, July 27-31 (Kluwer Academic Press) to appear. [21] Hayder, H. E., Hu, F. Q. and M. Y. Hussaini, 1999 : Towards perfectly absorbing boundary conditions for Euler equations, AIAA Journal, 37, no 8 , 912-918. [22] Hesthaven, J. S., 1998 : On the analysis and construction of perfectly matched layers for the linearized Euler equations, J. Comput. Phys., 142, 129-147. [23] Higdon, R. L., 1986: Absorbing boundary conditions for difference approximations to the multi-dimensional wave equation. Mathematics of Computation ,Vol 47, No. 176,437-459. [24] Holstad, A. and Lie, I., 2001: Transparent boundary conditions using a mixed finite element formulation of the shallow water equations, Norwegian Met. Inst. (DNMI), Research Report No. 120, 48 pp. 28

[25] Hu, F. Q., 1996 : On absorbing boundary conditions for linearized Euler equations by a perfectly matched layer, J. Comput. Phys., 129, 201-216. [26] Israeli, H. and S. A. Orszag, 1981 : Approximation of radiation boundary conditions, J. Comput. Phys., 41, 115-135. [27] Jones, R.G., J.M. Murphy and M. Noguer, 1995: Simulation of a climate change over Europe using a nested regional climate model. I: Assessment of control climate including sensitivity to location of lateral boundary conditions, Q. J. Roy. Met. Soc,121, 1413-1499 [28] Kallberg, Per, 1977: Test of a boundary relaxation scheme in a barotropic model, ECMWF Res. Dept Internal Report No.3, 21pp. [29] Kalnay, E., 1999: Numerical Weather Forecasting and Predictability, Cambridge University Press, Cambridge, MA. [30] Kar, K. S. and R. P. Turco, 1995 : Formulation of a lateral sponge layer for limitedarea shallow-water models and an extension for the vertically stratified case, Mon. Wea. Rev., 123, 1542-1559. [31] Karni, S., 1996: Far field filtering operators for suppression of reflections from artificial boundaries, SIAM J. Numerical Anal, 33, 1014-1047. [32] Kosloff, R. and D. Kosloff, 1986: Absorbing boundaries for wave propagation problems, J. Comput. Phys., 63, 363-376. [33] Lehmann, R., 1993: On the choice of relaxation coefficients for Davis lateral boundary scheme for regional weather prediction models, Meteor. Atmos. Phys, 53, 1-14. [34] Lie, I., 2001 : Well-posed transparent boundary conditions for the shallow water equations, Appl. Numer. Math.. To appear. [35] McDonald, A., 1997: Lateral boundary conditions for operational regional forecast models: A review HIRLAM Technical Report No 32, 31pp. [36] McDonald, A., 2001: A step toward transparent boundary conditions for meteorological models, Irish Met. Service, Technical Note No. 57, 22 pp. [37] McDonald, A., 2001: Well posed boundary conditions for semi-Lagrangian schemes: The two dimensional case, HIRLAM Technical Report No. 47, 38 pp. [38] McDonald, A. and J.E. Haugen, 1992: A two-time level three dimensional, semilagrangian and semi-implicit grid point model. Mon. Wea. Review,120, 2603-2621 [39] Metral, J. and O. Vacus, 1999: Well-posedness of the Cauchy problem associated with Berenger’s system, Comp. Rendus de L’Acad. des Sci. Serie I-Math., 328, 847-852. [40] Miller, M. J. and R. P. Pearce, 1974: A three dimensional primitive equation model of cumulonimbus convection, Quart, J. Roy. Met. Soc, 100, 133-154. 29

[41] Miller, M. J. and A. J. Thorpe, 1981: Radiation conditions for the lateral boundaries of limited area models.Quart, J. Roy. Met. Soc., 107, 615-628. [42] Oliger, J. and Sundstrom A., 1978 : Theoretical and practical aspects of some initial boundary value problems in fluid dynamics, SIAM J. Appl. Math., 35, no 3 , 419-447. [43] Perkey, D. J. and C. W. Kreitzberg, 1976 : A time-dependent lateral boundary scheme for limited area primitive equations models, Mon. Wea. Rev., 104, 744-755. [44] Rahmouni, A., 2000: A well-posed unsplit PML model for linearized Euler equations, Comp. Rendus de L’Acad. des Sci. Serie I-Math., 331, 159-164. [45] Tam, C. K. W., Aurialt, L. and F. Cambuli, 1998 : Perfectly matched layer as an absorbing boundary condition for the linearized Euler equations in open and ducted domains, J. Comput. Physics, 144, 213-244. [46] Taylor ,M., 1981: Pseudo-differential operators, Princeton University Press. [47] Tsynkov, V. S., 1998 : Numerical solution of problems on unbounded domains. A review, Applied Num. Math., 27, 465-532. [48] Turkel, E., 1983: Progress in Computational Physics, Computers and Fluids, 11, 121-144. [49] Turkel, E. and A. Yefet, 1998: Absorbing PML boundary layers for wave-like equations, Appl. Numer. Math., 27, 533-557.

30

A

Appendix

Propagation of fast moving gravity waves for the Gaussian = advection equation. For the sake of simplicity, let us consider the advection equation with c > 0: ut + cux = 3D0 where=20

for t > 0

u = 3De− ln 2( L∆x ) x

2

We discretize the equation using a simple centered = difference : dul c = 3D (−ul+1 + ul−1 ) dt 2∆x where l stands for the grid number. Using=20

(A.1) (A.2)

(A.3)

uj = 3Dei(ωt−jα∆x) ,

(A.4)

c sin α∆x ∆x

(A.5)

where=20 ω = 3D and where

α∆x ∈ [−π, π]

(A.6)

We find that the group velocity of the dispersion relation is = given by:=20 ∂ω = 3Dc cos(α∆x) (A.7) ∂α Let us now consider the range π 0 ≤ α∆x ≤ (A.8) 2 which implies cos(α∆x) > 0. This means that all wave number components of the iniπ travel to the right, i.e., in the positive tial Gaussian = which satisfy 0 < α∆x ≤ 2 x−direction. However for=20 π ≤ α∆x ≤ π, = 20 cos(α∆x) < 0 (A.9) 2 ∂ω < 0. These wave = numbers of the initial Gaussian travel to the left, i.e. in the i.e., ∂α = negative=20 x−direction. The fastest wave propagating to the right is given by=20 ∂ω = 3Dc cos(0) = 3Dc (A.10) ∂α for α∆x = 3D0 which is the long-wave number limit, whereas the fastest wave=20 propagation speed to the left is given by: ∂ω = 3D − c (A.11) ∂α for α∆x = 3Dπ that is all the grid to grid waves travel at speed of −c. Use of either a finer mesh or a higher order differencing method will yield a more accurate solution, with less short-wave numbers=20 travelling to the left. 31

Captions Fig. 1 The evolution of the L2 height error along the line x = 18 for various PML layer thicknesses for Um = 30m/sec. The error computing only characteristic or second order Engquist-Majda boundary conditions is given for comparison. Fig. 2 Height contours for the Gaussian profile propagating with convective mean velocity of Um = 30m/sec. The computed result is given at t = 90sec. for a PML layer thickness of 6 grid points. Fig. 3 Same as Fig. 2 but the computed result is at t = 720sec. Fig. 4 Same as Fig. 2 but the computed result is at t = 16200sec. Fig. 5 The evolution of the L2 height error along the line x = 18 for various PML layer thicknesses for Um = 40m/sec. The error computing only characteristic or second order E-M boundary conditions is given for comparison. Fig. 6 Same as Fig. 2 but for Um = 40m/sec. Fig. 7 Same as Fig. 6 but the computed result is given at t = 5400sec. Fig. 8 Same as Fig. 6 but the computed result is given at t = 9900sec. Fig. 9 Same as Fig. 1 but for mean velocity Um = 100m/sec. Fig. 10 Same as Fig. 2 but for Um = 100m/sec. Fig. 11 Same as Fig. 10 but the computed result is given at t = 2070sec. Fig. 12 Same as Fig. 10 but the computed result is given at t = 5400sec. Fig. 13 Same as Fig. 1 but for the non-convective case when Um = 0m/sec. Fig. 14 Same as Fig. 2 but for the non-convective case when Um = 0m/sec. and result is given at t = 90sec. Fig. 15 Same as Fig. 14 but the computed result is given at time t = 2700sec. Fig. 16 Same as Fig. 14 but the computed result is given at t = 3600sec. Fig. 17 The evolution of the L2 height error along the line x = 18 for various PML layer thicknesses for Um = 30m/sec but for PML layer parameters σN = 0.006 and λ = 3. The error computing only characteristic or second order Engquist-Majda boundary conditions is given for comparison. Fig. 18 Same as Fig. 2 but for σN = 0.006 and λ = 3 and Um = 30m/sec. Fig. 19 Same as Fig. 18 but the computed result is given at time t = 7200sec. Fig. 20 Same as Fig. 18 but the computed result is given at t = 16200sec. 32

Fig. 21 Same as Fig. 17 but for Um = 40m/sec. Fig. 22 Same as Fig. 18 but for Um = 40m/sec. and result is given at t = 90sec. Fig. 23 Same as Fig. 22 but the computed result is given at time t = 5400sec. Fig. 24 Same as Fig. 22 but the computed result is given at t = 9900sec. Fig. 25 Same as Fig. 17 but for Um = 100m/sec. Fig. 26 Same as Fig. 22 but for Um = 100m/sec. and result is given at t = 90sec. Fig. 27 Same as Fig. 26 but the computed result is given at time t = 2070sec. Fig. 28 Same as Fig. 26 but the computed result is given at t = 5400sec. Fig. 29 Same as Fig. 1, but for the inflow layer case x = −18.

33

L1=6,A=4,Um=30,dt=9,sigmaN=0.012,lambda=2,no shapiro,dx=dy=20km

1

10

0

10

−1

10

−2

||h−href||2

10

−3

10

−4

10

CHAR. BC06 2nd BC06 PML06 PML10 PML16 PML20 PML30

−5

10

−6

10

−7

10

0

0.5

1

1.5 t(Sec.)

2

2.5 4

x 10

Figure 1: The evolution of the L2 height error along the line x = 18 for various PML layer thicknesses for Um = 30m/sec. The error computing only characteristic or second order Engquist-Majda boundary conditions is given for comparison.

34

Figure 2: Height contours for the Gaussian profile propagating with convective mean velocity of Um = 30m/sec. The computed result is given at t = 90sec. for a PML layer thickness of 6 grid points.

35

Figure 3: Same as Fig. 2 but the computed result is at t = 720sec.

36

Figure 4: Same as Fig. 2 but the computed result is at t = 16200sec.

37

L1=6,A=4,Um=40,dt=9,sigmaN=0.012,lambda=2,no sharpiro,dx=dy=20km

1

10

0

10

−1

10

−2

||h−href||2

10

−3

10

−4

10

CHAR. BC06 2nd BC06 PML06 PML10 PML16 PML20 PML30

−5

10

−6

10

−7

10

0

0.5

1

1.5

2

2.5

t(Sec.)

Figure 5: The evolution of the L2 height error along the line x = 18 for various PML layer thicknesses for Um = 40m/sec. The error computing only characteristic or second order E-M boundary conditions is given for comparison.

38

4

x 10

Figure 6: Same as Fig. 2 but for Um = 40m/sec.

39

Figure 7: Same as Fig. 6 but the computed result is given at t = 5400sec.

40

Figure 8: Same as Fig. 6 but the computed result is given at t = 9900sec.

41

L1=6,A=4,Um=100,dt=9,sigmaN=0.012,lambda=2,no sharpiro,dx=dy=20km

1

10

0

10

−1

10

−2

||h−href||2

10

−3

10

−4

10

CHAR. BC06 2nd BC06 PML06 PML10 PML16 PML20 PML30

−5

10

−6

10

−7

10

0

0.5

1

1.5 t(Sec.)

2

Figure 9: Same as Fig. 1 but for mean velocity Um = 100m/sec.

42

2.5

3 4

x 10

Figure 10: Same as Fig. 2 but for Um = 100m/sec.

43

Figure 11: Same as Fig. 10 but the computed result is given at t = 2070sec.

44

Figure 12: Same as Fig. 10 but the computed result is given at t = 5400sec.

45

L1=6,A=4,Um=0,dt=9,sigmaN=0.012,lambda=2,no shapiro,dx=dy=20km

−1

10

−2

10

−3

||h−href||2

10

−4

10

−5

10

CHAR. BC06 2nd BC06 PML06 PML10 PML16 PML20 PML30

−6

10

−7

10

0

5000

10000 t(Sec.)

Figure 13: Same as Fig. 1 but for the non-convective case when Um = 0m/sec.

46

15000

Figure 14: Same as Fig. 2 but for the non-convective case when Um = 0m/sec. and result is given at t = 90sec.

47

Figure 15: Same as Fig. 14 but the computed result is given at time t = 2700sec.

48

Figure 16: Same as Fig. 14 but the computed result is given at t = 3600sec.

49

L1=6,A=4,Um=30,dt=9,sigmaN=0.006,lambda=3,no sharpiro,dx=dy=20km

1

10

0

10

−1

10

−2

||h−href||2

10

−3

10

−4

10

CHAR. BC06 2nd BC06 PML06 PML10 PML16 PML20 PML30

−5

10

−6

10

−7

10

0

0.5

1

1.5

2

2.5

t(Sec.)

Figure 17: The evolution of the L2 height error along the line x = 18 for various PML layer thicknesses for Um = 30m/sec but for PML layer parameters σN = 0.006 and λ = 3. The error computing only characteristic or second order Engquist-Majda boundary conditions is given for comparison.

50

4

x 10

Figure 18: Same as Fig. 2 but for σN = 0.006 and λ = 3 and Um = 30m/sec.

51

Figure 19: Same as Fig. 18 but the computed result is given at time t = 7200sec.

52

Figure 20: Same as Fig. 18 but the computed result is given at t = 16200sec.

53

L1=6,A=4,Um=40,dt=9,sigmaN=0.006,lambda=3,no shapiro,dx=dy=20km

1

10

0

10

−1

10

−2

||h−href||2

10

−3

10

−4

10

CHAR. BC06 2nd BC06 PML06 PML10 PML16 PML20 PML30

−5

10

−6

10

−7

10

0

0.5

1

1.5 t(Sec.)

Figure 21: Same as Fig. 17 but for Um = 40m/sec.

54

2

2.5 4

x 10

Figure 22: Same as Fig. 18 but for Um = 40m/sec. and result is given at t = 90sec.

55

Figure 23: Same as Fig. 22 but the computed result is given at time t = 5400sec.

56

Figure 24: Same as Fig. 22 but the computed result is given at t = 9900sec.

57

L1=6,A=4,Um=100,dt=9,sigmaN=0.006,lambda=3,no shapiro,dx=dy=20km

1

10

0

10

−1

10

−2

||h−href||2

10

−3

10

−4

10

CHAR. BC06 2nd BC06 PML06 PML10 PML16 PML20 PML30

−5

10

−6

10

−7

10

0

0.5

1

1.5 t(Sec.)

Figure 25: Same as Fig. 17 but for Um = 100m/sec.

58

2

2.5 4

x 10

Figure 26: Same as Fig. 22 but for Um = 100m/sec. and result is given at t = 90sec.

59

Figure 27: Same as Fig. 26 but the computed result is given at time t = 2070sec.

60

Figure 28: Same as Fig. 26 but the computed result is given at t = 5400sec.

61

L1=6,A=4,Um=30,dt=9,sigmaN=0.012,lambda=2,no shapiro,x=−18,dx=dy=20km

1

10

0

10

−1

10

−2

||h−href||2

10

−3

10

−4

10

CHAR. BC06 2nd BC06 PML06 PML10 PML16 PML20 PML30

−5

10

−6

10

−7

10

0

0.5

1

1.5

2

t(Sec.)

Figure 29: Same as Fig. 1, but for the inflow layer case x = −18.

62

2.5 4

x 10