Abstract

An atmospheric convective system may be modeled as an ensemble of discrete plume elements. A representation of decomposited plumes provides the basis for mass-flux convective parameterization. A dry version of such a prototype model is constructed in a two-dimensional horizontally periodic domain. Each discrete plume element is approximated by a horizontally homogeneous segment such that the whole system is given by segmentally constant approximations (SCA) in the horizontal direction for each vertical level in a nonhydrostatic anelastic model (NAM). The distribution of constant segments is highly inhomogeneous in space and evolves with time in a highly adaptive manner.

The basic modeling strategy from a physical point of view is to activate new segments vertically upward with time when a convective plume is growing and to deactivate segments when a plume event is over. The difference in physical values crossing segment interfaces is used as a criterion for numerically implementing this strategy. Whenever a large difference is found, the given interface is stretched vertically by subdividing an existing segment into two. In turn, when a segment interface difference is found below the threshold, the given interface is removed, thereby merging the two segments into one.

This nonhydrostatic anelastic model with segmentally constant approximations (NAM–SCA) is tested on an idealized atmospheric convective boundary layer. It successfully simulates the evolution of convective plumes with a relatively limited number of segments (i.e., high compression) and with a much scarcer distribution of segments over nonplume environments (i.e., extremely inhomogeneous distribution of segments). Overall, this method compresses the size of the model up to 5 times compared to a standard NAM with homogeneous grid distribution without substantially sacrificing numerical accuracy.

1. Introduction

The mass-flux convective parameterization originally proposed by Ooyama (1971) and Arakawa and Schubert (1974) is currently employed by a majority of global models, both those for operational forecasts and for climate research. The original conceptualization for this mass-flux convective parameterization may be traced back to Riehl and Malkus (1958). They proposed that an ascending branch of the Hadley–Walker circulation consists of two distinctive subcomponents: an ensemble of well-localized fast-ascending elements called “hot towers” or deep moist convective towers and a slowly descending background environment. They argued that the existence of localized hot towers is crucial for qualitatively understanding the tropical heat budget. They furthermore quantitatively estimated the contribution of hot towers in the mean heat budget.

Riehl and Malkus’ (1958) hot towers may be considered a special case of plume models for atmospheric convection. Since publication of the seminal works by Stommel (1947, 1951), thermal “plumes” have been considered elementary building blocks both for shallow and deep atmospheric convection (cf. Simpson 1983a,b; Lilly 1983; Morton 1997). Subsequent extensive analog laboratory experiments for atmospheric convection (e.g., Morton et al. 1956; Scorer 1957; Turner 1962; see also Turner 1986 as a review) greatly promoted this idea.

The basic idea of the mass-flux convective parameterization is to represent the subgrid-scale convective processes within a grid box by an ensemble of those hot towers or more generally, plumes embedded in a homogeneous environment (cf. Fig. 1 of Arakawa and Schubert 1974). This idea is realized by introducing two equations for each physical variable within a grid box separately describing the two subcomponents: hot towers and environment. A prototype for such a closed formulation was provided by Riehl and Malkus (1958) for diagnostic purposes (cf. Yano 2009).

The most formal starting point for developing such a formulation is to suppose that a grid box itself is described by a full physical system. This is the basic idea of superparameterization, which adopts a cloud-resolving model (CRM) as a full system (Grabowski and Smolarkiewicz 1999; Randall et al. 2003). A description of the ensemble of hot towers embedded in an environment is essentially obtained by discretizing a full system in terms of these discrete subcomponents.

The mass-flux parameterization achieves this goal by subdividing the gridbox domain into these subcomponents. Each subcomponent is approximated by a horizontally constant element surrounded by a horizontally homogeneous environment. We call this approximation the segmentally constant approximation (SCA). The idea of SCA may be traced back to Stommel (1947), who introduced a hypothesis of horizontal homogeneity in constructing his entraining plume model. The notion of the “top-hat” profile is explicitly adopted by Morton et al. (1956).

Arakawa and Schubert (1974) present the complete formulation for their parameterization by introducing further approximations and hypotheses (cf. Yano et al. 2005). However, it is hard to overemphasize that the very basic idea of this parameterization resides in subdividing the gridbox domain into various subcomponents under SCA. It is also important to emphasize that the idea of this subdivision under SCA is not necessarily limited to the one between hot towers (or plumes) and the environment; but can be generalized for various different processes such as convective downdrafts, or even various mesoscale subcomponents such as stratiform clouds.

The goal of the present paper is to present a model that is simply constructed by subdividing the gridbox domain into subcomponents. Unlike the standard mass-flux parameterization, this model however introduces no further approximations nor hypotheses. For this goal, the present study adopts a nonhydrostatic anelastic model (NAM) as a full system, because NAM is a formulation commonly adopted both for the large-eddy simulation (LES) and cloud-resolving modeling. The constructed model can be considered as a prototype for mass-flux convective parameterization.

The present work constitutes a natural extension of previous work by Yano et al. (2004a, 2005). Yano et al. (2004a) propose wavelet “compression” as a procedure for reducing the size of a full model such as NAM to a degree comparable to that of a standard parameterization. Yano et al. (2005) further generalize this idea by adopting mode decomposition as a procedure for developing a representation for subgrid-scale physical processes. Under this perspective, the mass-flux parameterization is based on segmentally constant mode decomposition. However, there is major difficulty in constructing a subgrid-scale representation on this basis.

Yano et al. (2005) emphasize the importance of completeness and orthogonality of decomposition modes to facilitate construction of a subgrid-scale representation from a full model. Once these conditions are satisfied, a given system can be truncated (compressed) extensively by following the principle of the multiresolution approach as in the case for wavelets (cf. Mallat 1998). This framework permits activating and deactivating modes in phase space along the evolution of the system in a flexible manner as extensively discussed in Yano et al. (2005). Unfortunately, segmentally constant decomposition does not satisfy completeness and orthogonality.

To circumvent this difficulty in the present study, activation and deactivation of constant segments are performed in the physical space in such a way that the total number of segments is always kept at a minimum. Although the developed methodology is still under the general framework of the multiresolution approach, the actual implementation turns out to be something more akin to adaptive mesh refinement (AMR; e.g., Berger and Colella 1989; Bell et al. 1994; see Baker 1997 as a review; see also, e.g., Dietachmayer and Droegemeier 1992; Skamarock and Klemp 1993; Hubbard and Nikiforakis 2003; Jablonowski et al. 2006; and St-Cyr et al. 2008 for atmospheric applications).

Thus, the present paper reports an initial development of a NAM under SCA (NAM–SCA), with adaptive activation and deactivation. A two-dimensional configuration for a dry version over a periodic domain is adopted. The vertical resolution is fixed in time, and refinement is performed only in the horizontal direction.

After introducing the basic model configuration, the formulation for NAM is reviewed in the next section. The formulation of the SCA system is presented in section 3. An algorithm for adaptive activation and deactivation of segments is introduced in section 4. As a preliminary test of the model, section 5 reports an application to a free-convection case proposed by Ayotte et al. (1996) for testing boundary layer parameterizations. Section 6 presents conclusions and details of numerical implementation are presented in the appendix.

2. A basic set of equations

a. Model configuration

The two-dimensional system is adopted. It is periodic in the horizontal direction x with a domain size (periodicity) L, and it extends vertically for a height H.

By following the basic principle of the mass-flux convective parameterization of Arakawa and Schubert (1974), we consider that the system consists of an ensemble of plumes,1 each of them with a constant value (i.e., SCA) at each vertical level. Such a geometrical configuration is schematically shown in Fig. 1. Here, six constant segment plumes are placed in the domain. The first two plumes (from left) horizontally extend for (ξ1, ξ2) and (ξ2, ξ3), respectively; and vertically extend from the surface to the height ζ6 and ζ4. However, the plumes are not necessarily originating from the surface. The last three plumes (from left) originate from the heights ζ3, ζ1, and ζ2, respectively, spanning for (ξ4, ξ6), (ξ6, ξ7), and (ξ7, ξ8). A plume may not have the same horizontal extent all the way to the top; the third plume from left spans for (ξ3, ξ5) from the surface to the height ζ3, then reduces its horizontal extent to (ξ3, ξ4) at z = ζ3, and then it continues upward to the top at z = ζ4.

Fig. 1.

Schematics for a geometrical configuration of the model. Physically, a system with an ensemble of plumes (shown by varying gray tones) is considered, which, in turn, is mathematically reinterpreted as a representation of the system under SCA.

Fig. 1.

Schematics for a geometrical configuration of the model. Physically, a system with an ensemble of plumes (shown by varying gray tones) is considered, which, in turn, is mathematically reinterpreted as a representation of the system under SCA.

In this manner, based on a plume picture, the whole system is decomposed into a set of horizontally homogeneous segments. Once such a geometrical decomposition is achieved, however, we no longer have to stay with an original plume picture any more. For example, the fourth segment from the left spanning ξ4xξ6 and ζ3zζ4 may be better considered (in a moist version) as a stratiform cloud. Also, keep in mind that each plume has a height dependence, although the schematics may indicate otherwise, as well as the whole environment indicated by white in Fig. 1.

The basic strategy of the present modeling is to perform the time integration of a full system under decomposition into horizontally constant segments under SCA. The system under SCA is derived in the next section. These segments are activated and deactivated with time by following the evolution of the ensemble of plumes, or the whole system, so that the geometrical configuration of the model is adjusted accordingly, as detailed in section 4. Before proceeding further, the basic formulation of the full system is first reviewed in the remainder of this section.

b. Basic formulation

The dry two-dimensional NAM is considered. All symbols with the subscript r refer to the reference background state depending on height z only. The pressure p and the potential temperature θ without subscript r refer to the deviations from these reference profiles, hence the total pressure is given by pr + p, for example. The reference density ρr is simply assumed a constant in the present study, although it will be included in the following formulations for a generalization in future studies.

The basic set of equations consists of the prognostic equations for the horizontal u and vertical w winds, and the potential temperature θ,

 
formula
 
formula
 
formula

and the mass continuity,

 
formula

with the acceleration of the gravity g. Here, (∂θ/∂t)flux represents surface flux specified later by Eq. (5.1). No diabatic heating is considered in the present study.

A thermodynamic boundary condition ∂θ/∂z = 0 is posed both at the surface (z = 0) and the top of the system (z = H). The rigid-lid boundary condition (w = 0) is assumed at the surface, whereas a free surface condition (p = 0) is assumed at the top. A radiative boundary condition proposed by Klemp and Durran (1983) reduces to the free surface condition in the limit of low wavenumbers. For this reason, we expect that the free surface condition loosely behaves like a top radiative boundary condition, but without a need to move to a wavenumber space as the original Klemp–Durran condition would require. This mitigates the necessity of introducing a sponge layer in combination with a limited segment number at the uppermost part of the model. No subgrid-scale parameterization, such as eddy viscosity, is included in the present version for the sake of presenting the SCA model in its simplest form (cf. Margolin et al. 1999).

c. Poisson problem

The evolution of the nonhydrostatic anelastic system may be evaluated by time integrating Eqs. (2.1a), (2.1b), and (2.1c) for the three variables u, w, and θ with a given pressure p. However, the main problem in performing this time integration is that there is no equation immediately available for evaluating p, thus an auxiliary equation must be introduced (see, e.g., Bernardet 1995). The desired equation is obtained by multiplying the linear operators (∂/∂x)ρr and (∂/∂z)ρr in Eqs. (2.1a) and (2.1b), respectively, and adding them together. The obtained divergence tendency must vanish because of the mass continuity (2.2), thus

 
formula

where v = (u, w).

In summary, the NAM system consists of a redundant set of Eqs. [(2.1a), (2.1b), (2.1c), (2.2), and (2.3)]. The standard procedure is to time integrate Eqs. (2.1a), (2.1b), and (2.1c) with the pressure p diagnosed by Eq. (2.3). On the other hand, NAM–SCA takes a rather unusual set consisting of Eqs. (2.1b), (2.1c), (2.2), and (2.3) for the reason explained in the next section.

3. SCA: Segmentally constant approximation

a. Model geometry2

We represent a physical system under an SCA (cf. Fig. 1), which consists of approximating a full system by, say, a n(z, t) set of horizontally constant segments defined over the intervals [xj−1/2,b(z, t), xj+1/2,b(z, t)] with an index j for the segments spanning for j = 1, 2, … n(z, t) (Fig. 2). Note the number n(z, t) of segments depend both on height and time. Being consistent with the periodic boundary condition, a cyclic condition is introduced on the interface positions, that is, xn+j−1/2,b(z, t) = xj−1/2,b(z, t) + L. As a result, any physical variable ϕ(x, z, t) is approximated as

 
formula

with a segment-mean value ϕj(z, t) defined by

 
formula

(cf. Fig. 2) and an indicator Ij(x, z, t) for the jth segment defined by

 
formula
Fig. 2.

Application of SCA to an arbitrary prognostic variable ϕ given as a function of horizontal coordinate x. Under SCA, the variable ϕ is defined by constant values ϕj over specified intervals (xj−1/2,b, xj+1/2,b) with j = 1, … , n.

Fig. 2.

Application of SCA to an arbitrary prognostic variable ϕ given as a function of horizontal coordinate x. Under SCA, the variable ϕ is defined by constant values ϕj over specified intervals (xj−1/2,b, xj+1/2,b) with j = 1, … , n.

b. Prognostic equations

The equation for each segment is obtained by averaging the original full system over the given segment. For this purpose, we first rewrite the prognostic Eqs. (2.1b) and (2.1c) in a general flux form,

 
formula

in terms of an unspecified prognostic variable ϕ (either w or θ) and forcing (source) term F. Equations (2.1b) and (2.1c) are recovered by appropriately specifying the corresponding term F. We integrate the system given by Eq. (3.3) over the jth segment that spans the interval (xj−1/2,b, xj+1/2,b). The derivation is essentially the same as when the zeroth-order finite volume method is applied (cf. Godunov 1959; section 5.6.1 of Durran 1999; LeVeque 2002).

The final expression obtained for the prognostic equation under SCA is

 
formula

Here,

 
formula

is a fractional length occupied by a jth segment, and variables (e.g., ϕj−1/2,b) with subscript b refer to those at the interface x = xj−1/2,b between the two segments j − 1 and j. The segmentally averaged forcing Fj is defined in an analogous manner to Eq. (3.2a). The vertical flux ()j has been approximated by wjϕj, following a standard mass-flux approximation [cf. Eq. (3.3) of Yano et al. 2004c].

The choice of the physical values at the segment interfaces remains arbitrary under SCA. By following a standard approximation introduced in plume-dynamics studies as well as in the mass-flux formulations (cf. Asai and Kasahara 1967; Arakawa and Schubert 1974), the upstream (upwind) approximation may be taken, hence

 
formula

This also corresponds to the simplest stable advection scheme available under the finite-volume approaches.

c. Horizontal velocity

Both the vertical velocity w and the potential temperature θ are evaluated prognostically by casting Eqs. (2.1b) and (2.1c) into the form of Eq. (3.4). On the other hand, the horizontal velocity u is treated differently as detailed in the present subsection.

By directly integrating the mass continuity (2.2) horizontally with a given segmented constant w, it is found that u is linearly segmented, for example,

 
formula

over a jth segment. By setting x = xj+1/2,b in this formula, we obtain a formula for evaluating the horizontal velocity at a next segment interface, j + ½, from a previous one, j − ½,

 
formula

A chain application of Eq. (3.6b) defines the horizontal velocity at every segment interface except for the domain mean u left undetermined.

As in any cloud-resolving modeling under periodic boundary conditions, there is no obvious way to determine u; thus it must be posed externally [cf. Eq. (8) of Grabowski et al. 1996; see also Eq. (7) of Grabowski 2004]. In the present study, we simply set u = 0.

d. Poisson problem

In the present work, the Poisson problem (2.3) is solved on a regular grid without applying SCA. This choice is consistent with the main goal of the present paper for presenting the new model, NAM–SCA, in its simplest setup. For this purpose, all the physical variables represented under SCA are transformed back to a “continuous” representation with the help of Eq. (3.1). The horizontal velocity is defined continuously by Eq. (3.6a).

Handling the discontinuities arising from SCA in evaluating the source term in the Poisson Eq. (2.3) is not as straightforward. We avoid this difficulty by assuming that all the physical variables are numerically defined at every grid point. A horizontally homogeneous grid with a resolution Δx is assumed for this purpose. Here, Δx is the minimum segment size (cf. section 4a). A standard finite-difference approach is then applied for evaluating the source term as detailed in the appendix. Although the procedure may not be justified in a mathematically rigorous sense, our tests repeating runs with the resolution Δx/3 in the Poisson problem show that the pressure calculations are not sensitive to the numerical resolution adopted.

e. Numerical implementation

A closed set of the system consisting of Eqs. (2.1b) and (2.1c) cast into a form of Eq. (3.4) is thus numerically integrated in time along with the diagnostic Eqs. (3.6) and (2.3).

The system is discretized in the vertical direction by Nz layers with a depth Δz = H/Nz. The bottom interface (full level) of the first layer is designated by a vertical level k = 0, and the top interface of the uppermost layer by k = Nz, corresponding to the model surface (z = 0) and the model top (z = H), respectively. In solving the Poisson problem (2.3), we assume that the physical variables are defined at a middle level of each vertical layer, which we designate as a half level with an index of a half integer, ½, 1 + ½, … , Nz − ½ (cf. Fig. 3). The half-level index is also used for identifying a model vertical layer in the following section. Additional numerical details are in the appendix.

Fig. 3.

Vertical configuration of the model. Full and half levels are indicated by full and dashed horizontal lines, respectively.

Fig. 3.

Vertical configuration of the model. Full and half levels are indicated by full and dashed horizontal lines, respectively.

4. Adaptive activation and deactivation of segments

As emphasized in the introduction, the present work has been inspired from the multiresolution approach (cf. Mallat 1998), which enables us to simulate the time evolution of a physical system under a strong truncation or compression by deactivating unnecessary modes and activating newly generating modes with time. The present work intends to perform such activation and deactivation with time by adopting a similar strategy in the real space.

As in the case for the wavelet-based compression (Yano et al. 2004a), standard deviations [cf. Eqs. (4.3) and (4.4) found later in the study] are adopted as measures for defining a threshold for activation and deactivation of constant segments. In the case of the wavelet compression, the standard deviation is known to provide an optimized threshold in the sense that, beyond this threshold, the rate of loss of variance by truncation becomes faster than the rate of enhancing compression rates. Qualitatively, a similar behavior is expected for the present activation and deactivation strategy.

In the present study, both a local criterion [Eq. (4.3)] and a global criterion [Eq. (4.4)] are adopted as thresholds. The local criterion is central because new segments must be added by following the ascent of a local plume from a lower layer toward the vertical layer of concern. Existence of a plume is identified by a marked difference of a physical value crossing two segments, which is strong enough in a given vertical layer. However, a global criterion is additionally required to prevent excessive activation of segments because of the existence of disturbances that are insignificant in a global sense, although they are identified to be significant solely based on a local criterion.

Table 1 lists the parameters for activation and deactivation approximately in the order of appearance in the text with the default values.

Table 1.

List of model parameter default values, which are given by categories roughly in the order introduced in the text.

List of model parameter default values, which are given by categories roughly in the order introduced in the text.
List of model parameter default values, which are given by categories roughly in the order introduced in the text.

a. Initialization and basic managements on vertical levels

The basic strategy is to sustain a continuous evolution of a convective system by maintaining a random noise at the lowest vertical layers. For this goal, a full resolution is introduced to the first few vertical layers of the model and then the system is initialized with a random noise in the first vertical layer. These perturbations work as continuous seeds for generating plumes.

More precisely (cf. Fig. 3), the model is initialized with a horizontally homogeneous distribution of Nx constant segments with the model’s finest horizontal resolution Δx = L/Nx from the surface z = 0 to the level z = kmΔz (i.e., for the lowest km layers). Above z = kmΔz, the number of homogeneously distributed constant segments is drastically reduced to Mx (<Nx), with a much larger segment size ΔX = L/Mx (>Δx).

The initial random noise along with an initial high resolution in the lowest km layers allows perturbations to develop at the lowest levels, which work as a source for plumes. To sustain the plume source with time, the high horizontal resolution Δx is always maintained at the lowest kb layers (up to z = kbΔz) throughout the simulations, where kbkm. Activation and deactivation of segments are performed over a vertical range from z = kbΔz to ktΔz (where kt < Nz) as detailed later.

We introduce the uppermost level ktΔz, where new segments may be activated. Above this level, the total number of segments is always kept at Mx as initially assigned. This last procedure, by averaging out an upward smaller-scale vertical flux into a larger scale, aims to introduce numerical dissipation without explicitly introducing a sponge layer at the uppermost layers.

b. Basic principles for activation and deactivation: Threshold measures

More precisely, at a very technical level it is not a segment itself but an interface between segments that is activated and deactivated with time. A new interface subdivides a preexisting segment into two, leading to a finer resolution, and a deactivation of an existing interface merges the existing two segments into one, leading to a lesser resolution. The basic principle is schematically shown in Fig. 4.

Fig. 4.

A schematic for the activation principle. A segment with a significant difference (indicated by a thick vertical bar) at the level k + ½ is extended both upward and downward for Δka layers (Δka = 2 in the schematic, as indicated by a thick dashed line). Deactivation is essentially a backward procedure to the activation with the parameter Δka replaced by Δkd (not shown). See the text for the details.

Fig. 4.

A schematic for the activation principle. A segment with a significant difference (indicated by a thick vertical bar) at the level k + ½ is extended both upward and downward for Δka layers (Δka = 2 in the schematic, as indicated by a thick dashed line). Deactivation is essentially a backward procedure to the activation with the parameter Δka replaced by Δkd (not shown). See the text for the details.

As a criterion for activation and deactivation of segment interfaces, a difference of physical values between two segments crossing an interface is used,

 
formula

where ϕj is a value of a variable ϕ at the jth segment, and lj−1/2 is an integer length of the shortest segment adjacent to this interface given in the unit of the maximum model resolution Δx, thus

 
formula

This qualitatively [Eq. (4.1)] provides an “energy” measure for ϕ.3 As a result, larger segments are preferably activated and smaller segments are preferably deactivated. The shorter of the two segments adjacent to an interface is used for the energy measure so that an activation of a shorter segment is always not favored.

If the difference (4.1) is above a threshold for any of the prognostic variables, a new segment interface is added (activated) above and below the given vertical layer (if an interface does not exist yet) over Δka vertical layers. On the other hand, when differences of “neighboring” interfaces are below thresholds over Δkd layers both above and below, the corresponding segment interface is removed (deactivated). When Δkd = 0, the test is performed only at the given level of the interface in concern.

The thresholds are defined relative to the standard deviation of a given variable. Both local and global measures are introduced. As a local measure, we take the standard deviation for a variable ϕ at a half vertical level k + ½ defined by

 
formula

with a fractional length σj,k+1/2 ≡ (xj+1/2,k+1/2xj−1/2,k+1/2)/L for the jth segment and a horizontal mean ϕk+1/2:

 
formula

Additionally, a standard deviation averaged over the whole depth of the model (total standard deviation),

 
formula

is considered a global measure.

In numerical implementation, deactivations are performed first and then activations follow. As a result, the deactivated interfaces may be immediately activated again at the same time step when the latter condition is less restrictive than the former, but with the physical values homogenized over these reactivated interfaces (cf. section 4e). In this sense, deactivation works, to a certain degree, like a “cleaner” or horizontal filtering of the field. Deactivation and activation are performed for every nd and na time steps, respectively.4

c. Activations

For the activations, the basic principle is to test at a given half level, say, k + ½, if there is an interface with a strong difference relative to the standard deviation defined at a half level immediately above or below, that is, or k − ½. Then, depending on the case, a new interface is added upward or downward (cf. Fig. 4).

Under this strategy, the following test is first applied at every half level, k + ½, from k = kb − 1 to kt − 1, at every interface j for every variable ϕ:

 
formula
 
formula

for the half levels immediately above and below, that is, k′ = k + 1 and k − 1, as long as the level stays within the range kbk′ ≤ kt − 1. Here, γa and γmin are the thresholds for local [Eq. (4.5a)] and global [Eq. (4.5b)] conditions, respectively, relative to the standard deviations defined by (4.3) and (4.4). When both conditions are satisfied at a level k′ + ½ with either k′ = k − 1 or k + 1, the jth interface located at xj−1/2,k+1/2,b is extended either downward or upward for Δka layers depending on k′ = k − 1 or k + 1. For example, when k′ = k + 1, the interface at xj−1/2,k+1/2,b is extended upward from the level k + ½ to .

d. Deactivations

Deactivation is essentially a backward procedure to the activation with the parameter Δka replaced by Δkd. The tests are performed over the 2Δkd + 1 layers. When a criterion is satisfied, the segment at the middle layer is removed.

More specifically, at every half level, k + ½, the following conditions are tested for every prognostic variable:

 
formula

Here, γd is a local deactivation threshold that is set equal to γa in a standard case (cf. Table 1), although sensitivity tests are performed with γdγa (cf. Table 2). The test is performed from k′ ranging from k − Δkd to k + Δkd as long as k′ remains within the range kbk′ ≤ kt − 1. If the test is passed (i.e., either of these two conditions is satisfied over all k′) for three consecutive interfaces at each vertical level, a middle interface indexed by, say, j is classified as a candidate for deactivation.

Table 2.

List of Ayotte test cases. List of Ayotte test cases used in the scatter plots of Fig. 10; there are four major categories shown by different symbols in the figure.

List of Ayotte test cases. List of Ayotte test cases used in the scatter plots of Fig. 10; there are four major categories shown by different symbols in the figure.
List of Ayotte test cases. List of Ayotte test cases used in the scatter plots of Fig. 10; there are four major categories shown by different symbols in the figure.

When Δkd is not zero, the following test is furthermore applied for k′ spanning from k − Δkd to k + Δkd but excluding k:

 
formula

The test is performed over the three interfaces: j′ = j(k + ½) − 1, j(k + ½), and j(k + ½) + 1. Here, the explicit argument, k + ½, is added in j(k + ½) to explicitly indicate the vertical level used for defining the segment index. The test is performed as long as these three interfaces continue vertically. Otherwise, the test is simply skipped. Once these tests are passed, the interface in concern is removed.

e. Assignments of new physical values

New physical values must be assigned to the newly reconfigured segments, in either a case of subdividing a segment into two by activation or a merger of two segments into one by deactivation. When a new interface subdivides a segment into two, the same physical values as those for the original single segment are assigned to the two new segments.

On the other hand, when the jth interface at the level k + ½, say, is removed, then the two original segments merge into one, thus values for all physical variables must be assigned for the merged segment. They are defined by taking an average of values for the two old segments weighted by their relative lengths, that is,

 
formula

with

 
formula

Note that the deactivation is performed only if the number of total interfaces is larger than the greater of 3 or Mx. The original interfaces assigned at the positions xb = iΔX with i = 1, … , M are never deactivated.

5. Numerical validation: Ayotte test case

a. Outline

Ayotte et al. (1996) provide a set of idealized dry large-eddy simulations (LESs) for testing the planetary boundary layer (PBL) parameterizations (see also Hourdin et al. 2002). Here, NAM–SCA is tested as one of the parameterization schemes by the procedure prescribed in Ayotte et al. (1996). Among the cases provided by Ayotte et al. (1996), the case 24F, which is a simulation of a free-convection PBL, is adopted in the present study as the simplest and easiest to test.

For comparisons, Ayotte et al. (1996) provide an experiment with an LES (Ayotte–LES), originating from Moeng (1984). This LES run is also adopted as a reference in the present study. Note that the adopted resolution in NAM–SCA (Δx = 50 m, Δz = 20 m; cf. Table 1) is comparable to the one adopted in the Ayotte–LES (cf. Ayotte et al.’s 1996 Table III). However, it is important to note that the present NAM–SCA, for the sake of simplicity, contains no subgrid-scale parameterization, unlike standard LESs.5

By following the test instruction provided by Ayotte et al. (1996), the domain-mean values of an equilibrium state obtained by running the Ayotte–LES for five large-eddy times, or 5τ, are adopted as the initial conditions. The domain-mean potential temperature at t = 5τ of the Ayotte–LES run is also adopted as a reference profile θr. On top of these domain-averaged values, a random perturbation is added to the potential temperature field with a standard deviation of 0.2 K to the first two layers of the model, following the standard procedure in LES initialization.

The Ayotte–LES is run for another 10τ, and the results averaged over the last 5τ are provided to the model comparison project participants for a comparison with parameterizations. A result after 7τ of a parameterization is prescribed to compare with this averaged data. The present NAM–SCA run is also compared with the Ayotte–LES results in this manner.

In the case considered (24F), the system is forced with a constant surface flux H defined by H/ρrCp = 0.25 K m s−1. Here, Cp is the specific heat of the dry air at constant pressure. In the present model implementation, the surface flux is added as a tendency for the potential temperature in the first two model layers as

 
formula

where hb is the top height of the second layer. The surface flux is redistributed over the first two layers instead of only the first layer so that the lowest-level vertical velocity, evaluated at the interface between the first two layers, can fully feel the development of buoyancy by the surface flux.6

In this test case, the large-eddy time is τ = 515 s, and the height of the capping inversion at t = 5τ is at zi = 1033 m. The potential temperature profile of the initial condition (t = 5τ of the Ayotte–LES run; long dashed curve in Fig. 5a) has a numerically, artificially generated, strong inversion capping, which is gradually eroded with time because of a constant heat flux from the surface.

Fig. 5.

The Ayotte test case. The vertical profiles of the domain-mean potential temperature obtained after 7τ by NAM–SCA with varying horizontal domain sizes. The horizontal domain size increases from left to right as 6.4 (standard case), 12.8, 25.6, and 51.2 km. The horizontal resolution is fixed to Δx = 50 m. Also shown are the initial condition (long dashed) and the result from the Ayotte–LES after 7τ (short dashed). (b) Solid curves as in (a), but for the vertical heat flux. The long dashed curve is the one from Ayotte–LES at t = 7τ. The short dashed curves are the same as the solid curves, except for an estimate (5.2) based on the potential temperature tendency over one large-eddy time. In both panels, the curves are shifted horizontally with a constant (+2 K) increment for displaying the nonstandard cases (L = 12.8, 25.6, and 51.2 km).

Fig. 5.

The Ayotte test case. The vertical profiles of the domain-mean potential temperature obtained after 7τ by NAM–SCA with varying horizontal domain sizes. The horizontal domain size increases from left to right as 6.4 (standard case), 12.8, 25.6, and 51.2 km. The horizontal resolution is fixed to Δx = 50 m. Also shown are the initial condition (long dashed) and the result from the Ayotte–LES after 7τ (short dashed). (b) Solid curves as in (a), but for the vertical heat flux. The long dashed curve is the one from Ayotte–LES at t = 7τ. The short dashed curves are the same as the solid curves, except for an estimate (5.2) based on the potential temperature tendency over one large-eddy time. In both panels, the curves are shifted horizontally with a constant (+2 K) increment for displaying the nonstandard cases (L = 12.8, 25.6, and 51.2 km).

b. Full-resolution case: The domain size dependence

The first case performed is an experiment with a full-resolution Δx found everywhere by setting Mx = Nx = 128, but it otherwise retains the standard parameters given in Table 1. A snapshot of this simulation at t = 7τ is shown in Fig. 6. A few well-organized upward convective plumes are identified in both in the vertical velocity (upper panel) and the potential temperature (lower panel) fields. It is also seen that intrusion of plumes into the inversion locally generates a negative potential temperature anomaly.

Fig. 6.

A snapshot of the Ayotte case NAM–SCA simulation with a full resolution at t = 7τ : anomalies (deviation from the reference state) (top) vertical velocity and (bottom) potential temperature.

Fig. 6.

A snapshot of the Ayotte case NAM–SCA simulation with a full resolution at t = 7τ : anomalies (deviation from the reference state) (top) vertical velocity and (bottom) potential temperature.

The domain-mean θ profile (θ) obtained by the full-resolution case at t = 7τ is shown as the leftmost solid curve in Fig. 5a. The curve is compared with the corresponding profile obtained by the Ayotte–LES, which is shown by the leftmost short dashed curve. Furthermore, Fig. 5b shows the vertical heat flux in a similar format: the long dashed curve shows the result from Ayotte–LES, whereas the leftmost solid curve shows the vertical heat flux wθ, where the prime designates a deviation from the horizontal mean. The total heat flux estimated from a potential temperature tendency over a large-eddy time τ as defined by Eq. (25) of Ayotte et al. (1996), but correcting their sign mistake,

 
formula

is also shown by the leftmost short dashed curve. The estimated total heat flux includes various numerical effects indirectly contributing to the vertical flux.

The first defect of the model to point out may be the tendency of an instantaneous heat flux wθ (the leftmost solid curve in Fig. 5b) to not constitute a straight line over a well-mixed layer, as should be. A closer inspection shows that the domain-mean potential temperature is not quite stationary with time either, even after 7τ. Unsteadiness of heat flux is attributed to a relatively small domain size adopted for the simulation, which can only contain a few plumes at any single moment within the given two-dimensional domain (cf. Moeng et al. 2004), as shown by a snapshot in Fig. 6.

To see the change with increasing numbers of available plumes at a single moment, the full-resolution run is repeated with varying sizes of the horizontal domain (with the fixed horizontal resolution). The results are shown by shifting the solid and the short dashed curves to the right with a constant (+2 K) increment in Fig. 5 for L = 12.8, 25.6, and 51.2 km. We see that the instantaneous vertical heat flux tends to straighten with increasing horizontal domain sizes, associated with increasing numbers of available plumes.

Because of the simplified physics of the current NAM–SCA, the full-resolution result suffers from two additional defects: 1) because of a lack of any surface scheme (not to mention the absence of the Monin–Obukhov condition), the surface layer remains strongly convectively unstable compared to the Ayotte–LES result; 2) above the surface layer, the NAM–SCA profiles are noticeably more stable than the Ayotte–LES within the well-mixed convective boundary layer.

To infer a reason for a stably stratified mixed layer, time evolution of both the domain-mean potential temperature θ and heat flux wθ are plotted in Figs. 7a and 7b, respectively, with an interval of 2τ. Over the first two large-eddy times (long dashed curve), the well-mixed layer is much destabilized (Fig. 7a) with an extensive heat supply to the lower half of the mixed layer, as manifested by a steeper slope in heat flux over the layer 0.4–0.6 km compared to 0.8–1.0-km levels (Fig. 7b). However, over the next four large-eddy times (short dashed and chain dashed curves), more heat is supplied to the upper levels accompanied by a rapid decrease of heat-flux gradient in the upper levels (Fig. 7b), and as a result, the stratification gradually turns from an unstable to a stable profile. At t = 8τ (double chain–dashed curve), the potential temperature profile comes to equilibrium, and from that point the potential temperature increases homogeneously with time in the vertical, maintaining the same weakly stable profile (not shown).

Fig. 7.

The time evolution of the full-resolution case: (a) the domain-mean potential temperature θ and (b) the instantaneous domain-mean heat flux wθ for t = 0 − 8τ with an interval of 2τ with varying curves. Note that θ monotonously increases with time except for a very shallow surface layer, and the same curve as that for θ is used for the heat flux at the same t (in the increasing order of time: solid, long dash, short dash, chain dash, and double chain dash).

Fig. 7.

The time evolution of the full-resolution case: (a) the domain-mean potential temperature θ and (b) the instantaneous domain-mean heat flux wθ for t = 0 − 8τ with an interval of 2τ with varying curves. Note that θ monotonously increases with time except for a very shallow surface layer, and the same curve as that for θ is used for the heat flux at the same t (in the increasing order of time: solid, long dash, short dash, chain dash, and double chain dash).

We believe that a tendency toward a weak stratification is due to very buoyant convective plumes that directly penetrate into the upper half of the mixed layer without much mixing in the lower half of the layer. These penetrating plumes are generated, after an initial transitional destabilization dominated by shallower plumes, because of an unusually strong unstable surface layer.

The plumes can transport heat upward even over a stably stratified layer because of their nonlocal nature of heat transport. In fact, such a nonlocal heat transport is the core of the hot-tower hypothesis by Riehl and Malkus (1958; see also Yano 2009). They recognize that such a nonlocal process must be invoked to explain an up-gradient heat transport against a typical equivalent potential temperature profile in the upper half of the tropical troposphere. A similar process is clearly going on in the present weakly stratified mixed layer, except that the stratification is much weaker and the layer is much more vertically extended.

This argument is also supported by examining the contribution of heat transport due to mass-flux parameterizations in convective boundary layers. It is typically seen that a heat-flux gradient due to mass flux monotonously decreases with height from positive just above the surface layer to negative above a middle level (e.g., Fig. 7d of Siebesma et al. 2007; Fig. 3c of Pergaud et al. 2009), as also seen in heat flux at equilibrium in the present case (cf. Fig. 7b). As a result, if the eddy-diffusive transport were simply removed, heat would be preferentially transported to an upper half of the mixed layer, resulting in a stable stratification.

c. Compression tests

An example of snapshots from a run with a relatively strong compression is shown in Fig. 8. The parameters used for this case are the default parameters (cf. Table 1), γa = γd = 1.0, na = nd = 10, etc. In the upper frame, where the vertical velocity field is shown, positions of the segment interfaces are also indicated by vertical bars.

Fig. 8.

As in Fig. 6, but under a strong compression with thresholds γa = γd = 1 at (a) t = 8τ and (b) t = 9τ. (top) Interfaces between the segments are also indicated by vertical bars.

Fig. 8.

As in Fig. 6, but under a strong compression with thresholds γa = γd = 1 at (a) t = 8τ and (b) t = 9τ. (top) Interfaces between the segments are also indicated by vertical bars.

At t = 8τ (Fig. 8a), a plume toward the left edge of a domain is well represented by a dense distribution of segments, whereas the quiescent region outside is represented by a much scarcer distribution of segments. One large-eddy time later (Fig. 8b), the main activity of plumes shifts rightward, and so does the region with dense distribution of segments by time-dependent activation and deactivation of segments.

The modifications of the vertical profiles with increasing compression rates are shown in Fig. 9 for potential temperature and heat flux in the same format as in Fig. 5. That is, cases with increasing thresholds— γa = γd = 0.2, 0.5, 1.0, 2.0— are shown by shifting the curves from the left to the right with a constant (+2 K) increment. The relative total number of segments compared to the full-resolution case, or the compression rate, of the system is 0.490, 0.284, 0.144, and 8.74 × 10−2, respectively. Compared to a degree of decreasing compression rates, deterioration of the vertical profiles both for θ and wθ compared to the full version (the leftmost curve in Fig. 5) appears to be relatively slow, although the deterioration of the heat flux stands out more than that for the θ profile.7

Fig. 9.

As in Fig. 5, but showing the change with decreasing compression rates. Four cases are shown by shifting the curves from the left to the right with a constant (+2 K) increment, γa = γd = 0.2, 0.5, 1.0, and 2.0. Note that the curves for γa = γd = 0.2 represent true values.

Fig. 9.

As in Fig. 5, but showing the change with decreasing compression rates. Four cases are shown by shifting the curves from the left to the right with a constant (+2 K) increment, γa = γd = 0.2, 0.5, 1.0, and 2.0. Note that the curves for γa = γd = 0.2 represent true values.

To better compare the runs, the relative error is defined as the ratio between the RMS error of a run and the RMS difference from the full resolution run at t = 7τ,

 
formula

Here, θ(z, t) refers to a domain-mean vertical profile of the potential temperature for a given run at time t, whereas θf(z, t) is the same but for the full-resolution case, shown by the leftmost curve in Fig. 5.

A relative CPU is also defined by the CPU of a given run divided by the CPU for the full-resolution run. Each run is performed for 25 large-eddy times. All the calculations are performed with a local Xeon machine (5 GHz). The CPU time for the full-resolution run is 1232.61 s.

Extensive sensitivity tests are performed by changing the model parameters in various ways. Tests are categorized into the four sets as listed in Table 2. The results are shown by scatterplots against the compression rate both for the relative RMS error and the relative CPU in Figs. 10a,b, respectively, in which each category is shown by different symbols.

Fig. 10.

Compression tests of NAM–SCA for the Ayotte case: the compression rate (relative number of segments to a full-resolution case) is scatter plotted against (a) the relative RMS error [Eq. (5.3)] and (b) the relative CPU. The list of the plotted cases is given in Table 2.

Fig. 10.

Compression tests of NAM–SCA for the Ayotte case: the compression rate (relative number of segments to a full-resolution case) is scatter plotted against (a) the relative RMS error [Eq. (5.3)] and (b) the relative CPU. The list of the plotted cases is given in Table 2.

Figure 10a shows that an increase of error with decreasing (stronger) compression rates is relatively slow above the compression rate 0.2–0.3, and only below this level does it begin to disperse enormously, but with the lowest error obtained still remaining comparable to those obtained under weaker compressions. This result shows that NAM–SCA can run with a relatively strong compression rate (down to ca. 0.2) with a relatively small relative error (ca. 0.2). This may be compared with a relatively slow increase of error under wavelet compressions as shown in Fig. 11c and 12c of Yano et al. (2004a). Qualitatively, the SCA works in a similar manner as wavelet compression, by conceptually following the spirit of multiresolution approach (cf. Mallat 1998).

Figure 10b shows that overall, the CPU time decreases linearly with decreasing (stronger) compression rates; the spread is relatively weak. It implies that the CPU overhead for actually performing activation and deactivation does not sensibly depend on the choice of the parameters, especially the frequency of activation and deactivation na and nd. The lowest bound of the scatter is approximately extrapolated to the relative CPU of 0.2, when the compression rate approaches zero. The nonvanishing CPU in this limit provides a measure of the numerical burden for solving the Poisson problem with full resolution in the current algorithm regardless of the compression rate of the other parts.

Extensive investigations are performed for identifying a particular combination of the parameters that can compress the system extensively while keeping a relatively small error. However, none is identified as a particularly strong choice, but stronger compressions always lead to larger errors regardless of the parameter choice, as indicated by Fig. 10a. Exceptionally fortunate cases in scatter provide good results (high compressions with low relative errors) mostly due to much weaker compression rates in earlier parts of the integrations. Recall that both the relative error and the compression rate are evaluated at a particular moment t = 7τ.

6. Discussion and conclusions

The core of mass-flux convective parameterization may be considered as subdividing a gridbox domain into various constant segments corresponding to an environment and convective elements, such as convective towers. The goal of the present paper has been to present a model that simply introduces such a subdivision of a model domain into a nonhydrostatic cloud-resolving or a large-eddy simulation model. The model constitutes a prototype for mass-flux convective parameterization, which is further developed by introducing additional assumptions and hypotheses.

From a numerical point of view, the constructed model is based on a finite-volume approach (cf. Godunov 1959; section 5.6.1, Durran 1999; LeVeque 2002). However, a major philosophical departure from the standard finite-volume methods is that the horizontal extent of a volume element (constant segments) is not necessarily small at all, but is intended to represent a convective plume, or even an entire environmental segment in the same spirit as the parameterization of convective mass-flux. As a result, the distribution of finite volumes could be much more inhomogeneous than in standard approaches. A degree of strong inhomogeneity is in keeping with the spirit of the multiresolution approach (cf. Fournier et al. 2005).

These finite-volume elements are added and removed (activated and deactivated) with time by following the birth and death of convective elements. The procedure may be considered a type of adaptive mesh refinement (AMR; e.g., Berger and Colella 1989; Bell et al. 1994). However, a major philosophical departure from standard AMRs must also be emphasized.

AMR is based on the idea of refinement, where the major goal is to describe local “singularities” associated with shocks and discontinuities (e.g., Berger and Colella 1989; Grauer et al. 1998; Rosenberg et al. 2006) by introducing a local refinement of grid or meshes. To make this procedure systematic, a hierarchical structure for mesh refinement is often developed so that a higher refinement of a singular structure is achieved simply by moving to a higher order of the hierarchy. The hierarchical framework also enables us to mathematically prove the convergence of a solution.

On the other hand, the present SCA approach is based on the idea of compression under the spirit of the multiresolution analysis, where the major goal is to reduce the total number of constant segments to achieve a higher numerical efficiency. In the traditional mass-flux framework, a main portion of the gridbox domain is considered to be occupied by a horizontally homogeneous environment. Under the same spirit, the environment (nonplume region) is intended to be represented by as few segments as possible in the present approach (cf. Fig. 8). Note that under this approach the system always has a fewer number of segments than a full horizontal resolution case, and it asymptotically converges to a full-resolution case as the compression is weakened.

The development of the nonhydrostatic anelastic model with segmentally constant approximations (NAM–SCA) has enabled us to successfully simulate the evolution of convective plumes within a boundary layer with a relatively limited number of segments; however, the number of segments required for representing a single plume turns out to be closer to a few rather than a single segment, as a quick inspection of Fig. 8 shows. Thus, the present adaptive refinement scheme suggests the difficulty of representing a single plume by a single segment as formally envisioned in the mass-flux formulation.

A moist version with simple microphysics has already been developed. Preliminary results for tests with deep moist convection appear to be more promising than the dry boundary layer case reported here.

Further improvements of the model are also anticipated. At the most technical level, the Poisson problem is solved on a regular grid with the present model. This is a major task to take on, because to the best of our knowledge, there is no ready-made package that can be adopted in the SCA configuration. Currently, possibilities are under consideration for reformulating the Poisson problem on an irregular grid under a conjugate gradient approach (cf. Smolarkiewicz and Margolin 1994; Skamarock et al. 1997). Such a procedure should scale the numerical efficiency of the Poisson solver linearly with the compression.

The numerical efficiency of the model as a whole may also be improved by further enhancing the compressibility of the model. The “renormalization” approach (Yano et al. 2003) would be one of the possibilities to investigate for this purpose.

Subgrid-scale eddy parameterization can also be incorporated into the model. Effects of fluctuations within a constant segment can be included by taking an analogous procedure to the standard moment expansion approach (cf. Mellor and Yamada 1974). It is especially important to note that the inclusion of vertical eddy diffusion is straightforward because SCA is applied only in the horizontal direction.

The most important implication suggested by the present study is that explicit modeling and parameterizations are not two distinctively separated approaches as is often considered, but they only constitute two extreme limits of available approaches for investigating subgrid-scale processes (cf. Yano et al. 2005). The present study suggests that the gap may be filled by reinterpreting the mass-flux formulation as a type of finite-volume approach, coined SCA. By applying SCA to a different degree to a NAM (LES and CRM), a hierarchy of models that fill the gap between the explicit model and parameterizations may be developed. This is another possibility currently under investigation.

Acknowledgments

The first author expresses particular gratitude to A. Fournier for the discussion during a long process searching for a best adaptive approach. The first author has also benefited from discussion with D. Rosenberg, A. St-Cyr, and W. Skamarock on adaptive approaches, and C.-H. Moeng on the interpretations of the boundary layer simulations. Discussion with P. Bechtold, P. Bernardet, J.-F. Gueremy, F. Guichard, J.-M. Piriou, J.-L. Redelsperger, and J. Stein on various occasions is also acknowledged. Thanks are also owed to C. Rio, P. Marque, and F. Hourdin for arranging the original LES simulation data from Ayotte et al. (1996). The text was carefully read by M. Freer.

REFERENCES

REFERENCES
Arakawa
,
A.
, and
W. H.
Schubert
,
1974
:
Interaction of a cumulus cloud ensemble with the large–scale environment, Part I.
J. Atmos. Sci.
,
31
,
674
701
.
Asai
,
T.
, and
A.
Kasahara
,
1967
:
A theoretical study of the compensating downdraft motions associated with cumulus clouds.
J. Atmos. Sci.
,
24
,
487
496
.
Ayotte
,
K. W.
, and
Coauthors
,
1996
:
An evaluation of neutral and convective planetary boundary-layer parameterizations relative to large eddy simulations.
Bound.-Layer Meteor.
,
79
,
131
175
.
Baker
,
T. J.
,
1997
:
Mesh adaptation strategies for problems in fluid dynamics.
Finite Elem. Anal. Des.
,
25
,
243
273
.
Bell
,
J.
,
M.
Berger
,
J.
Saltzman
, and
M.
Welcome
,
1994
:
Three-dimensional adaptive mesh refinement for hyperbolic conservation laws.
SIAM J. Sci. Comput.
,
15
,
127
138
.
Berger
,
M. J.
, and
P.
Colella
,
1989
:
Local adaptive mesh refinement for shock hydrodynamics.
J. Comput. Phys.
,
82
,
64
84
.
Bernardet
,
P.
,
1995
:
The pressure term in the anelastic model: A systematic elliptic solver for an Arakawa C grid in generalized coordinate.
Mon. Wea. Rev.
,
123
,
2474
2490
.
Dietachmayer
,
G. S.
, and
K. K.
Droegemeier
,
1992
:
Application of continuous dynamic grid adaption techniques to meteorological modeling. Part I: Basic formulation and accuracy.
Mon. Wea. Rev.
,
120
,
1675
1706
.
Durran
,
D. R.
,
1991
:
The third-order Adams–Bashforth method: An attractive alternative to leapfrog time differencing.
Mon. Wea. Rev.
,
119
,
702
720
.
Durran
,
D. R.
,
1999
:
Numerical Methods for Wave Equations in Geophysical Fluid Dynamics.
Springer, 465 pp
.
Fournier
,
A.
,
G.
Beylkin
, and
V.
Cheruvu
,
2005
:
Multiresolution adaptive space refinement in geophysical fluid dynamics simulation.
Adaptive Mesh Refinement—Theory and Application, T. Plewa, T. Linde, and V. G. Weirs, Eds., Lecture Notes in Computational Sciences and Engineering, Vol. 41, Springer-Verlag, 161–170
.
Godunov
,
S. K.
,
1959
:
A difference scheme for numerical computation of discontinuous solutions of equations in fluid mechanics.
Math. Sb.
,
47
,
271
306
.
Grabowski
,
W. W.
,
2004
:
An improved framework for superparameterization.
J. Atmos. Sci.
,
61
,
1940
1952
.
Grabowski
,
W. W.
, and
P. K.
Smolarkiewicz
,
1999
:
CRCP: A cloud resolving convection parameterization for modeling the tropical convective atmosphere.
Physica D
,
133
,
171
178
.
Grabowski
,
W. W.
,
X.
Wu
, and
M. W.
Moncrieff
,
1996
:
Cloud-resolving modeling of tropical cloud systems during phase III of GATE. Part I: Two-dimensional experiments.
J. Atmos. Sci.
,
53
,
3684
3709
.
Grauer
,
R.
,
C.
Marliani
, and
K.
Germaschewski
,
1998
:
Adaptive mesh refinement for singular solutions of the incompressible Euler equations.
Phys. Rev. Lett.
,
80
,
4177
4180
.
Haltiner
,
G. J.
, and
R. T.
Williams
,
1980
:
Numerical Prediction and Dynamic Meteorology.
John Wiley & Sons, 477 pp
.
Hourdin
,
F.
,
F.
Couvreux
, and
L.
Menut
,
2002
:
Parameterization of the dry convective boundary layer based on a mass-flux representation of thermals.
J. Atmos. Sci.
,
59
,
1105
1123
.
Hubbard
,
M. E.
, and
N.
Nikiforakis
,
2003
:
A three-dimensional, adaptive, Godunov-type model for global atmospheric flows.
Mon. Wea. Rev.
,
131
,
1848
1864
.
Jablonowski
,
C.
,
M.
Herzog
,
J. E.
Penner
,
R. C.
Oehmke
,
Q. F.
Stout
,
B.
van Leer
, and
K. G.
Powell
,
2006
:
Block-structured adaptive grids on the sphere: Advection experiments.
Mon. Wea. Rev.
,
134
,
3691
3713
.
Klemp
,
J. B.
, and
D. R.
Durran
,
1983
:
An upper boundary condition permitting internal gravity wave radiation in numerical mesoscale models.
Mon. Wea. Rev.
,
111
,
430
444
.
LeVeque
,
R. J.
,
2002
:
Finite-Volume Methods for Hyperbolic Problems.
Cambridge University Press, 578 pp
.
Lilly
,
D. K.
,
1965
:
On the computational stability of numerical solutions of time–dependent non–linear geophysical fluid dynamics problems.
Mon. Wea. Rev.
,
93
,
11
26
.
Lilly
,
D. K.
,
1983
:
Editor’s note.
Mesoscale Meteorology—Theories, Observations and Models, D. K. Lilly and T. Gal-Chen, Eds., Reidel, 447–450
.
Mallat
,
S.
,
1998
:
A Wavelet Tour of Signal Processing.
2nd ed. Academic Press, 637 pp
.
Margolin
,
L. G.
,
P. K.
Smolarkiewicz
, and
Z.
Sorbjan
,
1999
:
Large-eddy simulation of convective boundary layers using nonoscillatory differencing.
Physica D
,
133
,
390
397
.
Mellor
,
G. L.
, and
T.
Yamada
,
1974
:
A hierarchy of turbulence closure models for planetary boundary layers.
J. Atmos. Sci.
,
31
,
1791
1806
.
Mesinger
,
F.
, and
A.
Arakawa
,
1976
:
Numerical methods used in atmospheric models.
GARP Publication Series, No. 14, WMO/ICSU Joint Organizing Committee, 64 pp
.
Moeng
,
C-H.
,
1984
:
A large-eddy-simulation model for the study of planetary boundary–layer turbulence.
J. Atmos. Sci.
,
41
,
2052
2062
.
Moeng
,
C-H.
,
J. C.
McWillaims
,
R.
Torunno
, and
P. P.
Sullivan
,
2004
:
Investigating 2D modeling of atmospheric convection in the PBL.
J. Atmos. Sci.
,
61
,
889
903
.
Morton
,
B. R.
,
1997
:
Discrete dry convective entities I: Review.
The Physics and Parameterization of Moist Atmospheric Convection, R. K. Smith, Ed., Kluwer Academic, 143–173
.
Morton
,
B. R.
,
G.
Taylor
, and
J. S.
Turner
,
1956
:
Turbulent gravitational convection from maintained and instantaneous sources.
Proc. Roy. Meteor. Soc.
,
A234
,
1
33
.
Ooyama
,
K.
,
1971
:
A theory on parameterization of cumulus convection.
J. Meteor. Soc. Japan
,
26
,
3
40
.
Pergaud
,
J.
,
V.
Masson
,
S.
Malardel
, and
F.
Couvreux
,
2009
:
A parameterization of dry thermals and shallow cumuli for mesoscale numerical weather prediction.
Bound.-Layer Meteor.
,
132
,
83
101
.
Randall
,
D. A.
,
M.
Khairoutdinov
,
A.
Arakawa
, and
W.
Grabowski
,
2003
:
Breaking the cloud parameterization deadlock.
Bull. Amer. Meteor. Soc.
,
84
,
1547
1564
.
Riehl
,
H.
, and
J. S.
Malkus
,
1958
:
On the heat balance in the equatorial trough zone.
Geophysica
,
6
,
503
538
.
Rosenberg
,
D.
,
A.
Fournier
,
P.
Fischer
, and
A.
Pouquet
,
2006
:
Geophysical–astrophysical spectral-element adaptive refinement (GASpAR): Object-oriented h-adaptive fluid dynamics simulation.
J. Comput. Phys.
,
215
,
59
80
.
Scorer
,
R. S.
,
1957
:
Experiments on convection of isolated masses of buoyancy fluid.
J. Fluid Mech.
,
2
,
583
594
.
Siebesma
,
A. P.
,
P. M. M.
Soares
, and
J.
Teixeira
,
2007
:
A combined eddy-diffusivity mass-flux approach for the convective boundary layer.
J. Atmos. Sci.
,
64
,
1230
1248
.
Simpson
,
J.
,
1983a
:
Cumulus clouds: Early aircraft observations and entrainment hypotheses.
Mesoscale Meteorology—Theories, Observations and Models, D. K. Lilly and T. Gal-Chen, Eds., Reidel, 355–373
.
Simpson
,
J.
,
1983b
:
Cumulus clouds: Interactions between laboratory experiments and observations as foundations for models.
Mesoscale Meteorology—Theories, Observations and Models, D. K. Lilly and T. Gal-Chen, Eds., Reidel, 399–412
.
Skamarock
,
W. C.
, and
J. B.
Klemp
,
1993
:
Adaptive grid refinement for two-dimensional and three-dimensional nonhydrostatic atmospheric flow.
Mon. Wea. Rev.
,
121
,
788
804
.
Skamarock
,
W. C.
,
P. K.
Smolarkiewicz
, and
J. B.
Klemp
,
1997
:
Preconditioned conjugate–residual solvers for Helmholtz equations in nonhydrostatic models.
Mon. Wea. Rev.
,
125
,
587
599
.
Smolarkiewicz
,
P. K.
, and
L. G.
Margolin
,
1994
:
Variational solver for elliptic problems in atmospheric flows.
Appl. Math. Comp. Sci.
,
4
,
527
551
.
St-Cyr
,
A.
,
C.
Jablonowski
,
J. M.
Dennis
,
H. M.
Tufo
, and
S. J.
Thomas
,
2008
:
A comparison of two shallow-water models with nonconforming adaptive grids.
Mon. Wea. Rev.
,
136
,
1898
1922
.
Stommel
,
H.
,
1947
:
Entrainment of air into a cumulus cloud.
J. Meteor.
,
4
,
91
94
.
Stommel
,
H.
,
1951
:
Entrainment of air into a cumulus cloud II.
J. Meteor.
,
8
,
127
129
.
Turner
,
J. S.
,
1962
:
The ‘starting plume’ in neutral surroundings.
J. Fluid Mech.
,
13
,
356
368
.
Turner
,
J. S.
,
1986
:
Turbulent entrainment: The development of the entrainment assumption, and its application to geophysical flows.
J. Fluid Mech.
,
173
,
431
471
.
Yano
,
J-I.
,
2009
:
Deep-convective vertical transport: What is mass flux?
Atmos. Chem. Phys. Discuss.
,
9
,
3535
3553
.
Yano
,
J-I.
,
P.
Bechtold
, and
J-L.
Redelsperger
,
2003
:
“Renormalization” approach for subgrid-scale representations.
J. Atmos. Sci.
,
60
,
2029
2038
.
Yano
,
J-I.
,
P.
Bechtold
,
J-L.
Redelsperger
, and
F.
Guichard
,
2004a
:
Wavelet-compressed representation of deep moist convection.
Mon. Wea. Rev.
,
132
,
1472
1486
.
Yano
,
J-I.
,
R.
Blender
,
C.
Zhang
, and
K.
Fraedrich
,
2004b
:
1/f noise and pulse–like events in the tropical atmospheric surface variabilities.
Quart. J. Roy. Meteor. Soc.
,
130
,
1697
1721
.
Yano
,
J-I.
,
F.
Guichard
,
J-P.
Lafore
,
J-L.
Redelsperger
, and
P.
Bechtold
,
2004c
:
Estimations of mass fluxes for cumulus parameterizations from high-resolution spatial data.
J. Atmos. Sci.
,
61
,
829
842
.
Yano
,
J-I.
,
J-L.
Redelsperger
,
F.
Guichard
, and
P.
Bechtold
,
2005
:
Mode decomposition as a methodology for developing convective-scale representations in global models.
Quart. J. Roy. Meteor. Soc.
,
131
,
2313
2336
.

APPENDIX

Numerical Procedures

A basic FORTRAN code used for NAM–SCA is available online from an anonymous FTP site (available online at ftp://cnrm-ftp.meteo.fr/pub-moana/yano/crm-sca/with_activation_ayotte). An extensive user’s guide (available online at ftp://cnrm-ftp.meteo.fr/pub-moana/yano/crm-sca/with_activation_ayotte/ms.pdf) is also included.

a. Overview

From a numerical algorithmic point of view, the NAM–SCA system consists of the three major components. First, by applying SCA, the horizontal direction of the system is numerically already represented by a finite volume in Eq. (3.4). Second, in the vertical direction, Eq. (3.4) is defined in a fully continuous manner that must be discretized for a numerical purpose. On the other hand, the Poisson problem (2.3) is defined in a fully continuous manner in both the horizontal and vertical directions. Thus, discretization must be introduced in two spatial directions to numerically solve this problem. Finally, the horizontal velocity is evaluated by Eq. (3.6) analytically once the vertical-velocity divergence is given.

b. Discretizations

The system is discretized in the vertical direction for solving both the prognostic Eq. (3.4) as well as the Poisson problem (2.3). For this purpose, staggered grids are introduced in the same manner as in Mesinger and Arakawa’s (1976) C grid with a homogeneous vertical resolution Δz. Use of a finite volume also in the vertical direction is feasible, but the choice would make the numerical algorithm slightly more involved.

Under this discretization, the interfaces between the vertical layers are defined by

 
formula

for k = 0, … , Nz with Nz = Hz. The vertical velocity w is evaluated at these interfaces (full levels). In solving the Poisson problem (2.3), we assign the values of the other physical variables at half levels defined by

 
formula

for k = 1, … , Nz.

In solving the Poisson problem (2.3), the system is further discretized into a horizontal direction by assuming a homogeneously distributed grid with a constant resolution Δx over the whole length L of the domain. Under the given discretization, all the physical variables are defined at every grid point based on Eqs. (3.1) and (3.6a). The discretization enables us to evaluate the source term [the right-hand side of Eq. (2.3)] by applying standard centered finite differences in two spatial directions. The source term is transformed into the Fourier space in the horizontal direction, and the vertical derivative in the Laplacian is discretized according to the defined staggered vertical grids, leading to a tridiagonal form of the matrix for the Laplacian. This matrix can be inverted analytically (i.e., with a finite number of numerical operations), and then a Fourier-transformed pressure perturbation p is obtained.

The time integral is performed by a mixture of the two schemes. The Euler time step (cf. Haltiner and Williams 1980) is adopted both for the vertical and horizontal (lateral mixing) advections, whereas the second-order Adams–Bashforth scheme (Lilly 1965; Durran 1991) is adopted for the physical tendency (tendency due to the forcing F). Replacement of the integration scheme for the first two processes by the second-order Adams–Bashforth scheme was found not to improve the numerical efficiency of the calculations.

c. Identification of segment interfaces

To identify the vertical continuation of the segment interfaces, their positions are indexed in two complementary manners at a given vertical level zk, where k is a half integer at a half level—the position of the segment itself, given by xj−1/2,b(zk) for a jth interface,8 and the segment index j for a given discretized horizontal position, xi = iΔx, say jb(xi, zk). Note that not all of the horizontal positions xi are occupied by a segment interface. In case of an absence, we set jb(xi, zk) = −j, with j corresponding to the index of the closest interface to the left. The index jb(xi, zk) is used for identifying the vertical extent of a segment interface given at x = xi. This indexing system enables us to activate and deactivate segments with time without constructing a hierarchical structure.

d. Finite difference in the vertical direction

We integrate Eq. (3.4) in time for the vertical velocity w and the potential temperature θ, respectively, at full and half levels zk and zk+1/2. In these computations, for example, the vertical flux divergence for the jth segment at the vertical level k + ½ is evaluated by

 
formula

where the subscript k, for example, j(k + ½), indicates a value at the vertical level k averaged over a range defined by the jth segment at the vertical level k + ½. The flux is evaluated by an upstream formula ()j(k+1/2),k = wj(k+1/2),kϕj(k+1/2),k* in which the vertical level k* is set k* = k − ½ if wj(k+1/2),k > 0, and k* = k + ½ if wj(k),k < 0.

An evaluation of ϕj(k),k* is performed by

 
formula

where σj(k+1/2),j′(k*) is a fractional contribution of the j′th segment at the vertical level k* given by

 
formula

The pressure gradient is evaluated in a similar manner, but over a fixed segment at a height zk in concern,

 
formula
e. Evaluation levels and vertical interpolations

As already discussed in appendix section b above, the vertical velocity is defined at the full vertical levels, and the potential temperature and the horizontal velocity are defined at the half vertical levels, respectively. For this reason, both the physical tendency [tendency due to the forcing F in Eq. (3.4)] and the vertical advection are evaluated at the full and the half levels, respectively, for the vertical velocity and the potential temperature.

On the other hand, the lateral exchange [the second term in the left-hand side of Eq. (3.4)] is evaluated at the half levels as a layer between two neighboring full levels constitutes a “physical” layer. For this purpose, right before this calculation is performed, the vertical velocity is linearly interpolated onto the half levels. Finally, both prognostic variables are interpolated onto the full vertical levels after adding a tendency of lateral mixing.

Footnotes

Corresponding author address: Jun-Ichi Yano, CNRM, Météo-France, 42 Ave. Coriolis, 31057 Toulouse CEDEX, France. Email: yano@cnrm.meteo.fr

1

We call them plumes rather than hot towers, because they are dry by design.

2

Full arguments of a dependent variable are explicitly shown in this section for clarity. These arguments are omitted in the remaining part of the paper for an economy of the presentation.

3

To see this assume for simplicity that ϕ varies |ϕjϕj−1| over a distance lj−1/2. Then a contribution of the variance 〈ϕ2j (i.e., an energy measure) over this distance is estimated as 〈ϕ2jlj−1/2(ϕjϕj−1)2. The root mean square of this formula gives Eq. (4.1). See also discussion in section 3d of Yano et al. (2004b).

4

A possibility of ndna is considered to evaluate the relative importance of deactivation and activation in accomplishing an efficient compression of the system. As will be later shown by Fig. 10b, it turns out that such differentiation does not sensibly contribute to a computational efficiency.

5

Most notably, the present code does not satisfy the Monin–Obukhov condition at the surface as recommended by Ayotte et al. (1996). Though this modification is very straight forward, the authors rather decide to present the behavior of NAM–SCA in its simplest implementation.

6

Note that, thus, this is the minimum value possible for a parameter hb. The choice of this parameter is comparable to that by Margolin et al. (1999). They assume a subgrid heat flux exponentially decreases upward with a scale height of 25 m.

7

Recall that the estimate of the heat flux by Eq. (5.2) is based on a tendency over τ, thus a problem of statistical fluctuation is avoided.

8

An interface indexed as j − ½ is called the jth interface for convenience.