Dynamic simulation of rainfall kinetic energy flux in a cloud resolving model

We present the first simulation of rainfall kinetic energy flux in a cloud resolving model. This demonstrates the potential for conducting erosion prediction studies using a regional climate model. Soil erosion is already a major global problem, and there is growing pressure on our land to deliver food and stability. Rainfall kinetic energy flux is an important variable in erosion prediction but is generally parameterized from intensity due to measurement difficulties. We show that a cloud resolving model can be used to dynamically simulate the kinetic energy of rain from basic physics, using four microphysics schemes. Results are within the range of observations and also capture the observed variability in kinetic energy for a given intensity. Large drops are shown to contribute disproportionately to total kinetic energy compared with their number, suggesting that several existing relations between terminal velocity and size of raindrops are poorly suited for kinetic energy modeling.


Introduction
[2] Rainfall kinetic energy is an important variable for soil erosion prediction and is commonly used in empirical and physical erosion models. Since the earliest work attempting to quantify the relationship between precipitation and soil loss [Wischmeier and Smith, 1958], it has been accepted that raindrop impact is a crucial determinant of soil loss [Rosewell, 1986]. Raindrops impacting the surface deliver the kinetic energy needed to detach soil particles from the Earth. Kinetic energy has therefore been used in various indices of rainfall erosivity which describe the erosive potential of rainfall and are used to predict rainfall-driven soil loss.
[3] Difficulties in measuring rainfall kinetic energy, however, have led to the development of equations parameterizing kinetic energy from intensity in a wide range of locations worldwide, based upon measurements of raindrop size spectra. This approach has several drawbacks. First, there is evidence of significant variability in kinetic energy for a given intensity of rain, owing to variations in drop size spectra and/or deviations of raindrops from terminal velocity [Rosewell, 1986;Kinnell, 1981;Fox, 2004;Parsons and Gadian, 2000;Steiner and Smith, 2000]. Second, since energy-intensity equations are developed based on data at a given location, their applicability to other regions is questionable. Third, emphasis is generally placed on providing a good fit for these parameterizations at high rain intensities, as these generally constitute the most erosive events. This can result in quite unrealistic predictions of energy at low intensities, such as positive values at zero intensity in the Revised Universal Soil Loss Equation [Renard et al., 1997] or negative values at low intensities in the Universal Soil Loss Equation [Wischmeier and Smith, 1958]. However, there is evidence to suggest that low intensity rain can in fact contribute significantly to erosion [Kirkbride and Reeves, 1993;Marques et al., 2008;Baartman et al., 2012]; so the use of many energyintensity equations would therefore be inappropriate in these situations.
[4] Here we present a new approach to soil erosion research, using a cloud resolving model to dynamically simulate the kinetic energy of rain from basic physics. We show that the model is capable of producing energies within the range of observations from several studies reported in the literature. These methods may be easily applied in a full regional climate model (RCM) with a microphysics parameterization scheme.
[5] Aside from avoiding the shortfalls of assuming a fixed relationship with intensity that have been outlined above, this approach has other advantages. Erosion models require climate data at high temporal and spatial resolutions. These data are often unavailable but are easily acquired from climate model output. The potential of climate model data to be used as input for ecological models has been recognized by Zhang et al. [2010], who used statistically downscaled mean monthly precipitation data from present and future General Circulation Model simulations to force an erosion model that requires intensity rather than energy data. In addition, using an RCM to simulate rainfall kinetic energy offers the opportunity to perform a range of studies investigating the role of the climate in erosion and of factors that may influence this. The approach presented here is limited by computing power, as convection resolving resolution is needed, but advances in this area will allow for climate and climate change simulations of kinetic energy, which may provide useful information for policy makers, land management and agricultural practitioners.

Methods
[6] Kinetic energy can be calculated for rainfall given assumptions about the size and speed of raindrops. The raindrop size distribution (DSD) is typically approximated by either an exponential [Marshall and Palmer, 1948] or a gamma distribution [Willis, 1984], which can both be expressed in the general form: (1) where N(D) is the number of drops of diameter D millimeter per unit volume, N T is the total number density of drops, is the gamma function, and and are the shape and slope parameters of the distribution. With a shape parameter of zero, equation (1) reduces to an exponential distribution. The fall speed of drops, V(D), is assumed to equal their terminal velocity, defined by a balance between gravitational and drag forces.
[7] Adopting the common assumption that raindrops are spherical with volume cD 3 where c = /6, rainfall kinetic energy flux can then be calculated as where w is the density of water. Taking a relationship between velocity at ground level and drop size of the form where a, b, and f are empirically determined constants, the solution to equation (3) is [8] Rainfall kinetic energy flux at the surface was modeled according to equation (5) using a simple cloud resolving version of the Weather Research and Forecasting model (WRF version 3.4) [Skamarock et al., 2005], a widely used nonhydrostatic mesoscale numerical weather prediction model. The cloud resolving version used for this study includes the complete microphysics parameterizations of the full model but contains no radiation or land surface scheme (no surface heat and moisture fluxes or frictional effects) and neglects effects of the Earth's rotation (the Coriolis force is set to zero).
[9] Microphysics schemes are employed to parameterize the subgrid scale precipitation processes occurring within clouds. Bin-resolving schemes explicitly simulate the DSD but are computationally costly. Consequently, bulk parameterizations are normally employed, which approximate the functional form of the DSD and prognose moments of this distribution for each hydrometeor class: single moment schemes predict the third moment (mass mixing ratio), while double moment schemes also prognose the zeroth moment (number density).
[10] For this study, we performed simulations of an idealized supercell storm using four different bulk microphysics schemes: Thompson et al. [2008] (Thompson), Morrison et al. [2009] (Morrison), Lim and Hong [2010] (WDM6), and Yau [2005a, 2005b] (MY). These all assume a DSD of the form given in equation (1) but use different values for the shape parameter which may be fixed or user defined. The six simulations will be referred to hereafter as Thompson-0, Thompson-1, MY-0, MY-1 (with = 0 or 1 as indicated), Morrison (fixed = 0), and WDM6 (fixed = 1). All these schemes are double moment in rain but vary in their treatment of frozen precipitation. ƒ (equation (1)) is diagnosed from the predicted mass and number concentration. Equation (4) is a general form of the V(D) relation assumed in all these schemes. Only Thompson adopts a nonzero value for f; the remaining schemes assume a simple power law function. We will show later that this has important consequences for modeling the rainfall kinetic energy flux.
[11] The idealized supercell storm simulated was the default WRF case study. The storm was initialized with no topography in a domain of area 84 84 km 2 with a hori- zontal resolution of 2 km and 41 vertical levels between the surface and the model top height at 20 km. Initial soundings for potential temperature and water vapor mixing ratio were taken from Weisman and Klemp [1982], with the tropopause positioned at 12 km. Vertical wind profiles from Weisman and Rotunno [2000] were used, but shear was extended to 7.25 km instead of 6 km, with quarter-circle shear (in both x-and y-directions) up to 2.25 km, and shear only in the x-direction above this. The storm was initiated by a potential temperature perturbation of 3 K centered in the middle of the domain at 1.5 km above the surface. The perturbation had a radius of 10 km in the horizontal and 1.5 km in the vertical direction and tapered to zero at the storm edge [Morrison and Milbrandt, 2011].
[12] The relationship between modeled kinetic energy and intensity was approximated by a power law following the same methodology as Steiner and Smith [2000, hereinafter SS00]: The coefficients A and B were chosen for this case to minimize the root mean squared error (RMSE, given by equation (10) in SS00) by iterative variation of B.

Results
[13] Figure 1 presents modeled kinetic energy results against rainfall intensity for each of the six microphysics schemes studied. Data was averaged within 1 mm h -1 intensity bins and averaged, and was found to be insensitive to bin size. These are compared with a sample of energy-intensity parameterizations (E(I)) from the literature (39 kinetic energy-intensity equations from a total of 22 observational studies chosen by Salles et al. [2002]) .
[14] Simulated kinetic energy flux against rain intensity is shown in Figure 2, for the Thompson-1 model run (plots for other schemes were similar). The red lines represent the upper/lower range of the SS00 observations, taken from their Figure 3. The mean (median) RMSE for all model runs was 54% (56%), compared with a mean of 36% from observations in SS00 (only intensities greater than 0.1 mm h -1 were included to replicate the analysis by SS00). Variability differed among schemes, with higher RMSE values for the Morrison (83%), MY-1 (74%), and WDM6 (62%) schemes, and lower values for Thompson-0 (25%), , and MY-0 (50%).
[15] Figure 3 depicts the variation in mean-mass diameter of raindrops, D M , the number density of raindrops, N T , and the mass mixing ratio of rain water, Q, with intensity for each microphysics scheme, at the surface. The mean-mass diameter of the raindrops was chosen as a characteristic diameter, but an analysis of the median volume diameter revealed similar behavior.

Data Comparison
[16] It can be seen from Figure 1 that the model is capable of simulating rainfall kinetic energies within the range of observations. Modeled kinetic energies for all schemes lie mainly within the bounds of the literature values for a given intensity. The only exceptions to this are the Morrison and MY-0 schemes, which produced higher kinetic  energies at low and high intensities, respectively, but which nonetheless performed well overall. Correlation coefficients between bin-averaged modeled and literature energy values exceeded 0.97 for all schemes. Within a given scheme (Thompson, MY), using a gamma distribution ( = 1) instead of an exponential ( = 0) slightly lowered the kinetic energy in agreement with Fox [2004]. This is attributable to a narrowing of the drop spectra, resulting in fewer large drops.
[17] Parameterizing energy from intensity assumes that any variability around this fit is negligible. However, energy fluxes deduced from radar data by SS00 (see their Figure 3) demonstrate that there is considerable variability in energy for a given intensity, a claim also supported by observational studies [Rosewell, 1986;Kinnell, 1981], theory, and literature reviews [Fox, 2004;Parsons and Gadian, 2000]. Modeled energy fluxes (Figure 2,  compare favorably with those observations, showing a large variance in energy for lighter rainfall which becomes more constrained at high intensities. Modeled energy fluxes are greater than observed at low intensities; this is more pronounced in similar plots (not shown) for Morrison, MY-0, and MY-1 and less pronounced for WDM6. At high intensities, energy flux is greater than the observations in the MY-0 and MY-1 schemes, and less in the Thompson-1 and Thompson-0 schemes. These features are consistent with Figure 1, where the schemes display similar behavior relative to the mean literature energy-intensity curve. An interesting feature is the apparent bounding of observations between two limits above and below the model data. This is also present in the observations of SS00, and the authors suggest these limits are dependent on drop fall speeds and on upper/lower limits to raindrop size.
[18] Comparing RMSE statistics for the energy-intensity power law fit (equation (6)), average model variability is greater than observed by SS00 (mean RMSE of 54% compared with 36%) but differs among schemes. Morrison, MY-1, and WDM6 exhibited larger variances, while Thompson-0, Thompson-1, and MY-0 were closer to the observed value. This variability stems only from variations in drop spectra. Assumptions about drop terminal velocity, and any deviations from this, are additional sources of error (see section 4.3).

Analysis of Microphysics Schemes
[19] Closer examination reveals differing behaviors among the six microphysics options. As expected, rainwater mass increases with intensity ( Figure 3c). However, Figures 3a and 3b reveal two categories of behavior: schemes that distribute this rain mass among roughly the same number of drops of increasing size, or among an increasing number of raindrops of the same size. In WDM6, MY-0, and MY-1, as intensity increases, raindrops become larger, but not more numerous. By contrast, Morrison, Thompson-0, and Thompson-1 simply produce more raindrops of the same size. This contradicts observational evidence that drop size increases with intensity [e.g., Willis, 1984]. In these schemes, excessive size sorting of drops in the column is controlled imposing a maximum drop size (Thompson) or minimum (Morrison). At higher intensities, additional raindrops are formed by melting frozen hydrometeors.

Importance of Raindrop Velocity
[20] For any variable X, the flux at the surface is equal to the product of X and its rate of arrival at the surface (the fall speed of raindrops). Rainwater mass flux (surface precipitation rate) therefore depends linearly on fall speed. Kinetic energy, however, itself varies with the square of velocity, so its flux has a cubic dependence on speed (equation (3)). Consequently, an accurate representation of the fall speed of larger drops becomes crucially important when energy is considered.
[21] Several relations between terminal velocity and drop size, V(D), have been proposed in the literature. Those assumed in the microphysics schemes tested are compared with observations by Gunn and Kinzer [1948], whose results agree well with measurements in other studies, in Figure 4a. Evidently, some V(D) equations have emphasized a good fit for smaller drop sizes at the expense of larger drops. Only the Thompson scheme avoids this problem. In general, V(D) relationships which include an exponential term [Uplinger, 1981;Atlas et al., 1973] can account for the saturation of the velocity curve as the drag increases on larger drops. Power law expressions [Liu and Orville, [1969]; Atlas and Ulbrich, 1977] result in large overestimations of fall speed for heavier drops; van Dijk et al. [2002] made this same observation.
[22] The importance of larger drops is illustrated by the schematic in Figure 4b. This shows the mean drop size distribution at the surface and the resulting mass and kinetic energy flux distributions with size (results were similar for all schemes). The peaks of the curves occur at progressively greater sizes: larger drops constitute disproportionately to the total mass and kinetic energy flux relative to their size. Raindrops larger than 3 mm constitute only 0.4% of the total number of drops but contribute 32% to the mass and 52% to the kinetic energy flux on average. These numbers include the bias resulting from the use of straight power law V(D) relations in the MY, Morrison, and WDM6 schemes. Despite this, large drops still contribute significantly to total mass (27%) and kinetic energy flux (40%) in the Thompson runs, which use a more appropriate V(D) fit. It is also worth noting the significance of larger drops for mass flux and therefore for surface precipitation. The sphericity approximation also holds less well for large drops, which are more oblate.

Limitations
[23] The success of this method will depend on the ability of RCMs to produce realistic precipitation rates and rain DSDs. Differences between microphysics schemes (e.g., the number of moments, treatment of drop breakup, and ice phase microphysics) may influence these factors. Model deficiencies and development will therefore have important consequences for rainfall energy flux simulations. However, here we have shown that for a given precipitation intensity, results are within the range of observations. [24] We have shown that it is possible to dynamically simulate the kinetic energy flux of rainfall with a cloud resolving model, using a range of microphysics options. Results are mainly within the range of observations reported in the literature.

Conclusions
[25] The model broadly captures the observed variability reported in other studies, which results from variations in drop size distribution for a given intensity. Parameterizing kinetic energy from intensity ignores this variability.
[26] In addition, emphasis has been placed on obtaining good fits at high intensities for these parameterizations, which can result in unrealistic results for light rain. This is avoided by the proposed method.
[27] Our findings suggest that it is important to use a V(D) relation which fits the full range of drop sizes well, especially for kinetic energy flux, but also for surface rainfall. Several commonly used power law relations overestimate larger drop speeds.
[28] Each scheme revealed strengths and weaknesses in different areas. In all simulations, the energy flux-intensity relationship was reasonable, but the Morrison and MY-0 options produced higher energies at low and high intensities, respectively. WDM6 and Morrison exhibited too much variability in energy compared with the observations of SS00. Thompson-0, Thompson-1, and Morrison do not capture the observed increase in raindrop size for higher intensities but compensate for this with a rise in drop number concentration. Only the Thompson scheme uses a velocity-drop size relation that accounts for the saturation of speed as drop size increases.
[29] The method proposed in this study may be easily extended for use in a full regional climate model. This would capture the variability in rainfall kinetic energy that parameterizing from intensity ignores, as well as providing high temporal and spatial resolution data to force erosion models. However, direct observations would be needed for verification in real data case studies. This method requires convection resolving scales, which have recently been achieved in a 10 year simulation [Kendon et al., 2012]. As computing power improves, climate and climate change simulations of rainfall kinetic energy flux will become possible. This will allow projections of rainfall energy that may be useful for erosion prediction and management.