A semi-Lagrangian advection scheme is developed for falling hydrometeors in hopes of replacing the conventional Eulerian scheme that has been widely used in the cloud microphysics scheme of numerical atmospheric models. This semi-Lagrangian scheme uses a forward advection method to determine the advection path with or without iteration, and advected mass in a two-time-level algorithm with mass conservation. Monotonicity is considered in mass-conserving interpolation between Lagrangian grids and model Eulerian grids, thus making it a positive definite advection scheme. For mass-conserving interpolation between the two grid systems, the piecewise constant method (PCM), piecewise linear method (PLM), and piecewise parabolic method (PPM) are proposed. The falling velocity at the bottom cell edge is modified to avoid unphysical deformation by scanning from the top layer to the bottom of the model, which enables the use of a large time step with reasonable accuracy. The scheme is implemented and tested in the Weather Research and Forecasting (WRF) Single-Moment 3-Class Microphysics Scheme (WSM3).
In a theoretical test bed with constant terminal velocity, the proposed semi-Lagrangian algorithm shows that the higher-order interpolation scheme produces less diffusive features at maximal precipitation. Results from another idealized test bed with mass-weighted terminal velocity demonstrate that the accuracy of the proposed scheme is still satisfactory even with a time step of 120 s when the mean terminal velocity averaged at the departure and arrival points is employed. A two-dimensional (2D) squall-line test using the WSM3 scheme shows that the control run with the Eulerian advection scheme and the semi-Lagrangian run with the PCM method reveal similar results, whereas behaviors using the PLM and PPM are similar with higher-resolution features, such as mammatus-like clouds.
Most numerical weather prediction (NWP) models employ a bulk-type cloud microphysics scheme to account for grid-resolvable precipitation processes, based on the study of Lin et al. (1983). As compared to the continuous efforts given to the development of advanced microphysical processes in the cloud modeling area (e.g., Grabowski 2006), little effort has been made toward improving the accuracy of computing sedimentation processes for precipitating particles in the cloud modeling area. To our knowledge, for computing the sedimentation of precipitation, almost all cloud schemes in NWP models use a subtime step so that falling particles do not cross the vertical grid within the model time step. Aside from uncertainties in numerical accuracy, this procedure is much more time consuming as the model resolution is lower. Thus, in a low-resolution general circulation model (GCM) with a large time step, the efficiency of the Eulerian method becomes a serious computing burden, even though the simulated climatology is reasonable (e.g., Shimpo et al. 2008).
The immediate remedy to this problem is to replace Eulerian advection with semi-Lagrangian (SL) advection [see Staniforth and Cote (1991) or Smolarkiewicz and Pudykiewicz (1992) for a review of the SL method] with a large time step. Pellerin et al. (1995) and Grabowski and Smolarkiewicz (1996) used SL advection for a cloud model without considering the mass conservation. Kato (1995) introduced a box-Lagrangian scheme for raindrop advection with mass conservation; however, the box-Lagrangian scheme as a first order for interpolation produces a diffusive result. There are several higher-order interpolation schemes that have been used in SL advection considering mass conservation, such as those of Laprise and Plante (1995), Lin and Rood (1996), Zurroukat et al. (2002), Kaas (2008), and Juang (2008). These mass conserving SL schemes have been applied to horizontal advection for traces; however, most cloud schemes in NWP and GCMs employ the classic Eulerian advection scheme for the sedimentation of falling hydrometeors. For example, none of the 10 cloud microphysics schemes in the Weather Research and Forecasting (WRF) model (Skamarock et al. 2008), as of August 2009, employs this advanced numerical technique for handling sedimentation.
This study proposes an SL scheme with a higher-order interpolation in order to improve the numerical accuracy and computational efficiency of calculating sedimentation of falling precipitation drops. The WRF Single-Moment Microphysics (WSMMPs) 3-class (WSM3; Hong et al. 2004) scheme is the target model for this SL scheme. The details of the proposed SL advection scheme are illustrated in section 2. The fundamental characteristics of the proposed approach are evaluated in a theoretical test bed, used in past literature [i.e., Xiao et al. (2003), which is based on Yabe et al. (2001)], in section 3. In section 4, the proposed algorithm is further evaluated and compared to the current formula for sedimentation in the WSM3 scheme in 1D and 2D frameworks within the WRF model. Concluding remarks are given in the last section.
2. The proposed forward semi-Lagrangian (SL) scheme
In this section, we expand upon the ideas of Juang (2007, 2008) to introduce a one-step-forward SL scheme, and to illustrate a modification necessary to take care of a leading shock wave to avoid negative cell advection. We also propose an optional modification to yield mean path advection similar to the traditional SL method.
a. Forward semi-Lagrangian advection with mass conservation and positive definiteness
The sedimentation for q, the mixing ratio of precipitation such as rain, ice, snow, or graupel, can be written in flux form as
where all variables are standard meteorological expressions, ρ is the density of air, and VT is the corresponding terminal velocity. Equation (1) can be written in advection form as
Here, the first term in the right-hand side of Eq. (2) can be combined with the term in the left-hand side to form a Lagrangian total derivative in time; the vertical derivative of the terminal velocity in the second term in the right-hand side of Eq. (2) can be written as
which is a total derivative or Lagrangian form of the flux-form advection of Eq. (1). Because ρqΔ represents the mass of the precipitating hydrometeor, the Lagrangian form, Eq. (4), tells us that mass conservation can be obtained locally as
where subscripts A and D represent arrival and departure locations, respectively, following the terminal velocity in a given time step, from time n to time n + 1, and Δ can be replaced by Δz as the vertical cell thickness.
The result of mass conservation in Eq. (5) can be derived from Eq. (1) by the traditional flux-form computation, as demonstrated in Fig. 1a. Note from now on, the variable at the interface of any cell is given by symbol hat over the variable. Assuming ẑA+ and ẑA− are two edges of a model grid cell as in traditional upstream schemes, the local change of the sedimentation at this cell may be calculated using Eq. (1), which can be written as
The fluxes through the ẑA+ and ẑA− model edges from time n to time n + 1 in term of Δt can be expressed by
as shown in Fig. 1a by gray-shaded and line-shaded areas, respectively. The overlaid area, which is gray shaded and line shaded, is the common area that belongs to both shaded areas. The difference of these two fluxes remains the unoverlaid areas as
Substituting the above results into Eq. (6), with some manipulation, we obtain
which is the same conclusion as in mass advection form of Eq. (5), instead of computing the flux of each cell edge as in traditional flux-form computation by Eq. (7). Thus, we can say that this proposed method in Eq. (5) [or Eq. (9)] is much simpler and easier to implement than traditional flux-form computation in Eq. (7). Furthermore, the results of Eqs. (9) and (5) should be always true, whether the departure or arrival cell edges are model grid cell edges or not. In other words, Eqs. (5) or (9) are valid with either an upstream or downstream scheme. This is the main thesis and discretization proposed by Juang (2007, 2008).
Since the falling speed depends on the water mass along its path, it is a good approximation to use initial terminal velocity for advection. Thus, we propose to use the terminal velocity at departure point, which is the model grid point, for the advection path. To discretize Eqs. (5) or (9) we utilized a schematic diagram to illustrate this method repeatedly in Fig. 1b with more detail. Since terminal velocity VT is a function of mass at the model layer, the terminal velocity ŵk at the interface of the model layer k used for advection can be obtained from VT at adjacent model layers, which can be expressed as
where α is the interpolation coefficient, and n is related to the number of neighboring points for interpolation. Based on different choices of n and α, we can have different orders of accuracy for interpolation. In this study, we use linear (n = 1) or cubic interpolation (n = 2) by terminal velocities at adjacent levels with two (αk−1 = αk = 0.5) or four points, respectively. Then any given departure cell depth at k has two terminal velocities ŵD,k and ŵD,k+1 at two interfaces ẑD,k and ẑD,k+1. Consequently, the arrival points at ẑA,k and ẑA,k+1 can be expressed as
where all variables in the right-hand side are given, and sedimentation on the left-hand side can be solved. To complete one-time step advection, the values at arrival cells should be interpolated to model cells. Since Eq. (12) states a local mass conservation, the interpolations from all arrival cells to all model cells have to be mass conservative.
Figure 1c is given to illustrate the mass conservation in one time step in a six-layer model as an example. Since each of the interfaces for a given model layer (cell) is the interface of the neighboring model layers, the mass conservation in each model layer in this Lagrangian advection indicates that total mass is conserved during advection (see the arrows from departure edges at time n to arrival edges at time n + 1) by the following expression as
as in Eq. (12). It is noted that the arrival location has to be interpolated back to the model grid point since the locations of arrival points are not model edges, as mentioned, which are even possible below the surface of the ground. As long as the interpolation is monotonic with a positive definiteness, mass conservation as shown in Fig. 1c produces a precipitation mass at the arrival point, which can be expressed above the surface of the ground, and below ground as
And it can be easily proved that
as total mass conservation. All arrival points in Eq. (14b) are negative, which means they lie within the layer below the ground as in Fig. 1c, and are represented by mR as total precipitation reaching the ground. The interpolation for all segments of the arrival cells above ground into model grid cells, in Eq. (14a), confirms mass conservation with monotonic features, as given in the appendix, using either the piecewise constant method (PCM), piecewise linear method (PLM), or piecewise parabolic method (PPM). Finally, the mixing ratio of hydrometeors at the model grid points at n + 1 time step is given by
with the assumption of no change in the density during advection, as in pressure and temperature fields.
b. Modification to avoid negative cell advection
Since the departure mass is positive as is the departure cell in the model layer, the arrival mass, in the left-hand side of Eq. (12) should always be positive if the arrival cell is positive. However, it is possible that the difference in falling velocity at a given cell edge can generate a zero or negative cell depth at arrival time; the arrival mass can then be negative in order to satisfy Eq. (12). To illustrate this point in detail, we start with the definition of arrival cell depth from Eq. (11), which can be expressed as
The condition preventing negative arrival cell depth in Eq. (17) results in the following condition:
which is called the deformation CFL condition (DeCFL), as mentioned by Xiao et al. (2003) and Smolarkiewicz and Pudykiewiez (1992). To satisfy the condition in Eq. (18), we altered the value of the sedimentation velocity at the bottom cell edge by scanning from the model top to the surface, which can be expressed as
where c should be less than 1 and is a tunable experiential parameter. The value of c is set as 0.05 in all tests, similar to those in Xiao et al. (2003) for modification of falling displacement.
c. Determine the mean terminal velocity
The terminal velocity at the departure cell can be a function of mass as
which was found to be a good approximation to determine the falling path at a given time step. However, it is possible to improve the accuracy of the falling path by using the mean terminal velocity during the sedimentation as follows.
First, Eq. (20) is used as the first guess of mean terminal velocity as
then, using Eq. (10) to obtain velocities at cell edges with the modification to avoid the DeCFL condition with Eq. (19), then the terminal velocity at arrival time could be obtained, using Eqs. (12) and (20), by
which is based on Eq. (11). Because terminal velocity varies with the density of falling hydrometeors, the mean terminal velocity from the departure to the arrival cells is given as
For further accuracy, we can iterate the above process by giving the mean terminal velocity in Eq. (24) to terminal velocity in Eq. (21), then following the procedure to get Eq. (22) and then back to Eq. (23). To avoid oscillation in convergence, we apply mean terminal velocity for the next iteration by the averaged value of current mean terminal velocity and previous mean terminal velocity. From the cases in sections 3 and 4, we found that the percentage of the maximal difference between the first-guess and the first mean terminal velocities is about 4%, between the first mean and the first-iteration terminal velocities it is about 0.3%, between the first-iteration and the second-iteration terminal velocities it is 0.04%, and between the second-iteration and the third-iteration terminal velocities it is 0.01%. Thus, the first guess is a good approximation, and first mean is accurate enough, even though the increase of iteration number guarantees enhanced accuracy.
3. Theoretical tests without a cloud scheme
We follow the experimental setup of Xiao et al. (2003) and use their two theoretical cases to examine the proposed SL scheme as described in the previous section. The initial distribution of rainwater is defined (following Kato 1995) in Xiao et al. (2003) as
where ρq is in grams per cubic meter, zc is 1.0 × 104 m, and D is in meters. The computational domain is 0 ≤ z ≤ 14 000 m. An equally spaced grid with Δz = 70 m is used for all experiments in this section. Two conditions of terminal velocity are examined, one is a constant terminal velocity as 5 m s−1 and the other is a function of ρq as
Figure 2 demonstrates that the new SL scheme developed in this study reasonably reproduces the results as shown in Fig. 2 of Xiao et al. (2003). It is also visible that the higher order of the interpolation method, PCM, PLM, and PPM, in that order, results in less diffusion, as compared to the analytical solution. Note that our PPM with a maximum of 0.9 g is very close to the result obtained from the cubic method of Xiao et al. with a maximum of 0.8 g, but with a higher accuracy, indicating that our scheme is less diffusive than theirs.
Figure 3 shows the results from the second case with terminal velocity as a function of rainwater amount without the correction of bottom edge terminal velocity due to DeCFL, as mentioned in section 2. The figure shows that there are erroneous features when the time step is larger than 30 s, which violates the DeCFL condition. Shown are the plots for the PPM results; the same behavior is observed with other interpolation methods.
Figure 4 shows the same configuration as in Fig. 3, but with the bottom edge terminal velocity correction given by Eq. (19) with c = 0.05. It is clear that the scheme is stable with reasonably shaped profiles, even up to 120 s as a time step, which is a much larger time step than in Xiao et al. (2003). Figure 5 shows the results of a comparison between runs with mean terminal velocity ranging from no iteration and iteration with a bottom edge terminal velocity correction. The left panels are without iteration but with mean terminal velocity, the right panels are from two iterations, the top panel is from PLM, and the bottom panel is from PPM. The comparison shown in Fig. 5 demonstrates that the impact of the revised terminal velocity by iteration is not significant; however, the refinement of terminal velocity could alleviate unstable noise irrespective of accuracy in interpolating the Lagrangian cell back to the model grid point, as can be seen in a comparison between Fig. 4 with a time step of 120 s and Fig. 5.
4. Tests of the proposed SL scheme in the WRF model
The WRF model used in this study is the Advanced Research WRF version 3.1 model (Skamarock et al. 2008), which was released in April 2009. The WSM3 scheme used in this study is a simple ice-phase microphysics scheme using the revised ice process treatment of Hong et al. (2004) after Hong et al. (1998). In this study, the aforementioned SL scheme is implemented in the WSM3 algorithm and evaluated against the current Eulerian advection method with subtime steps.
a. Evaluation of sedimentation for falling hydrometeors
A 1D idealized test bed was designed as in the previous section, but with different initial precipitation, as shown in Fig. 6. A typical amount of rainwater with a maximum of 10 g m−3 was assigned at 5 km. The integration time step is 120 s, which is the maximum value possible for the cloud microphysics computation (Hong et al. 1998). Further details in the computational procedures for the WSM3 scheme are available in Hong and Lim (2006). Note that the terminal velocity is not updated during the fall-term substep in the current version of the WSM3 scheme. In other words, the terminal velocity for the first substep is used without updating, given the falling particles. The resulting precipitation mass can be unrealistically large near the cloud bottom, as shown in Fig. 6a. This problem almost disappears when the velocity is updated every subtime step for minor loops (Fig. 6b). It is, however, noted that there is a tendency to increase the maximum value with time. Realizing the fact that values greater than 10 g m−3 are due to inaccuracy in treatment of forward advection by the first-order-accuracy Eulerian scheme, slower traveling below the cloud center results in an exaggerated maximum value, which leads to the unrealistically large value at 8 min.
Such noiselike precipitation near the surface is diminished when the SL scheme is used (Figs. 6c,d). The results from the PCM and PLM runs are similar, with some wiggling features of rain at the cloud center. The results from the PPM run are very similar to those from the PLM run (not shown). A close inspection reveals a less-smoothed contour above the cloud center when the PCM scheme is introduced (cf. contour at 8 min in Figs. 6c,d). This is due to the less accurate treatment of ŵk in Eq. (10) in the PCM than in the PLM, as illustrated in section 2. Meanwhile, the wiggling features of the maximum mass at 8 min were alleviated by revising the velocity at the arrival points by one-step iteration. The resulting precipitation amounts at the surface are very similar to each other, except for the Eulerian computation, without updating the terminal velocity during minor loops (Fig. 7). Despite the erroneous features in the current formula in the WRF model, the total surface precipitations at 60 min are similar because the hydrometeors aloft eventually fall.
b. Idealized thunderstorm experiment
An idealized 2D thunderstorm case follows the study of Hong and Lim (2006). The grid in the x direction was composed of 201 points with a 250-m grid spacing. The model was integrated for 60 min with a time step of 3 s, thus, the number of the sedimentation loop is 1 for the control (CTL) experiment. Note that the diffusion options are changed to conform to a 3D turbulent kinetic energy (TKE) option, whereas the diffusion coefficients are prescribed as in Hong and Lim (2006). Thus, the shape of the hydrometeors was smoothed in Hong and Lim (2006).
Based on the results in the idealized test bed, in this thunderstorm test bed we chose the mean terminal velocity for advection in the SL scheme, without revising the velocity at the arrival point. The distribution of precipitation particles at the mature stage of the storm is realistic, and differences between the CTL, PCM, and PLM experiments are not significant (Fig. 8). One interesting feature is the shape of the rain at the leading edge of the storm, particularly below 4 km. It seems to be a mammatus, produced by the differential sedimentation of raindrops in the presence of a near saturated environment. A relatively small amount of precipitation is distinct near the freezing level at about 5 km behind the storm center in the case of the PLM, as compared to the results from the other two methods. Another visible difference is the width of the precipitation near the surface, which is smaller in the PLM run than in the CTL and PCM runs. The time series of the surface precipitation demonstrate that with up to 35 min of integration time, the results from the CTL and PCM are similar, whereas the PLM and PPM schemes produce similar evolutions of rainfall intensity (Fig. 9).
A comparison of the above results from different sedimentation schemes demonstrates that the accuracy of the sedimentation treatment for falling hydrometeors can influence storm morphology and surface precipitation evolution by modulating the microphysics terms related to the falling precipitation (e.g., the accretion of water by precipitation and evaporation of falling particles).
The proposed SL scheme is simple and practical. Although some researchers have developed diagnostic precipitation for computational efficiency (e.g., Ghan and Easter 1992), our method is physically based and computationally efficient because it allows for full interaction among hydrometeors with large time steps.
The results we obtained using a theoretical test bed demonstrate that the performance of our method is comparable to those of conventional Eulerian approaches and other SL methods. The third-order PPMs outperform the PCMs and PLMs, based on the results from the example with constant terminal velocity, which were even more accurate than those of Xiao et al. (2003) with their highest-order scheme. The modification of the falling velocity at the bottom edge of any given cell from model top to model bottom avoids the DeCFL condition, which allows our scheme to incorporate large time steps, up to 120 s, with a negligible influence on the resulting precipitation profile. The modification that we implemented allows our scheme to have mean terminal velocity reduces noise and improves accuracy by about 4% for the advection path, but iterative conditions may not be necessary because the differences between iterations is only 0.01%–0.04%. Thus, iteration may not be necessary to obtain more accurate estimates of mean terminal velocity.
It is noted that all other microphysics schemes in the WRF, version 3.1 (see Skamarock et al. 2008) employ the same method as that in the WSM3 scheme. Given the minimum grid interval of 100 m (typically near the surface), the subtime step is activated even in the case of the model time step of 10 s. This is because the terminal velocity for raindrops of 10 g kg−1 can be as high as 10 m s−1. Therefore, less physical sedimentation behavior can appear in the WRF model, even in the high-resolution experimental setup, since the terminal velocity is not updated inside the minor loop, as shown in Fig. 6a. The resulting precipitation may not be significantly affected because hydrometeors aloft eventually fall. However, it is noted that the source–sink terms related to the falling precipitation would be harmed in the current implementation when other microphysics terms are considered.
Updating velocity in every minor loop is an immediate remedy, but this demands increased computing resources. Thus, the proposed SL scheme can be the ultimate resolution maximizing the accuracy of both sedimentation and computational efficiency. In terms of computational efficiency on a global model test bed following Shimpo et al. (2008), compared to 100% of CPU time in the case of the Eulerian method, the total computing times for the SL scheme with the PCM, PLM, and PPM were 84.2%, 91.4%, and 92.5%, respectively.
Considering both the accuracy and efficiency of the proposed SL scheme within the WRF test bed, the PLM is a good choice when terminal velocity is averaged at the departure and arrival points, without revising the arrival value by iteration. Based on the results of current tests, future work involving implementing this proposed SL scheme into other regional models with other cloud schemes can be easily performed without further–extra theoretical tests.
This study was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (2009-0080658), and the Korean Foundation for International Cooperation of Science and Technology (KICOS) through a grant provided by the Korean Ministry of Science and Technology (MOST). Two journal reviewers are thanked for their constructive and invaluable suggestions.
Monotonic Interpolation and Remapping
In section 2, the proposed scheme conserves mass as long as the integration of the given profile within any given model layer equals the mean value times the length of the model layer or vertical grid cell. While the limiter is applied to the profile within the range of local given values, it is monotonic with positive definite interpolation. For the PCM, it uses a mean value for all points within any given grid cell, so it is first order after integration, and automatically monotonic with mass conservation and positive definiteness.
For the PLM, a linear function passing through the mean value is used for all points within any given grid cell, the slope of the linear function is based on the values of the adjacent grid cells, thus it is a second order after integration, and it requires a specific treatment to be monotonic. The following is a summary of the procedures for the PLM to be monotonic.
First, the slope at all cell edges can be expressed as
where the hat indicates the value at a cell edge, as mentioned, and Q = ρq. Then the slope of any given cell is obtained by averaging the two nearby slopes at the cell edges. However, the slope at the cell is reset to be zero when the values at the two nearby cell edges have opposite signs. Thus, the slope at any given cell k can be written as
We can then determine the values at cell edges for any given cell k as
where superscript + indicates the value at the cell edge of ẑk+1, and superscript − indicates the value at the cell edge of ẑk. The last check for a positive definiteness is to make sure that there are nonnegative values at cell edges, with the condition as
Furthermore, the profile for the integration of this linear function at a given cell k can be written by values at cell edges as
where t = (z − ẑk)/(ẑk+1 − ẑk), which is between 0 and 1 from one edge to the other within a cell. Note that each cell has its own Qk+ and Qk−, thus Q may not be continuous at the cell edge, which is in a sense, piecewise.
For the PPM (Colella and Woodward 1984), a parabolic function is used as the profile for any given point within the given grid cell. The parabolic function is based on values of the nearby grid cells as well, and it was a third-order scheme after integration of the second-order profile. We followed “relaxed monotonicity for PPM” as proposed by Lin (2004) in this study with some corrections and modifications. First, we define the monotonic difference at any given cell k as
where sign, max, min, and abs hereafter are the same usage as FORTRAN intrinsic functions, and
Then, we determine the values at cell edges as
For a positive definiteness, the following check needs to be performed on the values at cell edges as
Consequently, the profile for integration of this parabolic function at given cell k can be written by the values at cell edges as
where, again, t = (z − ẑk)/(ẑk+1 − ẑk) from 0 to 1 within any given cell.
No matter which method is used, PCM, PLM, and PPM, the piecewise vertical profile they construct always satisfies the condition that the total integral of mass must be constant no matter what grid points are selected to be integrated, thus the total mass is conserved between the total integral of arrival cells and the total integral of model layers in two different groups with different vertical grid locations.
Corresponding author address: Dr. Song-You Hong, Department of Atmospheric Sciences, College of Science, Yonsei University, Seoul 120-749, South Korea. Email: email@example.com