An efficient semi-Lagrangian (SL) numerical scheme for the advection of positive-definite microphysical variables in cloud-resolving models is proposed. The scheme is a linear combination of SL schemes of the first and second order. The proposed scheme is compared with two high-order, monotonic, nonlinear Eulerian schemes, in two 2D tests. In the first test, advection of a positive scalar field with large gradients is performed within complicated idealized steering flows. In the second test, a hail storm is simulated using the 2D Hebrew University Cloud Model, describing the interactions between 15 size distributions of different microphysical quantities (each distribution is described by 43 mass bins). The proposed scheme is computationally efficient, has a low numerical diffusion and a high level of mass conservation accuracy, and preserves the sum of multiple advected variables as well as the linear relationships between the advected scalar variables. Although the proposed SL scheme is not exactly positive definite, the negative values and their number, which may appear only in the case of very strong gradients, are negligible. The scheme produces results similar to those of nonlinear Eulerian advection schemes while reducing computation time by approximately a factor of 10.
a. Application of advection schemes in atmospheric models
Advection schemes play an important role in atmospheric models, where physical entities such as velocity, temperature, density, humidity, etc., have to be advected to large distances, especially in large-scale and mesoscale models, where grid spacing may be of the order of several to dozens of kilometers. In such models, there are strict requirements for the advection schemes to be conservative and to avoid nonphysical oscillations that may cause unrealistic dynamical effects and negative values of positive-definite variables such as density and humidity. Such oscillations are typical of linear advection schemes of order 2 or higher and may be anticipated according to the Godunov theorem (Godunov 1959). The emphasis in Godunov theorem is that a high-order, monotonicity-preserving scheme must be nonlinear. Most numerical atmospheric models use complicated nonlinear, high-order advection schemes to transfer meteorological quantities (Wang et al. 2009; Ahmad et al. 2005; Yamaguchi et al. 2011).
The appearance of cloud-resolving models (CRMs) with grid spacing ranging from several hundred meters to 1–2 km has given rise to explicit cloud modeling using different microphysical schemes (Khain et al. 2015; Morrison et al. 2016; Khain and Pinsky 2018).
There are two main difficulties related to the use of high-order, nonlinear advection schemes in CRMs. The first is that they are highly time consuming. In CRMs with bin cloud microphysics, there are several hundred fields that have to be advected at each time step (Khain et al. 2015; Morrison et al. 2016; Khain and Pinsky 2018).
When applying nonlinear schemes, the advection coefficient matrix has to be calculated for each field separately every time step, thus requiring a considerable computation effort. Accordingly, in simulations of mesoscale convective systems using the Advanced Research WRF Dynamics Solver (Sobhani and Del Vento 2017), the advection of these hundreds of variables takes significant fraction of the overall computer time (more than 50%) and other computer resources. In linear schemes, on the other hand, the same advection matrix is used for all fields and has, therefore, to be calculated only once at each time step. Thus, the time devoted to the computation of the advection, is reduced by an order of magnitude.
The fulfillment of two conditions on advection of microphysical variables is highly desirable: 1) the advected sum of several variables (say masses in different bins) should be equal to the sum of the advected variables (say total liquid water content), and 2) the linear relationships of any two variables before advection should remain the same after advection. Failure to comply with these conditions can lead to certain errors, for example, errors in calculating supersaturation at the cloud edges (Grabowski and Morrison 2008) or distorting the shape of the distribution functions. As shown by Ovtchinnikov and Easter (2009) nonlinear advection schemes do not obey requirement 1 that, according to Ovtchinnikov and Easter (2009), may cause a distortion of the size distribution functions used in bin-microphysical models. Ovtchinnikov and Easter (2009) proposed several algorithms requiring a correction of the advected variables over the entire computational area, using the value of advected sum of these variables. This approach further increases the computational cost of advection schemes. Morrison et al. (2016) proposed the so-called Scaled Flux Vector Transport approach, which retains linear relationships between the advected scalars and reduces the computational cost by 10%–15% while producing solutions similar to the traditional approach; this has been considered to be a significant advantage.
Since microphysical variables are positive, nonlinear advection schemes apply various techniques (limiters) to avoid spurious oscillations and negative values (e.g., Zalesak 1979; Smolarkiewicz and Grabowski 1990; Thuburn 1997; Bott 1989; Skamarock 2006; Wang et al. 2009; Balsara and Shu 2000; Blossey and Durran 2008; Durran 2010).
Along with the Eulerian numerical schemes, the semi-Lagrangian (SL) advection schemes are widely used. A detailed overview of SL schemes is beyond the scopes of the study. Note only that numerous upstream SL schemes do not conserve mass (or its integral analog). Consequently, in the past decades, multiple, mass-conserving SL schemes have been developed (Leonard et al. 1995; Nakamura et al. 2001; Xiao and Yabe 2001; Cotter et al. 2007; Kaas 2008; Zerroukat et al. 2009; Juang and Hong 2010; Hansen et al. 2011). Despite the differences between the schemes, all mass-conserving SL schemes are based on the advection of air volumes whose boundaries are calculated by interpolation from the Eulerian grid and on remapping the advected quantities after the advection step on the regular grid.
Since SL schemes allow the utilization of significantly longer time steps (up to factor of 25; Robert 1982) than the Eulerian schemes, they have been used in many numerical weather predictions [see Bates and McDonald (1982) and the reference in this study; Williamson (2007)]. It should be remembered, however, that such long time steps become meaningless if large gradients or fast changes of the advected field and velocity field are present or if time steps are limited by characteristic times of physical processes.
b. Requirements for the advection of microphysical variables in CRM
It is a usual practice to use the same or similar high-order and nonlinear advection schemes for solving the motion equation for the advection of velocity components and scalar microphysical variables.
We would also like to stress the fundamental difference between the requirements for advection and for conservation of momentum and global thermodynamic variables (e.g., temperature and humidity) in large scale processes and in CRM for cloud species within clouds.
Whereas the changes in meteorological quantities at the synoptic scale are largely determined by advection and have characteristic time scales of several hours, the main changes in microphysical variables are determined by fast and intense microphysical processes with characteristic time scales of a few seconds, that is, much shorter than synoptic ones (Khain and Pinsky 2018). The generation and sink terms in microphysical processes are dominating in cloud dynamics and thermodynamics, and residential time of raindrops, graupel, hail in clouds is low (say, less than 10 min). It means that errors related to advection of these hydrometeors may not accumulate to significant values, remaining lower than those during treatment of microphysical processes. Although Morrison et al. (2018) showed that high-order, nonlinear advection schemes lead to some artificial broadening of an initially narrow droplet size distributions (DSDs) in the 1D cloud column with constant updraft velocity, more power sources of the DSD broadening such as mixing, collisions, in-cloud nucleation, size sorting supposedly contribute more strongly to the DSD evolution. The validity of this assumption seems to be supported by the results of the present study (see below). The dominating role of microphysics allows many authors to compare results of simulations of the same meteorological phenomena using different models and to attribute the differences to microphysical schemes but not to the differences in advection schemes (Khain and Pinsky 2018).
As for mass conservation, even a small error in airmass conservation will result in a forecast distortion and possible physical instability. At the same time, a loss of a few percent of, say, graupel mass content during a hail storm will not produce any effect on the storm parameters, and cannot be identified either by a general analysis of the storm structure or by comparison with observations.
Consequently, mass conservation in the advection of microphysical variables should be as perfect as possible, but the requirement for a precise 100% mass conservation can be relaxed. The errors related to the loss of mass in the advection of microphysical variables are acceptable, as long as they are less than the errors introduced by uncertainties of microphysical processes and the imperfect methods of solving the microphysical equations. For instance, errors in rates of condensation/evaporation or deposition/sublimation can lead to changes in the masses of hydrometeors (liquid water content, ice content) substantially exceeding those caused by imperfect mass conservation by advection. Accordingly, the requirements for the accuracy of the advection schemes for microphysical species are not as strict as for the advection of synoptic-scale variables.
At the same time, the description of cloud microphysical processes in CRMs imposes several highly desirable requirements for the advection schemes:
The advection schemes should be computationally efficient.
The advection schemes of microphysical values should be positive definite (or at least this condition should be valid with high precision).
The advection schemes should retain the following conditions of the advection operator H:
Here H:A(x, t) → A(x + uΔt, t + Δt), where x is source point and x + uΔt is the target point reached by an air volume moving with the velocity u to the end of time period Δt. Note that it is not simple to obey the latter conditions exactly, because the points x + uΔt do not necessarily coincide with regular grid points of the model meshes.
Note also that high-order advection schemes require the utilization of multiple neighboring grid points, which means that some grid points will be located outside the clouds (especially if cloud grid points are close to the cloud edge) and may lead to nonphysical influences on cloud microstructure. In this case, the advantage of high-order advection schemes over low-order schemes is no longer obvious.
Taking into account these features and requirements, we propose a hybrid SL scheme that is a blend of first- and second-order linear schemes. The first-order scheme is a monotonicity-preserving scheme, which suffers from large numerical diffusion. The second-order scheme is neither positive definite (as may be anticipated from the Godunov theorem) nor conservative but has little numerical diffusion. The combination of these schemes substantially decreases the deficiencies of each component. The scheme is tested by comparison with two monotonic advection schemes—first, in a given complicated, idealized 2D steering flows using an advection solver and, second, in simulation of a hail storm using the 2D Hebrew University Cloud Model (HUCM) with bin microphysics.
2. Semi-Lagrangian schemes
For the sake of convenient illustration, we consider a 2D equation of advection. The final formulas will be suitable for the 3D case as well. The advection equation is
where C is a positive quantity hereinafter referred to as concentration and u and υ are components of the velocity vector u. We assume that quantity C is conservative (dC/dt = 0) and so does not change during its advection. This assumption is always true for mass concentration, which is a conservative quantity and is valid for volume concentration as long as the flow is incompressible (see below).
Equation (1) is solved numerically on a finite-difference grid with grid points (xi, yi) = (iΔx, jΔy). In the semi-Lagrangian scheme, we assume that quantity at the arrival grid point (xi, yi) at time step t = (n + 1)Δt is advected from the departure point (x, y) = (xi − uΔt, yj − υΔt) (Fig. 1). Therefore, the concentration is determined by the interpolation scheme used to approximate. In this section we review the two-dimensional, bilinear and biquadratic semi-Lagrangian schemes (Bates and McDonald 1982; Colella 1990; LeVeque 2002), and present the “hybrid” weighted combination of the two. Since we propose to apply the scheme to the advection of microphysical variables using small time steps, we will assume that velocity does not change during one time step.
a. First-order semi-Lagrangian bilinear scheme
In the x direction,
In the y direction,
Here u+ and u− are
and similarly for υ+ and υ−. The scheme described by Eq. (2) is known as the CTU (Colella 1990; LeVeque 2002, section 20.2). It incorporates four nodes for the interpolation, and it is positive definite, as long as the Courant–Friedrichs–Lewy condition is satisfied, because all of the terms in the interpolation expressions are positive. Therefore, unlike high-order flux methods, it is not necessary to eliminate negative values. Because of its linearity, the coefficients αp and βq are determined only once at each time step and are used for the advection of all variables. Thus, while a bin microphysics scheme has several hundred microphysical variables to advect, the advection in the proposed scheme is reduced to a few arithmetic operations.
The CTU scheme is superior to the well-known upwind scheme. It was compared with the upwind scheme using a simple 2D configuration in which an initially rectangle-shaped fluctuation of constant magnitude C = 1 was advected by a given background wind. The upwind scheme leads to the distortion of the originally rectangular domain. At the same time, the CTU preserved the original shape much better (not shown). In addition, the upwind scheme required a shorter time step and suffered from higher numerical diffusion.
b. Second-order semi-Lagrangian biquadratic scheme
To decrease the numerical diffusion typical of the CTU, we tested second-order SL schemes, which can be obtained from second-order interpolation schemes written for the entire mesh area shown in Fig. 1. The BiQuadratic SL scheme (hereinafter BiQ) incorporates nine nodes (the center node and eight nodes surrounding it). The BiQ scheme is written formally in the same way as the CTU scheme in Eq. (2), but with different coefficients given as
and are defined in Eq. (3).
c. A hybrid semi-Lagrangian scheme
BiQ is neither conservative nor positive definite and tends, like other second-order linear schemes, to exhibit oscillations near sharp changes in concentration (section 11 of LeVeque 1990), in accordance with the Godunov theorem. At the same time, its numerical diffusion is low. The CTU, on the other hand, is positive definite but suffers from large numerical diffusion. To eliminate oscillations in BiQ in the zone of high gradients, we propose using hybrid schemes. In the present study we use a weighted combination of (1 − γ)CTU + γBiQ in the form , where [Eqs. (4) and (5)] and [Eq. (7)]. This combination of the two schemes reduces the numerical diffusion as compared with the CTU scheme (see below) while still preserving most of its positive definiteness. Traditionally, hybrid schemes are used to solve the diffusion advection problem and are a combination of the central and the upwind schemes (Versteeg and Malalasekera 2007). Below, we describe the implementation of the hybrid scheme Hyb-05, with γ = 0.5, and Hyb-08, with γ = 0.8. It is necessary to stress that the proposed SL schemes are linear and do not contain flux-correction factors, which makes them computationally efficient and retains the necessary relationships between the advected variables. The analysis of stability of hybrid schemes, their amplitude and phase errors as functions of the Courant number and parameter γ, is presented in the appendix. Note here that the scheme is stable at Courant numbers of less than 1.
Although the precise Hyb-05 and Hyb-08 schemes do not appear to have been discussed previously in the literature, the idea of blending schemes of different orders is ubiquitous in the development of advection schemes (e.g., TVD schemes or WENO schemes). For instance, Priestley (1993) proposed a weighted combination of linear and cubic semi-Lagrangian scheme. The proposed scheme was monotonic and obeyed mass conservation property. Often the weights of the higher- and lower-order schemes vary in space and across the different scalars, making these hybrid schemes nonlinear. The schemes, like that proposed by Priestley (1993), were developed for the utilization in large-scale models, like climatic models. The advantage of these schemes is the possibility to use large time steps. The main novelty of our scheme is its linearity and computational efficiency at small time steps, which makes it useful for CRM with bin microphysics, where its low computational cost, especially for multiple tracers, is a key feature.
d. Accounting for divergent flows
Throughout the above analysis, we assumed that the concentration remains constant along the fluid path:
where is the position vector in space. The concentration considered here is the volume number or mass concentration C = Cυ, that is, the number, or mass, of the advected particles per unit volume of the surrounding fluid. This is the concentration assumed in the scheme developed by Bott (1989) and in the TVD schemes developed by Sweby (1984), which are written in flux form. If the flow converges or diverges, however, Eq. (8) is not satisfied. It follows from the continuity equation that in the absence of sources and sinks the advection equation is
Therefore, after a time interval Δt the concentration at the computation node differs from the concentration at the departure point:
3. Comparison of different schemes
a. Flow configurations and initial concentrations
We compared between different schemes using idealized flow configurations. The compared schemes included CTU, BiQ, the hybrid schemes (Hyb), the Bott scheme (Bott 1989) presently used in HUCM, and the TVD-Superbee scheme (Sweby 1984; LeVeque 1990; Thuburn 1997; Kemm 2011; Zhang et al. 2015). The Bott scheme is a conservative, positive-definite advection scheme with low numerical diffusion, written in the flux form. The advective fluxes at the boundaries of grid cells are calculated using a polynomial of order l (here we used l = 2). These fluxes are normalized and then limited by upper and lower values. The scheme is stable for Courant numbers of less than 1. The TVD-Superbee scheme is also conservative and has very low numerical diffusion. When testing 2D flows, we used Godunov splitting by applying 1D TVD in each direction sequentially (section 19.5 of LeVeque 2002). Many years of successful use of the Bott scheme in the calculations of clouds and storms allows us to consider the Bott scheme as a benchmark.
Different background velocity and concentration fields have been used to test the numerical schemes (e.g., Tremback et al. 1987; Bott 1989; Blossey and Durran 2008; Lauritzen and Thuburn 2012). Our tests were similar to those used by Tremback et al. (1987) and Bott (1989). The goal of the comparison in this study was to compare the proposed semi-Lagrangian schemes with well-known nonlinear schemes, using complicated, sharp-gradient steering flows and highly nonuniform concentration fields, that is, under conditions that may exist in the atmosphere at sharp wind and concentration gradients (e.g., at cloud edges). All of the tests were performed using 2D configurations.
The initial concentration fields to be advected were prescribed in three forms: (i) as a uniform concentration in a rectangle that occupies part of the domain, (ii) as a uniform layer with zero concentration values outside the layer, and (iii) as periodic functions in the form of squares of sines and cosines on a rectangle (see Table 1). Forms i and ii were chosen to examine the ability of the schemes to preserve maximum values and maintain the initially high gradients of the concentration fields. Form iii was used to check the ability of the schemes to preserve the correct relationships between the advected fields (i.e., the linearity of the schemes with respect to initial conditions).
The domain size was L = 150 km in the horizontal direction and H = 25 km in the vertical direction (perpendicular to x). In all tests (with one exception; see below), the grid resolution was Δx = 300 m and Δz = 100 m. Two configurations of the flows were used. In the first configuration, a nondivergent circular flow was used:
where the vertical wavenumber mz = 0.5 and the horizontal wavenumbers are nx = 2. The velocity amplitude was A = 24 m s−1. The streamlines of the flow with nx = 2 are shown in Fig. 2. This configuration is similar to deformational flow (Tremback et al. 1987; Smolarkiewicz 1982). In the second configuration, an isolated circular flow was superimposed on a background flow (u, w) = (10; 0) m s−1 and the isolated circular flow was nondivergent:
The effect of the divergent flow was examined by eliminating the last term on the rhs of the equation for w. In Eq. (12) we use nx = 2 and mz = 0.5; x0 is the point of maximum u, that is, [nx(x0 + ϕ)]/L = 1/4 + N, with N = 1. The width of the cell was d ≈ L/20; the velocity amplitude was A = −20 m s−1. The phase ϕ controls the location of the cell. The streamlines of the flow corresponding to the cell are shown in Fig. 3. Details of different configurations used in the simulations are presented in Table 1.
We performed simulations both with no limiters (designated by “no bound”: nb) and with simple limiters. No special limiters were used for the Bott and TVD schemes, because the limiters are included in the formulation of the schemes. Two types of simple limiters were used: 1) bounding the solution from below by zero at each time step, (i.e., a low bound, designated lb) and 2) bounding the solution both from below by zero and from above by 1 or 2 with a low-up bound (designated lub). The limiters were introduced in to eliminate the oscillations that may possibly occur in zones of very high gradients in Hyb-05, Hyb-08, and BiQ. The following sections present the set of simulations with the different flow configurations.
We introduced two nondimensional measures, F and D, to characterize the properties of the numerical solutions; F is a measure of the oscillations of the numerical solution and is defined as
Measure F is zero for uniform functions. D represents the deviation from linearity and as such serves as a measure of the diffusivity of the scheme. It is defined as
Measure D is zero for linear functions. For a fully Lagrangian scheme, neither F nor D change in time (provided the flow is nondivergent); D decreases with increasing numerical diffusion, and it plays the same role as the entropy measure (Lauritzen and Thuburn 2012) but has the advantage of being nondimensional and normalized.
b. Circular flow with strong velocity gradients
The stream lines for this configuration are shown in Fig. 2. Figure 4 shows the fields of concentrations in the circular flow with mz = 0.5 and nx = 2, obtained using all of the schemes after 2 h. The maximum velocity was ±24 m s−1. This is a case of high-velocity gradients.
The CTU (Fig. 4a) demonstrates, on the one hand, the largest numerical diffusion, but on the other hand, the least oscillations. BiQ produces significant negative and positive values if no limiters are applied (Figs. 4b,c), but with the upper and lower limiters it looks similar to the Bott and TVD (Fig. 4d). Hyb-05 is less diffusive than the CTU and does not produce values beyond the maximum and minimum of the initial conditions, that is, larger than 1 and smaller than 0, even when no limiters are applied. Hyb-08 is less diffusive than Hyb-05, but produces negative values when no limiters are applied. The Bott (Fig. 4k) and TVD-Superbee (Fig. 4l) look much the same, with minor numerical diffusion.
The properties of the different solutions are shown in Table 2. In the table, m/m0 is the ratio of the integral of the concentration in the domain (for simplicity, referred to here as total mass) to its initial value; τ is the CPU time per one step; τb is the CPU time per one step for the Bott scheme, which serves as a reference; Δt is the time step; and Δtb is the time step for the Bott scheme. The Bott produces several values that are too high. The mass loss after 2 h is ~10% for CTU, ~5% for Hyb-05, and ~4% for Hyb-08, and all of the other schemes lose less. These mass losses are related to the strong gradients of the scalar field and high flow gradients. Mass losses of a few percent during 2 h can be considered as acceptable for the advection of cloud variables, especially if we take into account the uncertainties in description of many for microphysical processes.
The TVD-Superbee requires about ⅔ of Bott’s CPU time for one step, but the time step is one-half of the time step of all of the other schemes (for the sake of stability). At the same time, the CPU time required by the SL schemes is less by a factor of 7–9 when compared with the time required by the Bott.
Note that there exist linear flux-form schemes that are mass conservative. For instance, Tremback et al. (1987) proposed one linear advective scheme and a set of two linear flux-form schemes that are mass conservative and are of orders from 2 to 10. For the tests of the schemes, “deformational flow” is similar to the “circular” flow used in the present study. His flow field was composed of four cells in each direction, whereas ours consisted of four cells in one direction and one cell in the other.
Tremback et al. (1987) showed that the advective scheme was unstable. For the two-flux form with Courant = 0.7, the maximum amplitude was 18.04 [i.e., almost 2 times the initial amplitude (=10)]. The minimum was −0.92 (i.e., −10% of the initial amplitude). At similar time instance (2 h) in our Hyb-0.5 scheme computation the maximum amplitude was 1 (i.e., equal to the initial value) and the minimum was −0.02 (i.e., −2% of the amplitude). The mass loss was 6%. Note that mass loss in our SL scheme increases linearly with time. At times on the order of hydrometeor residential time in clouds the mass error in the SL scheme is below 1%. The higher-order schemes considered by Tremback indicated even larger amplitude errors than the second-order scheme. We, therefore, believe that for the purpose of advection of cloud microphysical fields the advantage of our scheme, with the low oscillations, overcomes the disadvantage of slight loss of mass.
c. Circular flow—isolated cell
The flow field for this configuration is shown in Fig. 3. It consists of a background flow of 10 m s−1 overlapped by a circular flow of 20 m s−1 defined in Eq. (12). The initial concentration was C = 1 for z ≤ 12.5 km and C = 0 for a larger z. Such sharp gradients of different variables can be observed between different layers both in the atmosphere and ocean. Figure 5 compares the fields of concentration obtained using the Bott, TVD-Superbee, Hyb-05, and Hyb-08 after 30 min.
The hill of C = 1 above the original layer (see Fig. 5) is formed by the rising air, and the trough of C = 0 is formed by the downdraft air with C = 0. The trough is swept along the background flow, but a small tongue of it is stretched to the left and upward by the negative horizontal flow near the bottom below the hill (Fig. 5). TVD and the Bott perform similarly, except for the narrow region of C > 1 in the Bott. Both Hyb-05 and Hyb-08 develop a very narrow band of concentration C > 1 in the interface between the main layer and the trough, with C = 1.05 for Hyb-05 and C = 1.2 for Hyb-08. Similarly, both hybrid schemes produce thin layers of negative concentrations around the hill and in the trough, with C = −0.070 for Hyb-05 and C = −0.169 for Hyb-08. These negative values are easily eliminated by using limiters, without creating any other changes in the results.
Table 3 shows the performance of the schemes in case of isolated circular flow. The mass conservation is satisfactory for all the schemes. The TVD scheme is the least oscillatory, with a slightly overestimated maximum. The Bott produces the maximum concentration, which is larger than the exact value by 10%. Among the two hybrid schemes, Hyb-05 is less oscillatory and more diffusive. The Bott and TVD are on the whole very similar. Again, the CPU time required by the semi-Lagrangian schemes is less than that of the Bott scheme by approximately a factor of 10.
To illustrate the effect of the flow divergence, we eliminated the last term on the rhs of Eq. (12) for w. As stated above, the Lagrangian schemes should be modified by a correction term, Eq. (10). The Bott and TVD do not need any correction because they are written in flux form. Figure 6 shows the concentrations after 30 min, obtained by the Bott, TVD-Superbee, Hyb-05, and Hy-08 (both modified with the correction term). One can see that the concentration fields of all the schemes are very similar (i.e., application of the correction term enables one to obtain correct results using the semi-Lagrangian schemes, even in a divergent flow). It is apparent that the correction in Eq. (11) leads to results similar to those produced by Bott.
d. Advection of sum of variables
As discussed in the introduction, advection by nonlinear schemes leads to the fact that the advected sum of variables is not equal to sum of the advected variables. In this section we are checking the accuracy of the condition on the advection operator H(A + B) = H(A) + H(B).
The first two examples that we consider are 1D problems with different shapes of initial concentration fields. In the third example we consider a 2D velocity field.
In example 1, the velocity is uniform in the x direction, (u; w) = (4; 0) m s−1, and the Courant number is equal to 0.89. We define two functions:
We define two concentration fields, C1 and C2, and examine their sum C3 with the following initial conditions:
In example 2, we use the same background wind field but the initial concentration consists of a set of unit step functions:
where xmax = 150 km. We decompose the initial concentration into two components, C1 and C2, and examine the sum C3:
Figure 7a shows the values C1 + C2 − C3 for the schemes tested. In the simulations Δx = 300 m and Δt= 67 s. In Fig. 7a, the first row shows results of example 1 [Eq. (15)] and the second row shows the results of example 2 [Eqs. (17)–(18)]. In all examples, the difference C1 + C2 − C3 for Hyb-05 is on the order of 10−15, that is, the order of the roundoff error. For the Bott scheme it reaches about 0.3% of maximum concentration for example 1 and about 15% (for the TVD scheme 10%) for example 2. So, the linear scheme Hyb-0.5 preserves the sum much better than the Bott scheme.
These errors of the Bott and TVD (not shown) schemes are located at the points of large second-order derivatives of the concentrations. Despite the lack of local conservation of sums by the Bott and TVD schemes, globally the sums are conserved as seen from Table 4 (example 2). For sake of convenience, Fig. 7b shows the values of C1, C2, and C3 of example 2 in the simulations using Hyb-0.5 and the Bott schemes at t = 2 h. The functional shape of these fields is very complicated, but the similarities between the two schemes are obvious. Bott scheme preserves the shape of C3 more accurately than Hyb-0.5, in which small overshoots and undershoots developed. Yes, Hyb 0.5 succeeds in preserving the sum C1 + C2.
In example 3, an example of 2D flow, we performed an experiment in which we simulated transport of three different initial concentrations on a rectangle, whose edges were located at L/6, 5L/6, H/6, and 5H/6 (L and H are the domain length and height; see Table 1) as follows:
Outside the rectangle, the concentration is zero. The flow was given in the form given by Eq. (11). Similarly to examples 1 and 2, we checked how well the linear relation C2 + C3 − C1 = 0 is preserved throughout the domain after advection using three schemes: Bott, TVD-Superbee, and Hyb-05.
TVD schemes are semilinear transformations; that is, for any constants a and b and a semilinear transformation T, the following relation holds (Lauritzen and Thuburn 2012):
If the relation C2 = 1 − C1 were valid over all of the computation domain, then we would have a = −1 and b = 1 (=C1) in Eq. (20) and the sum C1 + C2 = C3 = 1 would be preserved by the TVD scheme. However, Eq. (19) holds on the partial rectangle only, outside of which C1 = C2 = C3 = 0, so that the above relation does not hold globally. Furthermore, the initial concentration is either 0 or 1, but at a later time it may assume any intermediate value because of a small numerical diffusion that is present in all schemes.
Figure 8 shows the sum of fields C2 + C3 − C1 calculated in example 3 by the Bott, TVD, and Hyb-0.5 with and without application of low limiters, after 2 h of simulation. In the simulations Δx = 300 m, Δz = 100 m, u_max = w_max =16 m s−1, dt = 6.18 s, and the Courant number is equal to 0.98.
It is seen that the maximum values of the difference are produced by the Bott scheme. The maximum difference (i.e., error) can reach ±25% in some regions. The TVD scheme conserves the relationship between advected fields better, so the difference C2 + C3 − C1 is about ±10%. The linear scheme Hyb-05 produces no differences between C1 and C2 + C3 (only of order 10−15), in accordance with our expectation from a linear scheme.
Note that application of limiters makes the Hyb-05 nonlinear. To test the effect of the limiters a lower limiter lb = 0 was applied in an additional simulation. The results are shown in Fig. 8: the maximal value of C2 + C3 − C1 is 5 × 10−4, that is, negligibly small as compared with the values produced by the Bott and TVD schemes. Thus, the limiter lb = 0 applied to the Hyb-0.5 does not affect its linearity properties anyhow significantly. We attribute this effect to the fact of extreme rare appearance of negative values in this scheme.
e. Preservation of linear ratio between advected fields
In this test we use the same concentrations as given by Eqs. (15) and (16). Here we check to what extent the ratio between C1 and C2 is preserved by the two advection schemes Hyb-05 and Bott. To this end we define the function g(x, u, t) as
Function f2 depends on x − ut, because it should preserve the relation between C1 and C2 as they are advected by the uniform flow. Figure 9 presents the function g(x, u, t) at t = 800 s, t = 6000 s, and t = 7200 s. In the calculations u = 4 m s−1, dx = 300 m, and dt = 67 s; the Courant number is equal to 0.89.
One can see that g(x, u, t) is very close to unity for both schemes. Hyb-0.5 follows Bott very closely, except for the edges, but shows small oscillations about it. The deviation of Hyb-05 near the edges is due to the more diffusive nature of the scheme that introduces small values near the edges of the main structures. These small values are amplified when presented in terms of the ratio between C1 and C2, where C2 is small. For the Bott scheme g(x, u, t) is closest to unity at t = 7200 s, which is close to the period of f1 in Eq. (15) (7500 s).
4. Test of the SL schemes using HUCM
a. Brief model description and case study
HUCM is a 2D cloud model with spectral bin microphysics [see Ilotoviz et al. (2018), and the references therein]. The microphysical scheme describes the interactions between nine types of hydrometeors: aerosols, water drops, three types of ice crystals (plates, columnar, and dendrites), aggregates (snow), graupel, hail, and freezing drops. To describe the processes of melting and freezing, the liquid water mass within each bin of size distributions of freezing drops, graupel, snow, and hail is calculated. In addition, the rimed mass within each bin of snow size distribution is calculated. As a result, the microphysical scheme contains 15 size distributions described using 43 mass bins each (i.e., more than 600 positive-definite variables that must be advected at each time step).
Cloud microphysical processes such as diffusion growth, collisions, etc., are highly nonlinear, with multiple feedbacks, including positive feedbacks. For instance, formation of large hail leads to the acceleration of liquid water accretion by hail, accompanied by latent heat release during drop freezing. As a result, vertical velocities increase, which leads to further intensification of the accretion process and hail growth (Ilotoviz et al. 2016, 2018). Such nonlinearity and feedback loops make the results (cloud evolution) sensitive to both physical factors (such as aerosols, turbulence, instability of the atmosphere, and fall velocity of particles) and to the methods of solving microphysical equations. This high sensitivity of the scheme should be taken into account while comparing the results of applying different advection schemes.
Simulations were performed within the computational area of 300 km × 20 km, with a horizontal grid spacing of 300 m and a vertical resolution of 100 m. The dynamical time step was equal to 5 s. Microphysical processes of diffusion growth/evaporation and deposition/sublimation were calculated using microphysical substeps, which varied from 0.4 to 1 s. The velocity field, temperature, and humidity were calculated using the conservative Arakawa scheme (Arakawa 1966).
To test the applicability of the SL schemes, a severe hail storm was chosen that took place in Villingen-Schwenningen, in southwestern Germany, on 28 June 2006. The meteorological conditions of this storm were described by Khain et al. (2011) and are characterized by strong vertical wind shear, with northerly to easterly winds near the surface and a strong wind from the southwest at levels above 3000 m. The maximum reported size of hailstones was about 5 cm, the vertical velocities sometimes exceeded 40 m s−1, and the radar reflectivity reached 70 dBZ. The storm was successfully simulated using HUCM (Ilotoviz et al. 2016, 2018). The proposed SL schemes were applied to the advection of size distributions of hydrometeors and aerosols, and the results were compared with those obtained using the Bott and TVD.
b. Comparison of different advection schemes using HUCM
Since the results are highly sensitive to small fluctuations, it is difficult to expect snapshots of the same time instances to be similar when different advection schemes are used. Nevertheless, Figs. 10 and 11, which show the fields of liquid water content (LWC) at t = 30 min and at t = 60 min, obtained using different advection schemes, indicate a significant similarity.
At t = 30 min, the similarity between the fields is very high. The higher numerical diffusion of the CTU (the first-order Lagrangian scheme) is evident in a slightly smoother LWC field with a lower value of the LWC maximum. A significant similarity also remains at t = 60 min. However, there are larger differences between the snapshots. The CTU with maximum numerical diffusion now has significantly lower values of LWC (minimum rainwater content). While we consider the Bott to be the most precise scheme and TVD the next best, it is interesting that the difference between Hyb-05 and Hyb-08 turned out to be less than between the Bott and TVD.
Figure 12 characterizes the differences in the dynamics of the simulated storm caused by applying different advection schemes. Figure 12a shows the time dependence of the maximum values of the vertical velocity. One can see that after the spinup time, the averaged vertical velocity is about 30 m s−1 in all simulations. BiQ shows, however, maximum amplitude fluctuations, while the CTU indicates minimum vertical velocity, likely due to the highest numerical diffusion. Hyb-05 produces a Wmax closest to that in the Bott scheme. Hyb-05 even shows a good agreement with the Bott in the phase of fluctuations of the maximum vertical velocities. Figure 12b shows the vertical profiles of horizontally averaged velocity maximum. Again, the CTU shows the minimum value of the peak at z = 8 km (26 m s−1), whereas BiQ shows the highest maximum updrafts above z = 10 km. Hyb-05 provides vertical profiles of updrafts that are in the best agreement with the Bott.
Figures 13 and 14 characterize the warm (no ice) microphysical structure of the storm. Figure 13a shows the time dependence of the domain-integrated LWC. Again, Hyb-05 and Hyb-08 produce results closest to those of the Bott. Figure 13b shows the time dependence of LWC accumulated over the entire domain. One can see that BiQ and the CTU show the maximum and minimum values, respectively. The results of Hyb-05 and Hyb-08 actually coincide with those predicted by the Bott scheme. Figure 14a shows the vertical profiles of horizontally averaged cloud LWC, with the CTU producing the minimum values. Again, Hyb-05 and Hyb-08 produce results closest to those of the Bott. Figure 14b, showing the vertical profiles of horizontally averaged in-cloud droplet concentration, indicates that BiQ produces values larger than the other schemes. The CTU produces droplet concentration that is closest to that of the Bott. Hyb-05 also shows a good agreement with the results obtained by the Bott.
Figure 15 characterizes the ice microphysical structure of the storm. Figure 15a shows the time dependence of the domain-accumulated ice water content (IWC), which includes all the ice species. BiQ and TVD show the maximum and minimum values, respectively. The integral IWC simulated by the CTU is also lower than that of the Bott. Again, Hyb-05 and Hyb-08 produce results that are closest to those of the Bott. Figure 15b, showing the vertical profile of cloud-averaged IWC, indicates that BiQ again overestimates the IWC as compared with benchmark. The results of Hyb-05 and Hyb-08 are close to each other and are the closest to those of the Bott above z = 6 km.
To quantitatively characterize the rate of closeness of the simulations’ results to those obtained by the benchmark Bott scheme, we used the normalized square difference (NSD), determined as
Figure 16 shows NSD histograms calculated according to Eq. (18) for results shown in Figs. 13a, 13b, 14a, 15a, and 15b, with respect to the Bott results. One can see that in all of the cases (except Fig. 13b), the Hyb-05 results are the closest to the Bott, which agrees with the evaluation of this scheme (section 3). It is of interest that the difference between the results obtained by the Hyb schemes and by the Bott is less than the difference between the Bott and TVD-Superbee. This demonstrates that 1) even the schemes that are considered highly accurate can lead to different results being applied for solving such nonlinear phenomena as clouds and 2) the Hyb schemes, especially Hyb-05, can replace the current Bott scheme with tolerable loss of accuracy.
We conclude the comparison of the schemes by presenting snapshots of the radar reflectivity fields calculated using a polarimetric operator (Ryzhkov et al. 2011) that takes into account the existence of the liquid water fraction within hail and ice particles during their wet growth, as well as the rimed fraction in snowflakes (Ilotoviz et al. 2018). Figure 17 shows that at first sight the snapshots produced by different schemes look quite similar. However, a more detailed analysis shows that the CTU produces minimum radar reflectivity, while BiQ produces maximum reflectivity, which reaches the highest levels. It is remarkable that the snapshots are so similar at 1.5 h of the simulations performed with 5-s time increments. The results indicate that 1) even the schemes that are considered to be highly accurate can lead to different results being applied for solving such nonlinear phenomena as clouds (the Bott and TVD-Superbee) and 2) the Hyb schemes, especially Hyb-05, can replace the current Bott scheme with tolerable loss of accuracy while requiring approximately a factor-of-7 less CPU time.
Since the versions of the SL scheme differ by the value of weighting factor, determining the level of numerical diffusion, the simulations with HUCM (as well as the idealized tests) can be considered as the evaluation of the effects of numerical viscosity of advection scheme on concentration fields (in the idelized tests) and on dynamics and microphysics of hail storm in the HUCM simulation. In the case γ = 0, corresponding to the largest viscosity, the maximum values of vertical velocity were by 2–5 m s−1 lower and the integral values of LWC and IWC were about 10%–20% lower than in cases of lower viscosity (higher γ): Hyb-0.5, Hyb-0.8, BiQ, and Bott. Statistics (NSD) presented in Fig. 16 allowed us to choose the advection scheme that produces the results closest to those of benchmark (the Bott scheme).
c. Comparison of simulations with similar computational cost
It was of interest to compare results of benchmark scheme (Bott) with those of Hyb-0.5 and Hyb-0.8 in cases when the Bott scheme is used at the mesh with larger grid spacing to get approximately similar computational time in all cases. Accordingly, the storm simulation was performed using HUCM with the Bott scheme, but the horizontal grid spacing of 600 m and to 1200 m (instead of 300 m in benchmark simulations). Figures 13 and 15 show results of HUCM simulations with grid spacing of 600 m. One can see substantial deviation of the results from all those performed with grid spacing of 300 m. Utilization of grid spacing of 1200 m resulted in larger deviations. These results show that it is impossible to get similar results for similar CPU time as in case of application SL schemes using benchmark scheme of cruder resolution. The advantage of Hyb-0.5 and Hyb-0.8 is obvious. This result shows that a moderate change in the model resolution affects the results much stronger than the utilization of different versions of the SL scheme instead of the benchmark nonlinear complicated advection schemes.
A linear SL scheme, which is a combination of SL of the first order (CTU) and second order (BiQ), was constructed (Hyb-05 and Hyb-08). The multiple detailed tests of advection of a highly inhomogeneous concentration field within highly sheared flows show that Hyb-05 has a diffusivity lower than CTU and does not demonstrate artificial oscillations as BiQ does. Analysis the Hybrid SL scheme shows that it is stable at Courant number of ≤1. Analysis of the values of amplification factors and phase velocity errors shows that both amplitude and phase errors are not large and one can expect that the corresponding errors can be acceptable for description of microphysical processes.
In addition, Hyb-05 is nearly positive definite (without applying limiters) and has adequate (but not exact) conservative properties. Being linear, the hybrid SL obeys the property that the advected sum of variables is exactly equal to the sum of these advected variables. Utilization of the limiters from below (lb = 0) does not affect this property because of negligible number of negative values and their smallness in case of no limiters. Nonlinear schemes like that of Bott do not obey this property. The 1D tests show that both SL hybrid scheme and the Bott scheme that conserve well the linear ratios between advected variables. While in the case of good spatial resolution the accuracy of Hyb-05 was found to be close to that of nonlinear positive-definite advection Eulerian schemes (the Bott and TVD), its computational cost is lower by approximately a factor of 10.
Comparison with linear flux-form schemes, which are mass conservative (Tremback et al. 1987), shows that amplitude errors of these schemes substantially higher than those of the SL hybrid schemes discussed in the study. We, therefore, believe that for the advecting cloud microphysical fields the advantage of our scheme, with the low oscillations, overcomes the disadvantage of slight loss of mass.
In the present study we test the effects of different advective schemes on concentration fields in idealized flows and on dynamics and microphysics of a hail storm in simulations with HUCM. We found that the variation of parameter γ (weighting factor; see section 2c), which is a measure of diffusivity of the advection scheme, from zero (comparatively high diffusivity) to unity (low diffusivity) leads to changes of the integral values of LWC and IWC by about 10%–20% and maximum values of vertical velocity by 2–5 m s−1. These values are comparatively not large for such kinds of problems.
A detailed comparison of the advection schemes in simulations of a hailstorm using HUCM showed that the results produced by Hyb-05 and Hyb-08 are close to those produced by Bott’s scheme, which is taken as our benchmark. A good agreement was observed in time-averaged and space-averaged values, and even in particular snapshots. This agreement is quite remarkable if one takes into account the high nonlinearity of the microphysical processes, which make the results very sensitive to multiple factors. Note that in some cases the Hyb-05 results were closer to those obtained by the Bott scheme than the results obtained by the TVD-Superbee scheme. Remarkably, even the CTU, which is the simplest scheme of the first order, produces results close to those of the benchmark Bott scheme. It is noteworthy that the linear (in our case linear semi-Lagrangian) schemes possess very useful properties H(A + B) = H(A) + H(B) and H(A)/H(B) = A/B (where H is the operator of advection). These properties should minimize distortion of size distributions in the course of advection.
The conclusions support, therefore, our hypothesis that accurate results of cloud modeling can be obtained using a combination of simple linear SL schemes that require much less CPU time than high-order nonlinear schemes. The test simulations showed that the proposed SL schemes require less CPU time by a factor of 7–10 when compared with the nonlinear positive-definite advection Bott scheme. Taking into account that advection of positive-definite variables (especially in the bin microphysics schemes) consumes 70%–80% of the overall computer time, such a reduction of CPU time is of crucial importance. Therefore, we conclude that Hyb-05 complies with all the requirements of the schemes of advection of microphysical variables as formulated in section 1b above.
The supplemental storm simulations using HUCM show that more exact conservative schemes (like the Bott scheme) produce worse results than the SL hybrid schemes in simulations of similar computational cost, because the SL schemes can be run at higher resolution for the same computational cost as the Bott scheme at coarser resolution. This result demonstrates that a moderate change in the model resolution affects the results much stronger than the utilization of different versions of the SL scheme instead of the benchmark nonlinear complicated advection schemes.
In studies by Leonard et al. (1995), Lin and Rood (1996), and Skamarock (2006), different flux-form schemes were proposed that in the case of application of Strang time splitting may turn out to be also very computationally efficient. These schemes are not positive definite and indicate, sometimes, significant oscillations. Application of limiters may sometimes lead to serious conservation issues. Advantages and disadvantages of these schemes as compared with those of hybrid SL schemes require a special investigation.
There is another potential way to decrease computational time needed for advection. Advection of cloud microphysical quantities can be conducted only within the cloud area and its nearest surroundings. Since, clouds typically cover a small fraction of the computational domain, the area where advection is to be performed should be in principle much smaller than that of the entire computational area.
As the next step in the investigation, we intend to perform a test of the semi-Lagrangian schemes in the 3D framework.
This research was supported by the Israel Science Foundation (Grants 1393/14 and 2027/17) and the Office of Science (BER). Partial support for this work comes from Grants DE-SC008811, DE-SC0014295, and ASR DE-FOA-1638 from the U.S. Department of Energy Atmospheric System Research Program. We express our gratitude to Dr. M. Ovchinnikov for useful comments and remarks.
Analysis of Amplitude and Phase Errors
a. Linear stability analysis
Consider the one-dimensional case. Assume that at the grid point xk the velocity is positive, that is, uk > 0. Define the Courant number ε:
The BiQ contribution to the n + 1 time step is [see Eq. (7)]
Analysis of stability is performed using the von Neuman method; velocities are assumed to be constant in space and time. Concentration in point xk is represented by a Fourier series as
is the complex frequency, M = L/Δx, βm = πm/L, m = 1, …, M and L is the domain size. Because of the linearity there is no interaction between the different modes, and it is therefore enough to analyze the amplification of one Fourier mode. Substituting a single mode in Eq. (A7) to Eq. (A6) and defining the complex amplification factor
Using Eq. (A12) and considering its real and imaginary parts, one can get
The condition of the scheme stability is |Gm|2 ≤ 1. Let us denote δ = cosβmΔx; then −1 ≤ δ ≤ 1. Writing the module of the amplification factor in terms of γ, ε, and δ, one can get
We now show that ε ≤ 1 is a necessary and sufficient condition for the stability of the scheme for all γ. To find δ at which |Gm|2 reaches its maximum, we differentiate Eq. (A15) with respect to δ:
Taking into account that −1 ≤ δ ≤ 1, 0 ≤ 1 − δ ≤ 2, 1 − γ(1 − ε) = 1 − γ + γε ≥ γε, 1 + γε(1 − δ) ≥ 1, and −γεδ ≥ −γε, we get A0 ≥ −γε + (−γε) = 0. This inequality is valid at any ε. Assuming ε ≤ 1, we have from Eq. (A16):
Therefore |Gm|2 is maximal at δ = 1. Substituting δ = 1 into Eq. (A15), we have
Assuming now ε > 1, we get similarly to Eq. (A17) that
Therefore, |Gm|2 is minimal at δ = 1, so that |Gm|2 > 1; that is, the scheme is unstable at ε > 1.
Tables A1–A5 show the amplification factors |Gm| for the Fourier modes with wavelengths from to 2Δx to 40Δx for M = 40 grid points as functions of the Courant number and γ. Remember that γ = 0 corresponds to CTU and γ = 1 corresponds to the BiQ scheme. The amplification factor exhibits a weak monotonic increasing dependence on γ, indicating the decrease in numerical diffusion. The values of the amplification factor are very close to unity indicating small amplitude errors of the scheme.
b. Phase errors
The relative phase speed error is the ratio between the phase speed of each Fourier component of the numerical solution and the analytical phase velocity equal to u. It follows from Eq. (A9) that
Equating the imaginary part with that of Eq. (A12) we get the numerical dispersion relation
Since the numerical phase velocity is , the phase error can be evaluated as the ratio
Tables A6–A8 show the phase velocity ratio for the Fourier modes λ = 4Δx, 8Δx, and 16Δx and for M = 40 grid points as functions of the Courant number and γ. For λ = 2Δx, the numerical phase velocity vanishes for all Courant numbers and for all γ; for λ = 16Δx, normalized phase velocity is close to unity. For Courant number 1 the phase velocity ratio is 1 for all γ. For γ = 0.5, the phase velocity ratio is a monotonic increasing function of the Courant number and achieves the maximal value of 1 at Courant number 1. Analysis of the values of amplitude factors and phase velocity errors shows that both amplitude and phase errors are not large and that one can expect that the corresponding errors can be acceptable for the short time periods that are typical of microphysical processes.