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.
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.
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 ξ4 ≤ x ≤ ξ6 and ζ3 ≤ z ≤ ζ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 θ,
and the mass continuity,
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
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
with a segment-mean value ϕj(z, t) defined by
(cf. Fig. 2) and an indicator Ij(x, z, t) for the jth segment defined by
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,
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
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 (wϕ)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
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,
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 − ½,
A chain application of Eq. (3.6b) defines the horizontal velocity at every segment interface except for the domain mean u left undetermined.
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
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.
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.
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 kb ≤ km. 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.
As a criterion for activation and deactivation of segment interfaces, a difference of physical values between two segments crossing an interface is used,
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
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
with a fractional length σj,k+1/2 ≡ (xj+1/2,k+1/2 − xj−1/2,k+1/2)/L for the jth segment and a horizontal mean ϕk+1/2:
Additionally, a standard deviation averaged over the whole depth of the model (total standard deviation),
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
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 ϕ:
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 kb ≤ k′ ≤ 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 .
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:
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 kb ≤ k′ ≤ 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.
When Δkd is not zero, the following test is furthermore applied for k′ spanning from k − Δkd to k + Δkd but excluding k:
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,
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
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
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.
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.
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,
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).
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.
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
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τ,
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.
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.
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.
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.
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.
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
for k = 0, … , Nz with Nz = H/Δz. 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
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
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 (wϕ)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
where σj(k+1/2),j′(k*) is a fractional contribution of the j′th segment at the vertical level k* given by
The pressure gradient is evaluated in a similar manner, but over a fixed segment at a height zk in concern,
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.
Corresponding author address: Jun-Ichi Yano, CNRM, Météo-France, 42 Ave. Coriolis, 31057 Toulouse CEDEX, France. Email: firstname.lastname@example.org
We call them plumes rather than hot towers, because they are dry by design.
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.
To see this assume for simplicity that ϕ varies |ϕj − ϕj−1| over a distance lj−1/2. Then a contribution of the variance 〈ϕ2〉j (i.e., an energy measure) over this distance is estimated as 〈ϕ2〉j ∼ lj−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).
A possibility of nd ≠ na 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.
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.
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.
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.
An interface indexed as j − ½ is called the jth interface for convenience.