## Abstract

Atmospheric motions are generally characterized by a wide range of multiple length and time scales, and a numerical method must use a fine grid to resolve such a wide range of scales. Furthermore, a very fine grid requires an extremely small time step in order to keep explicit time integration schemes stable. Therefore, high-resolution meteorological simulations are very expensive.

A novel multiscale modeling approach is, therefore, presented for simulating atmospheric flows. In this approach, a prognostic variable representing a highly intermittent multiscale feature is decomposed into a significant and a nonsignificant part using wavelets, where the significant part is represented by a small fraction of the wavelet modes. The proposed multiscale methodology has been verified for simulating three cases: Smolarkiewicz’s deformational flow model, warm thermals in a dry atmosphere, and the dynamics of a vortex pair with ambient stable stratification. Comparisons with benchmark simulations and with a reference model are evidence for the convergence and stability of the proposed model. The comparison with the reference model has revealed that about 93% of the grid points are not necessary to resolve the significant motion in a warm thermal simulation, saving about 96% of the CPU time. Moreover, the CPU time varies linearly with the number of significant wavelet modes, showing that the present fully implicit adaptive model is asymptotically optimal for this simulation. These primary results point toward the benefit of constructing multiscale atmospheric models using the adaptive wavelet methodology.

## 1. Introduction

A quantitative understanding of the atmosphere requires an accurate representation of various multiscale meteorological processes such as cold fronts, rising thermals, moist convection, breaking of gravity waves, and turbulence, which has long been a challenge for numerical atmospheric models (e.g., see Carpenter et al. 1990). Recently, the atmospheric modeling community is exploring several approaches to address the scale gap between the large-scale physics, which can be resolved explicitly from first principles, and the smaller-scale processes that need to be parameterized. These approaches include the use of massively parallel high-performance computing (HPC) methods, adaptive mesh refinement (AMR) techniques, semi-Lagrangian advection schemes, Lagrangian stochastic models, and large-eddy simulation (LES) methods to name a few (e.g., Moeng 1984; Skamarock et al. 1989; Wilson and Sawford 1995; Behrens 1996; Rosenberg et al. 2006; Wehner et al. 2008). As a result, the resolution of current state-of-the-art operational weather prediction and climate models has significantly increased thanks to the latest HPC resources (e.g., see chapter 9 in Hamilton and Ohfuchi 2007). However, owing to the enormous range of length and time scales, one rarely has the luxury of using HPC resources to explicitly resolve many important phenomena (e.g., the mesoscale circulation that organizes moist convection). Therefore, much effort has also been paid on parameterizing the effects of unresolved physics. The research in two methodologies—the increase of the spatial resolution with more powerful HPC algorithms and the search for the best parameterization method—has been a current goal of the atmospheric modeling research community (Behrens 1998, 2006; Jablonowski et al. 2006; Bacon et al. 2007). The present paper aims to explore a novel multiscale technique for high-resolution atmospheric modeling, which uses an optimal number of adaptive grid points to capture isolated multiscale features.

Indeed, there are several common techniques in this direction, including nested grids, stretched grids, and AMR methods. Recently, there has been growing interest in the AMR approach for simulating isolated multiscale features in the atmosphere. Skamarock et al. (1989) presented an adaptive atmospheric model using the AMR method of Berger and Colella (1989), where a Richardson method was used for a posteriori estimation of truncation error. Tomlin et al. (1997) studied how an adaptive grid chemical transport model can reveal new features of plume concentration profiles. Blayo and Debreu (1999) verified that this approach gains in CPU time by about a factor of 3 in the context of simulating a barotropic modon. Behrens (1998) proposed a dynamically adaptive shallow-water model in the Cartesian *x–y* plane. Bacon et al. (2000) and Boybeyi et al. (2001) presented a nonhydrostatic model, the Operational Multiscale Environment Model with Grid Adaptivity (OMEGA), using unstructured, triangulated grids that can be dynamically or statically adapted to localized features of interest. OMEGA has further been verified for multiscale simulation of hurricanes (e.g., Bacon et al. 2007). Jablonowski et al. (2006) verified that a spherical 2D AMR model detects the chosen features reliably, and helps preserve the shape and amplitude of the transported field while saving computational efforts. Lauter et al. (2007) presented a barotropic AMR model of the atmosphere. Behrens (2006) discussed in more technical detail the AMR approach and parameterization in the context of an AMR model. For a review on the development of classic AMR models for global atmospheric chemistry/transport or operational cyclone forecasting, see the work of Nikiforakis (2005) and Bacon et al. (2003). Most classic AMR algorithms are designed to reduce local errors estimated according to a Richardson criterion, where a coarse-grid solution is compared with a fine-grid solution.

In contrast, a principal objective of the present study aims to explore a novel multiscale approach, the adaptive wavelet collocation methodology (AWCM), which that represents the most energy-containing, intermittent motion using a small fraction of the wavelet modes, where the residual motion corresponds to a large proportion of them. Hence, the AWCM tracks the multiscale energy contributions associated with intermittent features. In contrast to the classic AMR method, where an adaptive grid aims to capture localized features of a flow, the AWCM represents localized features in terms of multiscale wavelet modes, where a flow is decomposed into a significant (i.e., energy containing) and a nonsignificant (i.e., residual) part. The significant part is calculated on an adaptive grid that is generated naturally as a result of discarding the nonsignificant part. In other words, the nonsignificant part represents numerical error if the significant part contains all wavelet modes, resolving a certain desired dynamics.

Despite the fact that AWCM has advantages over the classic AMR approach, there are some challenges that need to be addressed for extending the AWCM to atmospheric applications. First, a locally refined mesh of the AWCM will also constrain the time step globally for an explicit scheme satisfying the CFL criterion. To address this disadvantage, Alam et al. (2006) proposed to solve the governing equations in the simultaneous space–time domain as if the time is a spacelike variable. Clearly, this approach is not practical for atmospheric flows. The present work has proposed to use an implicit scheme for the advection terms as well as for other terms. However, solving a large nonlinear system at each time step requires the construction of the Jacobian matrix of the discretized nonlinear system, which will also be multiplied with the error vector *υ*, and hence the computational overhead would grow like , where is the number of grid points. To address this overhead, the use of a multilevel solver has been proposed. Second, the AWCM provides an adaptive methodology only, which does not eliminate the need for subgrid-scale parameterization. As discussed, for example, in chapter 9.2 in Behrens (2006), the parameterization techniques must also be redesigned in order to take full benefit of adaptive methodologies. The present study has addressed the first of the above two.

In section 2, the proposed multiscale modeling approach is presented. The convergence of the proposed model is verified in section 3. Numerical simulations of warm thermals and that of vortex-pair dynamics are presented in sections 4 and 5, respectively. The model is applied to simulate a flow over an isolated topography in section 6. Finally, the proposed development has been summarized in section 7 with a brief discussion on further research in this direction.

## 2. A wavelet-based atmospheric model

### a. Background

Over the past 25 years, wavelet theory generated a tremendous amount of interest in many areas of research, ranging from applied mathematics (e.g., Grossmann and Morlet 1984; Urban 2002; Cohen 2003) to ocean modeling (e.g., Jameson and Miyama 2000). The introduction of the second-generation wavelet theory by Sweldens (1995) accelerated in the development of sophisticated computational fluid dynamics (CFD) techniques for simulating complex geometry flows (Kevlahan and Vasilyev 2005). Depending on specific applications, wavelet methods can use a Galerkin approach, where an average residual error is minimized in the entire domain, or can use a collocation approach, where the residual error vanishes on a set of collocation points. A comprehensive review on recent developments of wavelet methods for CFD applications can be found in Schneider and Vasilyev (2010), where both the Gelarkin and the collocation-based methods are also reviewed. The present development has been originated from the second-generation AWCM for time-dependent (e.g., Vasilyev and Bowman 2000) and elliptic (e.g., Vasilyev and Kevlahan 2005) problems, and is an extension of the wavelet methodology, studied by Alam (2006), to the field of atmospheric modeling. This AWCM uses a near-optimal degree of freedom to simulate, for example, homogeneous isotropic turbulent flows with appropriate subgrid-scale (SGS) parameterization (e.g., De Stefano and Vasilyev 2010). To adapt the spatial and temporal step sizes to the locally varying solution of a turbulent flow, Alam (2006) proposed a simultaneous space–time AWCM for solving nonlinear partial differential equations (PDEs; e.g., Alam et al. 2006). Kevlahan et al. 2007) compared the solution of the vorticity equation between the space–time AWCM and a Fourier collocation method, and found that only about 3% of the space time modes are sufficient to compute a flow with high temporal intermittency. Mehra and Kevlahan (2008) used spherical wavelets to extend the AWCM for solving PDEs on the sphere.

In summary, the AWCM (e.g., Vasilyev and Bowman 2000; Alam 2006) as well as other wavelet-based methods (see Schneider and Vasilyev 2010) provide an efficient framework for solving PDEs or simulating fluid flows, where the solution has important small-scale features only on a fraction of the domain. In the following section, the AWCM that is used in the present work is described briefly.

### b. The adaptive wavelet collocation methodology

For a brief introduction of the AWCM, let us consider a multiscale decomposition of any prognostic variable:

where and are sets representing values of the index *k*, is a set of coefficients for the wavelet basis , *j*_{min}, *l*, *μ* are integer numbers, *d* stands for dimension, and **x** = (*x*, *y*, *z*). Since this is a fixed time representation, *t* is omitted for simplicity. For the present work, we have *d* = 2 and **x** = (*x*, *z*). The decomposition of (1) has a significant part

at level *j*, which is given by wavelets such that , and a nonsignificant or residual part that is defined by

where *ε* is a tolerance that dictates the maximum wavelet level *j* and the symbol ‖*u*‖* _{p}* represents a certain measure (e.g., rms value) of

*u*. Note, the coefficients are scaled with a threshold

*ε*‖

*u*‖

_{p}. In fact, the significant part is a nonlinear approximation to

*u*, which does not oscillate at a frequency larger than 2

*(e.g., Cohen et al. 2001). The maximum error for a piecewise smooth function*

^{j}*u*(

**x**) can be estimated as

if wavelets with *n*-vanishing moments are used (Donoho 1992; Vasilyev and Kevlahan 2005; Alam 2006). Clearly, a sufficiently small *ε* implies a large and the residual part can be treated as a numerical error if is sufficiently large. For all simulations in this paper, we have used *n* = 6. A rigorous mathematical analysis of wavelet theory, which is not the subject of this article, can be found in standard text books (e.g., Mallat 1999; Cohen 2003). Primary questions for extending the AWCM to atmospheric modeling include the method of (i) calculating coefficients [(section 2c(1)], (ii) approximating spatial derivatives [section 2c(2)], and (iii) time evolution for prognostic variables (section 2g). The following development uses the lifting scheme (e.g., Sweldens 1995) and interpolating wavelets (e.g., Donoho 1992; Vasilyev and Bowman 2000).

### c. Calculation of wavelet coefficients and numerical differentiation

#### 1) Calculation of wavelet coefficients

Let represent a given *u*(**x**) on the *k*th grid point at level *j* [i.e., ]. First, the dataset is grouped into two: with those having an even index and with those having an odd index: . For example, , , and . Denoting the *j*th level as high resolution, a low-resolution data is defined by . Second, the missing details between high- and low-resolution information (i.e., between and ) are predicted using polynomials.

Consider polynomials of degree *n* based on each and its *n* neighbors, predict the value by , and define . The symbol means that a polynomial is based on and its *n* neighbors, and is evaluated on the corresponding grid point. Third, low-resolution data s are updated by , which ensures that the corresponding wavelets will have a zero mean (Sweldens 1995). This completes a two-scale decomposition of (the ‘′’ is dropped from *k* for simplicity), and is known as the second-generation interpolating wavelet transform.

A forward wavelet transform can be described as

and a multiscale decomposition is obtained by applying the above wavelet transform recursively. The original data can be recovered by the following inverse wavelet transform:

at any level *j*.

A two-level decomposition of sin(2*πx*) is presented in Fig. 1a, showing that sin(2*πx*) is well approximated locally if the grid point at *x* = 0.5 is discarded, where is small. Furthermore, large errors occur near *x* = 0.2 and *x* = 0.8, which can be minimized by adding new grid points.

#### 2) Approximation of derivatives

A spatial derivative is calculated using a collocation method that is based on the multiscale decomposition of the derivative, and is computed from coefficients of (1) by differentiating polynomials analytically (e.g., see Alam 2006). Let be a differential operator, represent the coefficients in (1), and **c** represents the function evaluation on the corresponding grid. Then, a weighted residual collocation method implies that the residual vanishes [i.e., ], at each grid point . For example, knowing the functional form of the wavelet basis, can be obtained by differentiating (2). Then, following Alam (2006), we have , where is the resulting differentiation matrix, , and is the wavelet transform matrix (i.e., the lifting scheme operations as described above). First, **d** is obtained by taking the forward wavelet transform of (4) of a given **c**, and then the inverse wavelet transform of (5) of **d** at each level *j* recursively results into a polynomial representation of , which is differentiated to find derivatives. The computational cost of this approach is approximately equal to that of calculating the wavelet coefficients, where neither nor is explicitly formed, thanks to the lifting scheme. It can also be shown that the maximum error of calculating *q*th-order derivative of *u*(**x**) is *O*(*ε*^{1−q/n}), where *n* is the order of the polynomial (e.g., Vasilyev and Kevlahan 2005; Alam 2006).

### d. Grid adaptation methodology

For the present work, the decomposition in (1) is associated with a multilevel grid, where there is a one-to-one correspondence between the grid points and wavelets. An example of a uniform three-level grid is shown in Fig. 1b, where a grid at level *j* is composed of a coarser grid and a complementary grid such that (or ). In other words, we have

where denotes an arbitrary coarsest level grid, containing *n _{x}* ×

*n*points such that

_{z}*n*= 2

_{x}

^{j}^{−1}

*m*and

_{x}*n*= 2

_{z}

^{j}^{−1}

*m*. For the two-dimensional example grid shown in Fig. 1b, we have

_{z}*j*=

*j*

_{min}+ 2 and

*m*=

_{x}*m*= 3. There are two advantages with the nonlinear approximation in (2). First, the numerical error associated with the approximation (2) is

_{z}*O*(

*ε*). Therefore, one has the flexibility and guarantee of prescribing a priori error tolerance. Second, a fraction of the grid points at level

*j*is used for calculating and reducing

*ε*by a factor of 2 does not imply doubling the number of grid points [e.g., (3)]. A more detailed theoretical discussion of this nonlinear approximation is given by DeVore (1998). Moreover, for a piecewise smooth function

*u*(x), the given tolerance

*ε*determines the maximum level

*j*or the minimum scale 2

^{−j}. Therefore, the wavelet methodology provides a natural framework for grid adaptation, where an arbitrary coarsest level grid is extended by adding fractions of complementary grids that are associated with the nonlinear approximation based on the given tolerance

*ε*. The grid adaption methodology is now outlined assuming that the solution u is exactly known. Let

*j*=

*j*

_{min}.

The significant solution and the associated adaptive grid are formed by discarding all nonsignificant points for which .

The refinement of follows by constructing an adjacent zone that contains only those points from the complementary grid , which belong to a set of nearest neighbors of only significant points from .

The adapted grid is formed by tuning the adjacent zone so that wavelet transforms and their inverse are reproducible. This is known as perfect reconstruction process.

If then steps 1 to 3 are repeated with

*j*←*j*+ 1. Otherwise, is taken to be the final adapted grid.

The above procedure results into an adaptive grid that contains only points. It has been verified that this is much smaller than the total number of grid points *N* = 2* ^{jd}* on the uniformly refined grid at scale 2

^{−j}for functions with highly localized features such as the velocity or the vorticity field (Vasilyev and Kevlahan 2005; Alam et al. 2006).

A principal objective of this study aims to develop a multiscale methodology for atmospheric modeling. The set of governing equations that form the dynamical core of the proposed model is now presented.

### e. Governing equations—General form

At the heart of an atmospheric model, a set of conservation principles forms a coupled set of PDEs that must be solved simultaneously. A system of prognostic equations describing an atmospheric motion can be compactly written as

where **u** is a two- or three-dimensional velocity vector, **Ψ** is a vector of *d* state variables including **u**, and is a vector that represents all forces, feedbacks, sinks, and sources (e.g., subgrid-scale or pressure gradient terms). Typically, the components of **Ψ** are density, velocity, potential temperature, and/or trace gas concentration (e.g., see chapter 2 in Pielke 2002).

The actual components of **Ψ** and depends on the specific application. The system in (6) forms the dynamical core of the model, which is currently being tested along with Boussinesq approximation for the purpose of comparison with similar models, where *ρ* denotes density and 0 < *η* < 1 is a small number (e.g., Defant 1951; Martin and Pielke 1983; Alam and Lin 2008). For example, Alam (2010) verified the present model for a 3D sea-breeze circulation system. The present paper examines the computational benefits of this development using 2D simulations.

### f. Boundary conditions

The model is designed for implementing three types of physical boundary conditions—a periodic condition [e.g., ], a Dirichlet condition [e.g., **Ψ**(**L**) = **Ψ**_{0}], or a Neumann condition [e.g., ], where **L** = (*L _{x}*,

*L*),

_{z}**x**= (

*x*,

*z*), and

**n**denotes the outward unit normal corresponding to the boundary. A nonreflecting boundary condition is used to eliminate nonphysical waves reflected from a solid boundary (e.g., Lane 2008).

### g. The numerical scheme

#### 1) Time integration

First, the system in (6) along with Boussinesq approximation is discretized using a Crank–Nicolson scheme:

which is second-order accurate in time, unconditionally stable for a linear equation, and there is no numerical dissipative error (e.g., chapter 4 of Tannehill et al. 1997). Note that the fully implicit treatment of advection is a distinct feature of the present development. For simplicity, we write (7) using a compact notation .

Second, pressure gradient forces **∇***P ^{n}*

^{+1}in (6) have been treated implicitly in (7) so that the time evolution of the velocity satisfies the divergence free condition

**∇**·

**u**

^{n}^{+1}= 0, such that

where **u**^{*} is the solution of (6) for which **∇** · **u**^{*} ≠ 0. The Crank–Nicholson time integration in (7) (e.g., for advection), requires solving a nonlinear system [(7)] at each time step *n* + 1. A multilevel solver is proposed in this paper, where denotes a discretization of (7) at level *j*. The pressure equation is solved with a multigrid solver.

#### 2) The FAS–JFNK methodology

For the dynamical core of the present model, benefits of the Jacobian-free Newton–Krylov (JFNK) method (e.g., Knoll and Keyes 2004) are combined with that of the full approximation scheme (FAS; chapter 8.3 in Wesseling 2004) so that fastest waves are either eliminated or damped with a fully implicit method (FAS–JFNK; e.g., Reisner et al. 2001; Bernsen et al. 2010). Theoretical details of both JFNK and FAS can be seen in cited references (e.g., Knoll and Keyes 2004; Wesseling 2004, and the references therein), and the computational power of these algorithms are well established such as in CFD. Neither of these two approaches got a combined usage in CFD, nor are their benefits well established in atmospheric modeling. The FAS–JFNK methodology is now presented briefly.

Let **u*** ^{k}* be an approximate solution of (7). The present JFNK method solves the linear model of (7) using a Krylov method with iterations [e.g., a generalized minimal residual algorithm (GMRES); e.g., Saad and Schultz (1986)], and completes a Newton step,

**u**

^{k}^{+1}=

**u**

*+*

^{k}**s**, until a stopping condition is satisfied. Instead of calculating the Jacobian matrix explicitly, the matrix-vector product is approximated via a Fréchet derivative [e.g., Eq. (2) in Reisner et al. 2001]. In the present context, JFNK is used as a relaxation method at intermediate levels of FAS iterations, where the iteration is stopped when either the condition is satisfied or steps of Newton iterations are completed, whichever occurs the first. In all presented experiments, Krylov steps to 5,

*γ*= 0.5, and Newton steps or 3 are used during a relaxation sweep. Note that not preconditioning matrix is constructed, which is a significant difference the way actual JFNK method would be used.

Denoting the JFNK method symbolically by , let us now present the FAS-JFNK solver below. In the following algorithm, the symbol means the truncation of the multiscale decomposition in (1) of *u* at level *j*. Following Wesseling (2004), the FAS equations are defined by at all levels <*j*, and is used at level *j*:

begin at level

*j*;if

*j*=*j*_{min}, then call ;otherwise call on the coarser level

*j*− 1, where ;end FAS-JFNK (…) if tolerance is satisfied.

Numerical experiments have verified the complexity—a linear increase of the CPU time with an increase of the number of grid points —of the proposed multiscale modeling approach that uses the AWCM and FAS–JFNK. Note that the use of a JFNK-type Krylov method in smoothing sweeps is a distinct feature of this development compared to the classic FAS method. In all numerical experiments, the present FAS solver have not used any anisotropic coarsening or refinement to accelerate convergence, which is one important advantage of using the Krylov method instead of the Jacobi iteration. Let us now present numerical results verifying the proposed model.

## 3. Model verification

### a. Analytical verification

For an analytical verification, the right-hand side of (6) is defined by . Then, the Crank–Nicolson method solves the nonlinear system:

at each time step *n* + 1, where **f** is the information from time step *n*. This is an ideal toy model for numerical verification because the dynamical core of the proposed model uses a FAS-JFNK algorithm to solve the above system at each time step. To examine the convergence, we define a Reynolds number , where *U* = 1 m s^{−1}, *L* = 1 km, and a CFL number , where Δ*x* is the finest resolution. The convergence was examined for 10 ≤ Re ≤ 10^{5} and 1 ≤ CFL ≤ 7, as well as for varying *ε*.

For this toy model, **f** is defined exactly by substituting into the above equation a given **Ψ**^{n}^{+1} = (*u ^{n}*

^{+1},

*υ*

^{n}^{+1})

^{T}that represents a Lamb–Oseen vortex (e.g., see chapter 13.1 in Saffman 1992) with a maximum velocity of 1.6 m s

^{−1}near the center of the domain that extends 10 km in each directions. Results in Fig. 2 correspond to Re ~10

^{5}and CFL ~ 6.6. The numerical solution

*u*

^{n}^{+1}, the exact solution , and the corresponding adapted grid is presented in Fig. 2, where Fig. 2c shows that the grid is refined locally, and is adapted to the solution (e.g., Fig. 2a) so that high-velocity gradients are resolved. The simulation was started with a uniformly distributed 5 × 5 grid with Δ

*x*= Δ

*z*= 2.5 km at level

*j*= 1. The grid was refined locally until the level

*j*= 9 is reached, and all wavelet coefficients with a measure less than

*ε*= 10

^{−4}were discarded. We have found that the number of points in the final adaptive grid is , where the finest local resolution has Δ

*x*≈ Δ

*z*≈ 9.7 m near the center of the domain, where the local CFL number, based on the maximum velocity 1.6 m s

^{−1}, is about 6.6. Such a high resolution without adaptivity requires a uniformly refined 1025 × 1025 grid or a total of

*N*= 1 050 625 grid points. Clearly, the adaptivity reduces this huge number of grid points by a factor of about 625 (or ). However, the quality of the numerical solution on the adaptive grid is comparable with that of the exact solution on the full grid as depicted in Figs. 2a,b, where the relative rms error is 1.13 × 10

^{−4}. Furthermore, running the model for 10

^{−6}≤

*ε*≤ 10

^{−1}confirms that the relative error remains proportional to the tolerance

*ε*, and the CPU time remains proportional to the number of grid points (e.g., Fig. 2d).

The present test confirms that (i) the Crank–Nicolson module calculates the solution accurately at each time step and (ii) a high resolution can be achieved by reducing the number of grid points drastically. However, this analytical test does not solve an equation that governs the dynamics of an atmospheric motion and, therefore, we need to test our model with a more acceptable case that is often used as a benchmark problem.

### b. Comparison with the Smolarkiewicz’s model

Smolarkiewicz’s deformational flow (Smolarkiewicz 1982) is one of the standard tests for the verification of advection schemes (Sykes and Henn 1995), where a set of counter-rotating cells distorts an initial scalar distribution. In this flow, the time evolution of the scalar field is a filamentary spiral, where the number of spiral turns increases and the width of the filamentary arm decreases (e.g., Staniforth et al. 1987). Therefore, as time proceeds, the grid must be refined locally so that the filamentary spiral of the scalar field is resolved.

The present simulation used a domain that extends 2 km in both horizontal and vertical directions. The deformational velocity field **u** = (*u*, *w*) is obtained by solving the nonlinear advection problem with an initial condition that is defined by the streamfunction **Ψ**(*x*, *z*) = *A* sin(*kx*)cos(*kz*) (e.g., Staniforth et al. 1987), using *A* = 10^{−3}*L*^{2}/2*π*, *k* = 2*π*/*L*, and *L* = 1 km, where and . These parameters result in a velocity field with a maximum of 1 m s^{−1} in both directions. The time evolution of the potential temperature *θ*(*x*, *z*, *t*) is obtained by solving with *κ* = 10^{−1} m^{2} s^{−1} that is adjusted from numerical experiments so that the solution remains stable for all time. This value of the thermal diffusion coefficient *κ* is, however, 100 times smaller than that was used by Martin and Pielke (1983) in a sea-breeze model. Note that the numerical solution of the Smolarkiewicz’s model with *κ* = 0 on a uniformly refined Eulerian grid remains stable only for a short time (Smolarkiewicz 1982).

The potential temperature at *t* = 0 min and its time evolution at *t* = 33 min are presented in Figs. 3a,b, where the minimum resolution is 3.9 m in both directions. As depicted in Fig. 3c, the grid is dynamically adapted to resolve the filamentary structure of the potential temperature. The number of points in the grid at *t* = 0 min is 7808 and that at *t* = 33 min is 25 344. Clearly, the number of grid points increased in 33 min by about a factor of 3 so that the small filamentary structures of the spiral scalar field are resolved. Figure 3d shows that increases to control the error dynamically [e.g., (3)] by adding–deleting new grid points. This adaptive simulation result agrees with those results that are presented by other researchers (e.g., Staniforth et al. 1987; Bacon et al. 2003), thereby providing a qualitative cross check on the model’s accuracy.

## 4. Dry thermal simulation

### a. Governing equations

To verify the ability of the proposed model for meteorological simulations, the time evolution of warm thermals in a dry atmosphere is considered, which offers much insight into more complicated atmospheric dynamics such as mixing and redistribution of heat and constituent species (Manton 1978). The dynamics of dry thermals with no condensation or evaporation can be modeled with the following equations (e.g., Lane 2008):

where all symbols are defined in Table 1. Similar to Lane (2008), diffusion terms are used for numerical modeling of subgrid-scale effects. These equations are solved using periodic conditions in the horizontal direction for all variables, and homogeneous Dirichlet conditions in the vertical direction for all velocity components. For the Buoyancy field, conditions in the vertical direction are *b* = 0 at the bottom boundary and at the top boundary.

### b. Comparison with known thermal dynamics

The present two-dimensional dry simulation in a 20 km × 10 km domain is similar to that of Bryan and Fritsch (2002). Assuming a stationary flow initially (i.e., *u* = 0 and *w* = 0), a warm thermal of radius 2 km is placed along the center line, *x* = 0, of the domain at 3 km above the bottom boundary (i.e., *z* = 3 km). The maximum buoyancy is 9.8 × 10^{−2} m s^{−2} at the center of the thermal, which decreases exponentially to 0 m s^{−2} outside the edge of the thermal. This buoyancy anomaly is equivalent to a potential temperature perturbation 3 K, which is similar to thermals in a boundary layer (Lane 2008). The evolved state of the thermal and that of the vertical velocity at *t* = 12 min is presented in Figs. 4a,b. The spanwise vorticity field at *t* = 12 min in Fig. 4c shows that the region of steep buoyancy gradient is associated with that of large velocity gradient (e.g., Carpenter et al. 1990; Lane 2008). Figure 4d shows that the grid has been adapted to local features of the flow, which is the achievement of this development, and its potential advantage for more realistic systems. These results have a good qualitative agreement with those that were predicted by various other numerical models (e.g., see, Carpenter et al. 1990; Wicker and Skamarock 1998; Bryan and Fritsch 2002; Lane 2008, and the references therein).

As discussed in Carpenter et al. (1990), the buoyancy gradient steepens as a result of nonlinear advection, thereby forming small-scale coherent structures that require a fine grid. Vertical distributions of the buoyancy field at *t* = 0, 4, 8, 12 min along the center of the domain are presented in Fig. 5a. Clearly, a steep vertical gradient of the buoyancy field occurs in different vertical locations at various times. It is thus necessary that the refinement or the coarsening must be dictated by the time evolution of local features. Clearly, statically refined or nested grid models are not suitable for this problem (e.g., Jablonowski 2004). The AMR model described in Skamarock et al. (1989) would compare solutions at two different resolutions before refining or coarsening the grid (e.g., Jablonowski 2004; Blayo and Debreu 1999). In contrast, the present model uses a coarse-grid approximation of the coherent thermal dynamics that is sufficient to identify the portion of the grid that must be coarsened or refined. In other words, grid points are rearranged, deleted, and/or added so that the buoyancy field is resolved at each time step using a near-optimal CPU time.

In Fig. 5b, the elapsed CPU time (s) at each time step and the number of adaptive grid points that are used at each time step are plotted against the simulation time *t* (min), where and CPU time are normalized by the number of grid points and elapsed CPU time corresponding to the last time step at *t* = 12 min respectively. We see that both the CPU time and the number of grid points increase approximately at the same rate with respect to the simulation time *t* (min). This indicates that, if the buoyancy field steepens as depicted in Fig. 5a, the model adds new grid points to resolve the buoyancy gradient, thereby increasing both and CPU time. To explain this increase in and CPU time (s), we have also compared the elapsed CPU time (s) at each time step with , showing a linear relationship, as depicted in Fig. 5c. These results verify that the performance of the model is asymptotically optimal for a dry thermal simulation.

The accuracy of this simulation can be quantified by calculating kinetic energy, , and potential energy, , where *dA* is an area element, and the buoyancy is related to the potential temperature by *θ* = *θ*_{0}(1 + *b*/*g*). The total energy is defined by *E*(*t*) := *E _{k}*(

*t*) +

*E*(

_{p}*t*). Relative deviations of these energies with respect to, and normalized by their initial values

*E*(0) and

_{k}*E*(0) are presented in Fig. 6, showing that the maximum total energy error is

_{p}*O*(10

^{−3}). The similar accuracy in energy errors between the present model and the piecewise parabolic model (PPM) of Carpenter et al. (1990) suggests the promise of multiscale approximation used here.

### c. Comparison with a reference finite-difference model

For a quantitative analysis of the proposed wavelet model’s performance, a finite-difference method is considered as a reference model, which employs a finite-difference scheme using a uniformly refined grid without any adaptivity. We have used the same time integration scheme for the reference and the wavelet model. The comparison between the result of the wavelet model and that of such a reference model would assess the improvement in accuracy and CPU time that are achieved by adaptivity. Let us now compare the performance between the wavelet and the reference model with respect to the simulation presented in section 4b.

First, we compare the number of grid points . The wavelet model uses at *t* = 0 min, which increases to at *t* = 12 min, so that steep buoyancy gradients are resolved using a horizontal resolution Δ*x* = 19.5 m and a vertical resolution Δ*z* = 9.8 m. The reference model requires at *t* = 12 min for resolving the steep buoyancy gradient, and the same grid was required at every time step. Clearly, a drastic reduction (e.g., at least about 93%) of grid points is possible with the wavelet model due to adaptivity. Second, we compare the CPU time. We have found that the wavelet model takes about 96% less CPU time compared to the reference model. Both models adopted a similar unconditionally stable time integration scheme for this simulation with a CFL number 10. This large CFL number confirms that the high CPU time required by the reference model was due to the fine spatial resolution.

## 5. Application to vortex pair dynamics

In a stably stratified atmosphere, interactions between turbulence and buoyancy affects the merging of a like-sign vortex pair—a primary vortex interaction associated with the inverse energy cascade from small to large scales (e.g., Hooke and Jones 1986; Iida et al. 2002; Brandt and Nomura 2007). The effect of stable stratification on the dynamics of a vertically oriented vortex pair was considered in most studies (e.g., Hardenberg et al. 2000), whereas that of a horizontally oriented vortex pair was, to the best of knowledge, not well investigated (e.g., Brandt and Nomura 2007). High-resolution numerical simulations of Brandt and Nomura (2007) using a 2048 × 2048 grid indicated that the generation of intermittent filamentary vortex structures is a direct result of stable stratification via baroclinically generated torque.

In this section, we study the performance of the proposed multiscale model for simulating the vortex interaction phenomena presented in Brandt and Nomura (2007). The numerical simulation was done using the set of equations in (8)–(11), where the initial velocity field consists of a corotating Lamb–Oseen vortex pair. The boundary conditions are same as those used for thermal simulation except periodic conditions are used in both directions for velocity. Each initial vortex has a circulation Γ = 1 km^{2} s^{−1}, and the initial distance between two vortexes is *b*_{0} km with *a*/*b*_{0} ~ 0.157, where the radius of each vortex is *a* km. With an initial separation distance *b*_{0} = 1 km, the initial rotational velocity is about 159 m s^{−1}. The computational domain extends 24 km in both horizontal and vertical directions.

First, the evolution of the vortex pair in a neutral environment, Fr = ∞; that is, the buoyancy frequency *N _{b}* is 0 has been investigated, where the Froude number is As depicted in Fig. 7 (top row), the dynamics of the vortex pair in a neutral flow is in good qualitative agreement with previously known experimental and numerical results, thereby suggesting the convergence of the method. Second, an ambient stable stratification is imposed to the model that is characterized by Fr < ∞, where the time evolution of the vortex pair—presented in Fig. 7 (bottom row)—has a good agreement with that is presented in Brandt and Nomura (2007).

The effect of stable stratification is studied by simulating the vortex merger at various Fr = 0.7, 1.5, 3.0, and ∞. In Fig. 8 the vorticity field and the corresponding adapted grids are presented. Clearly, all grids are sparsely populated, in other words, grid points are dynamically concentrated in the region of coherent motion, where strong velocity gradient occurs. Using a reference model as described in section 4c, we have found that the wavelet model saves about 98% computational work for those results that are presented in Fig. 8. This extreme reduction of computational work is a result of the high spatial intermittency of the solution. This example represents two important phenomena simultaneously: the transfer of energy from small to large scales via vortex merger and that in the opposite direction via buoyant interaction of vortexes. The result indicates this model’s potential for more realistic atmospheric simulations.

## 6. Simulation of flow over isolated topography

The numerical simulation of stratified flows over topographic obstacles is a significant challenge, where the terrain-following coordinates are commonly used in most atmospheric models (e.g., see Schär et al. 2002). In this section, we study how the proposed model can be extended to simulating stratified flows past an isolated topography.

The boundary conditions on the topography are imposed by parameterizing the body force exerted by the topography using the volume-averaging technique proposed in Whitaker (1996). The approach assumes that the topographic obstacle is a porous media having a permeability **K**. A brief derivation of the body force term is presented in the appendix. The body force term implements the homogeneous Dirichlet boundary condition approximately on the surface of the topography.

An isolated topography is defined by the Witch of Agnesi curve:

where the mountain height is *H* = 400 m and the width is *a* = 1000 m in the model domain, 40 km × 10 km. The flow is initially isothermal with a temperature of *θ*_{0} = 300 K and is driven by the horizontal wind speed *U* = 9 m s^{−1}. Horizontal boundary conditions are periodic and a vanishing vertical gradient is imposed on the model top for all variables. The inverse Froude number, based on the mountain half-width is , which means that nonhydrostatic waves are expected to propagate.

Numerical verification is also done by comparing results with the advection test of Schär et al. (2002) and with the immersed boundary approach of Lundquist et al. (2010) with an equivalent choice of parameters, but the results are not presented. Note, however that the topography of the present test differs from that of the advection test of Schär et al. (2002) and the governing set of equations also differ from both the experiments of Schär et al. (2002) and that of Lundquist et al. (2010). This comparison aims to understand an approximate range of values for **K**, and has found that the error for imposing boundary conditions on the velocity field remains within the wavelet threshold *ε*—the error tolerance for adaptivity—for 10^{−5} ≤ **K** ≤ 10^{−7}. This body force term has also been treated by the Crank–Nicolson method.

Nonhydrostatic wave propagation past an isolated topography is simulated, and the vertical velocity is presented in Fig. 9, showing that the present model tracks the topography and wave propagation by placing grid points adaptively, where it is necessary. This numerical study presents how a topographically driven flow can be simulated using the proposed multiscale approach. Although the results presented in Fig. 9 exhibits the expected nonhydrostatic wave propagation, more verification is necessary to quantify the numerical error. However, this result explains how the present development can be extended to simulate orographically driven flows.

## 7. Summary and future developments

### a. Summary

This paper has investigated the possibility and advantages of a novel multiscale methodology for atmospheric modeling such that the solution is represented by only a fraction of wavelet modes. These wavelet modes are associated with the significant dynamics that must be resolved in a numerical model. All wavelet modes, having coefficients below a threshold, are associated with the dynamics that must be parameterized, or can be discarded, depending on the specific application. The resolved significant part represents localized small-scale features on a large-scale background, and results into an adaptive mesh. To the best of knowledge, this is the first attempt proposing such a multiscale modeling approach in the field of atmospheric modeling, where wavelets are used to model multiscale features. The paper has also investigated the numerical simulation of coupled meteorological equations, where the time evolution of nonlinear dynamics has been treated implicitly with the aid of two advanced algorithms: the Krylov method and the full approximation scheme. The implicit algorithm using an optimal CPU cost is a clear step forward in the field of atmospheric modeling.

Presented numerical experiments were used previously by other researchers, and these examples (i) exhibit isolated small-scale phenomena, (ii) demand for a large amount of computing power, and (iii) represent multiscale processes that are also challenging for numerical atmospheric models (e.g., Carpenter et al. 1990; Brandt and Nomura 2007). We have verified numerically that the computational cost of the proposed model is much less than that of a classic finite-difference model for simulations presented in this paper. This high reduction of CPU time is a result of discarding a large number of wavelet modes associated with the nonsignificant motion. We have also verified numerically that the number of grid points in the adaptive grid is proportional to the CPU time that is required for calculating the solution. This means that local features are computed with a locally refined mesh without overburdening the global computational cost. The extreme computational performance of this model with chosen examples encourages us that advanced computational algorithms may add further benefits to the field of atmospheric modeling.

### b. Future developments

The results of this work suggests some useful future research directions. Since the present two-dimensional results show that the wavelet methodology has a smart ability of resolving localized dynamics, an extension of the present work for 3D simulations is a clear next step. The wavelet methodology can also be implemented on massively parallel architectures to speed up 3D simulations. Since atmospheric phenomena spans over a wide range of length scales, it may not be possible to resolve the solution at all scales, particularly in 3D. One must still parameterize certain small-scale features. An understanding of a suitable parameterization technique adapted to the present methodology is an important future development. An appropriate parameterization scheme combined with the proposed adaptive computational model may provide acsolution to some of the challenges facing state-of-the-art atmospheric models.

The set of governing equations for present experiments is an approximation to the original anelastic set of Ogura and Phillips (1962), neglecting vertical variation of potential temperature except in the leading-order contribution to the buoyancy term (Durran and Arakawa 2007). This approximate form is commonly used for modeling mesoscale phenomena such as sea-breeze circulation, dry convection, etc. (e.g., Defant 1951; Martin and Pielke 1983; Alam and Lin 2008; Lane 2008), or modeling the atmospheric boundary layer (Moeng 1984). The choice of compressible, anelastic, or approximate anelastic form depends on specific applications or target atmospheric conditions. The dynamical core of this model consists of the system in (6) that is written in the compressible form. The proposed implicit time integration scheme aims at solving the system (6) in its compressible form. However, a concise understanding of the computational advantages, using the Boussinesq approximation also allows comparisons with many other models. The verification of this new development using its compressible form in (6) is a potential future development of this study.

## Acknowledgments

This work was partly supported by a grant from the Industrial Research and Innovation Fund (IRIF) of the Government of Newfoundland and Labrador and a grant from the Natural Sciences and Engineering Research Council of Canada. This work was made possible by the facilities of the Atlantic Computational Excellence Network (ACEnet; available online at http://www.ace-net.ca/wiki/ACEnet). The author thanks the three anonymous reviewers for their useful suggestions.

### APPENDIX

#### A Technique for Modeling the Topographic Effect

Let us assume that the topography is a submerged porous media with a permeability **K**, where only a fraction of an elementary volume is occupied with fluid. The procedure starts by defining two averages of a fluid property *ψ*:

such that

Hence, the effect of the submerged topography can be modeled—by taking the average 〈·〉 of the momentum equation over an elementary volume containing fluid and porous media—with the addition of a body force term:

where and denote fluctuations of pressure and velocity within the averaging volume , and *A* is the interface between solid and fluid. Clearly, *A* = 0 in the fluid region, where topography is absent, and hence **f** = 0. Thus, the body force works only on the topography.

The volume-averaged (〈·〉) Euler equation is obtained by dropping the viscous terms from the volume-averaged Navier–Stokes equation (VANS), which is derived in (Whitaker 1996), assuming that the variation of *λ* is negligible. Note that the viscous term is also known as the Brinkman correction to porous media flow. Without reproducing the detailed derivation of VANS, we can write the volume-averaged Euler equation:

Clearly, the body force, **f**, has been included to model the effect of the solid–fluid interaction. This formulation is extended to the present development by assuming that the bottom topography is a porous material.

In the present development, a constant eddy viscosity is used to parameterize the subgrid-scale stress term. Following (Whitaker 1996), the body force term is parameterized by