Tracing the Rain Formation Pathways in Numerical Simulations of Deep Convection

Quantifying the microphysical process contributions to surface precipitation in numerical simulations can be challenging. This is due to the fact that many microphysical processes contribute to the formation and depletion of rain drops and there is almost always a spatial/temporal mismatch between where/when rain is formed and where/when it strikes the surface. In this work, we develop a tracing method that tracks the sources and sinks of raindrop mass and number as they are advected by the Weather Research and Forecasting model. Applying the method to an idealized squall line confirms that convective precipitation is dominated by warm rain processes (autoconversion and accretion) while stratiform precipitation is dominated by the melting of rimed and unrimed ice crystals. Sensitivity experiments in which the prescribed cloud drop number concentration is increased confirm the conventional wisdom that weakened autoconversion increases the fraction of raindrops originating from cold rain processes. The method also reveals that when applied to deep convection the Khairoutdinov and Kogan autoconversion scheme produces an excessive number of raindrops which are subsequently clipped in P3 microphysics to keep the rain size distribution within prescribed limits. This problem can mostly be mitigated by increasing the assumed radius for raindrops created by autoconversion.

2 of 14 on other predicted properties. Processes can then be individually parameterized in an attempt to capture the real-world complexity of evolving cloud systems. The fidelity of simulated precipitating systems is often assessed against observations of precipitation at or near the surface (e.g., Varble et al., 2011). The surface is not only the vantage point of ground-based observations but is the location where rain rate matters most for practical purposes of weather and climate simulation. In a numerical simulation, the parameterization of individual processes determines the surface rain rate. As a result, understanding the model's behavior necessitates linking precipitation to processes. The time lag between precipitation formation and fallout and the highly dynamic nature of precipitating systems makes this seemingly simple linkage challenging.
Here we present an approach that we call "process tracing" which keeps track of the fraction of rain mass and number created and removed by each microphysical process even as those rain properties are advected to the surface. The method is similar to previous applications of Eulerian tracing (e.g., Eiras-Barca et al., 2017) but is unique in that it is applied to processes within the cloud microphysics scheme. In Section 2 we present the method and its application to Weather Research and Forecasting (WRF). Section 3 shows results from the application of the method to simulations of an idealized 3D squall line. Section 4 concludes the paper.

Rain Process Tracing
We begin with the continuity equation of a cloud microphysical prognostic variable: is the hydrometeor mixing ratio, ( ) is the 3D horizontal wind vector, is the air density, ( ) is the hydrometeor's sedimentation velocity, ∆ * is the sub-grid scale mixing operator, and is the index of some microphysical process (e.g., accretion).
For every microphysical process , we define a process tracer such that: A process tracer is therefore a prognostic variable that's only microphysical process rate is and is sedimented and advected with the same velocity operator as the prognostic variable it is tracing ( ) . As long as process tracers and the prognostic variable they trace are initialized at the same value everywhere (typically zero), the sum of the former should balance the latter: Ovtchinnikov and Easter (2009) have shown that advection of tracers by higher order schemes can lead to a divergence between the right-and left-hand sides of Equation 3 implying it may not be exactly satisfied in the WRF model. However, we have found the in our application, associated advection errors do not exceed ∼1% (more on this in Section 2.3).
Our application is reminiscent of work conducted by Kumjian et al. (2015) who partitioned the rain category into rain originating from warm processes, melting of ice, and shedding from hailstones. However, the method is distinctly different in its goal. While the Kumijan method attempts to improve the simulation of rain drop size distributions via incorporating multiple rain categories, the process tracing method attempts to track process contributions without impacting the rain solution.
For this study, Equation 2 is applied to the rain mass and number mixing ratios ( and respectively). The model of choice is the WRF model (Skamarock et al., 2008) run with the Predicted-Particle-Properties (P3) microphysics scheme (Morrison & Milbrandt, 2015). The procedure for creating a process tracer is as follows.
First the process tracer is defined like any tracer within the WRF variable registry. The newly defined variable is 3 of 14 then passed into the cloud microphysics scheme where the appropriate process rates are applied to it. The tracer's sedimentation uses the same sedimentation flux computed for the prognostic variable it is tracing. Since WRF routinely advects and diffuses the tracer array, no coding modifications are needed for the interaction of the tracer with the model dynamics and parameterized turbulence.
We define a total of 6 mass tracers and 7 number tracers. For mass, we track contributions from autoconversion, accretion, riming/freezing, evaporation, melting of rime mass, and melting of pristine ice. The last two exploit P3's definition of riming fraction to partition the melting tendency: is the ice mixing ratio and rime is the riming fraction. For number, we track contributions from autoconversion, self-collection of rain, melting, riming/freezing, and evaporation. In most bulk microphysics schemes, the properties of the hydrometeor size distribution are periodically checked to ensure the parameterized maximum and minimum hydrometeor diameters do not fall out of the prescribed limits. When they do, the number is either clipped or supplemented. In this work, we treat these corrections as process rates (since they effectively are) and dedicate additional process tracers to track positive and negative corrections bringing the total to 7 number tracers.

Idealized Squall Line With the WRF Model
The test case is a 3D idealized squall line with a 0.5 km horizontal resolution, 250 m vertical resolution, and initialized with the analytical profile from Weisman and Klemp (1982). The domain extends 1,200 km zonally and 200 km meridionally with zonal wind linearly increasing from 0 to 13 m/s in the first 5 km. Lateral boundary conditions are open in the zonal direction and periodic in the meridional direction. An elevated warm column (+3 K) is placed at the center of the domain to initialize convection. The column is elevated to 1.5 km above the surface, has zonal and vertical widths of 3 km, and extends throughout the whole meridional direction. Random perturbations to the temperature (within the column) are invoked to break symmetry. Simulations are run for 13 hr with the last 6 hr used for the analysis in this study. All temporal averaging utilizes that timeframe. The large domain size ensures the squall line reaches a quasi-steady state in which its mean propagating properties become relatively constant ( Figure S2 in Supporting Information S1 depicts the squall line's evolution through hourly snapshots). To facilitate analysis, we define the leading edge of the squall line as the first zonal point in which the anomaly of surface potential temperature drops below −1 K. Squall line leading edges are associated with a sudden drop in temperature due to their cold pools (Rotunno et al., 1988). The first zonal point for which the temperature drops 1 K below the domain mean is thus used as a proxy for the squall line's leading edge. Defining this point as the beginning of the squall line (rather than some specific point along the zonal direction) helps decrease meridional asymmetry and thus mitigate the impact of meridional averaging. We refer to this technique as frontal realignment.
All sensitivity simulations in this study revolve around some alteration of the autoconversion parameterization developed by Khairoutdinov and Kogan (2000) and used as the default option in P3: Where , , and are parameters retrieved from fitting autoconversion rates to observations of stratocumulus clouds, th is the radius a raindrop is assumed to have upon being converted from a cloud drop, is the density of water, and is the density of air. The default value of th is 25 μm . Note that even though the parameters in Equation 6 were derived from stratocumulus observations, they are the default option in the P3 scheme which has 4 of 14 been developed to handle different cloud regimes. That said, Kogan (2013) derived an alternative set of parameter values for convective clouds which will be tested in this work. The cloud drop number concentration is prescribed in these idealized simulations. Because is a crude proxy for cloud concentration nuclei and therefore aerosol effects, we perform sensitivity studies with different values.

Diagnosing the Surface Precipitation Rate
While the tracers can quantify the process contributions to the rain mixing ratios at any height, this work will focus solely on the rate at which rain mass and number hit the surface. Since tracer mixing ratios sediment at the same sedimentation velocity as the prognostic variable they trace, surface mass precipitation and number precipitation as a function of specific process contribution rates could be evaluated as: is the surface mass precipitation rate, is the surface number precipitation rate, is the air density, and the mixing ratios along with the sedimentation velocities are evaluated at the lowest model level. A summary of the process rate contributions to and is provided in Table 1. Figure S1 in Supporting Information S1 compares instantaneous values of and to the sum of their traced process contributions. It can be seen that the tracer sums are virtually identical to and , mitigating any concerns about advection/diffusion errors. While this is expected for -because the residual tracers absorb some of the advection/diffusion error-the near perfect correspondence between and its associated tracers verifies adequate conservation for the purposes of this study.
Microphysical parameter changes can impact precipitation through a dynamic feedback or through an impact on precipitation efficiency. The former can be quantified through diagnosing total cloud condensation. Process tracing allows additional insight by categorizing microphysical processes into a warm rain conversion efficiency , a pristine ice to rain conversion efficiency , and a rimed ice to rain conversion efficiency from , defined as: is the column integrated condensation rate, is the column integrated rate of vapor deposition unto ice, and overbars denote domain averages.
With the conversion efficiencies, could be written in terms of condensation, vapor depositional growth onto ice, various conversion efficiencies, and evaporation: Since all our simulations use saturation adjustment, is effectively set by the dynamics. Perturbations to some cloud microphysics parameters may induce dynamic feedbacks that can alter surface precipitation or one of its contributions. To help separate out whether a precipitation contribution is brought about by a dynamic feedback or an adjustment to the inner workings of the cloud microphysics scheme we interpret changes in through as changes due to feedbacks, changes through one of the terms as changes due to microphysics, and changes through or evap as changes due to some combination of microphysics and dynamic feedbacks. Figure 1 shows a precipitation snapshot of the squall line after it reaches quasi steady state (9 hr into the simulation) without (a) and with (b) frontal realignment. Note that units of kg.m −2 .s −1 are converted to mm∕hr for all results presented. Along the zonal direction, precipitation starts off as convective (>10 mm/hr) before tapering to stratiform (<10 mm/hr) across a more extensive area. The 10 mm/hr threshold is based on conventional classification of convective versus stratiform rain rates (e.g., Leary & Houze, 1979). By contrasting panels (a) and (b), we cam see that frontal realignment reduces the squall line's asymmetry. That said, a good deal of meridional variability persists.

Results From Control Experiment
To visualize the meridional variability (and contrast it with zonal variability), Figure 1c ranks grid cells from highest (top of plot) to lowest (bottom of plot) along each zonal cross section before temporaly averaging. Temporal averaging is done to ensure the structural trend is robust. The ranking results in higher precipitation rates are displayed at the top of the plot while lower precipitation rates at the bottom. The wider top is a result of the ranking process itself, where grid cells exhibiting high precipitation rates preferentially occupy the top section of the plot. Note the same color scheme as Figures 1a and 1b is used thus yellow grid cells are convective while grid cells along the color gradient are stratiform. We can see that the convective region is made almost exclusively of convective cells (rather than some combination of convective and stratiform/non-precipitating cells). Similarly, the stratiform region is dominated by stratiform cells. Furthermore, there is a sharp transition between zonal cross sections dominated by convective precipitation to those dominated by stratiform. This result suggests that meridionally averaged precipitation rates are-at least qualitativelyrepresentative of local conditions. In other words, a cell is far more likely to be convective between the leading edge and ∼40 km behind the front and stratiform elsewhere. That said, meridional averaging averages out values of high precipitation. For this purpose, Section 3.4 includes a precipitation spectrum analysis where tracer contributions can be quantified on a per precipitation rate basis. All other results rely on meridional averaging after frontal realignment.
General features of the squall line vertical structure are depicted in Figure 2. As established in previous studies, vertical ascent results from buoyant air being shoved upward by a well-defined cold pool (Rotunno et al., 1988). The cold pool is driven in large part by rain fallout. A rear inflow jet brings dry air from aloft which enhances rain evaporation thus further reducing the cold pool's buoyancy. Hydrometeor contours superimposed on a 2D vertical velocity contour plot illustrate the coupling between squall line dynamics and the various condensed particle types. The ascending branch of the squall line is cloud dominant due to a high rate of condensation during ascent-driven cooling of moist air. A large anvil structure comprised of ice forms aloft and is advected horizontally as the ascending motions diverge. Rain forms and falls behind the ascending branch where it perpetuates the cold pool. Note that due to meridional averaging, mean vertical velocities -as depicted in Figure 2-are significantly weaker than local values.
The squall line's clear physical layout allows us to qualitatively deduce that rain originates from warm phase processes at the leading edge and becomes more dependent on cold phase processes further away from the leading  Figure 3a, where we show the temporally and meridionally averaged process contributions to the surface mass and number precipitation rates. At the leading edge, warm rain processes dominate through accretion. The contribution to mass by autoconversion is negligible and thus not shown. As expected, a fraction of the rain mass produced is evaporated as raindrops approach the surface. As we move away from the squall line's leading edge (i.e., toward the stratiform region), the contribution by melting of rimed ice starts to pick up but it is accompanied by equally large negative contributions from rain freezing and riming unto ice. This is a sensible result as we would expect active cold processes within a deep convective region to simultaneously rime rain while producing it by melting frozen hydrometeors. As we move into the zonally extensive stratiform region, we can see that precipitation is dominated by melting in terms of formation and evaporation in terms of removal. Some key differences between the stratiform and transition regions are the uptick in the relative contribution by pristine ice and the absence of a negative contribution by rain riming/freezing in the former. These findings are consistent with our understanding of stratiform precipitation in which we expect that ice growth by vapor deposition play a critical role in the formation of precipitation.
The number rates in Figures 3b and 3c show that autoconversion dominates positive contributions at the leading edge while self-collection significantly modulates the very large contribution by the former. Evaporation is another extensive negative contribution compared to the net number rate. As with the mass rate, cold rain contributions begin to pick up well into the deep convective zone. An interesting behavior that could be seen here is the fact that in and around the peak of precipitation, cold rain processes have a net positive contribution in mass but a net negative contribution in number. This is another interesting verification of intuition: the cold sector of deep convective clouds acts to simultaneously provide rain mass through melting while removing rain number though riming. That said, our quantitative analysis reveals that these numerical simulations have a very minor dependence on cold rain processes when it comes to convective precipitation.
The positive and negative correction tracers are alarmingly large-especially the latter. We hypothesize that mass growth rates are unable to keep up with an excessive production of rain number through autoconversion leading to inadvertent clipping by the P3 correction subroutines. We found a simple fix: increase the threshold radius within Equation 7 from 25 to 40 μm . The motivation for doing so stems from our suspicion that excessive number clipping is due to a copious production of rain drop number. We tested out this fix and show the updated rates alongside the control simulation in Figure 3. The new rates indicate the larger threshold radius results in far less worrisome correction contributions while having a relatively small impact on the resultant mass and number rates. We therefore retain the 40 μm fix for the rest of our sensitivity simulations. Meridionally ranked 2D distributions (as in Figure 1c) are shown in Figures S3 and S4 in Supporting Information S1 and confirm that changes along that axis of variability are also minimal. Table 2 summarizes key quantities of interest for the control and sensitivity experiments (which will be discussed in upcoming sections). The first two rows contrast the control experiment with and without the threshold radius adjustment. We can see the impact of the adjustment on the final rain rate is small and is driven primarily by a small decline in cold rain contributions. Interestingly, the contribution to total rain by warm rain is unchanged indicating that despite reducing the number autoconversion rate, accretion remains more or less the same.

Sensitivity to Autoconversion Parameters
The control simulations established the dominance of warm rain processes and thus their fundamental role in setting the shape of the convective rain zonal distribution. Kogan (2013) showed that the autoconversion parameters of the Khairoutdinov and Kogan (2000) scheme are not adequate for convective clouds and proposed a new set of parameters that effectively increase the sensitivity of autoconversion to cloud mass and number.
A sensitivity experiment in which the default autoconversion parameters are swapped for the ones proposed by Kogan (2013) (henceforth referred to as K13) is contrasted to the control simulation in Figure 4. Impacts on the rain mass rates are minimal. Cold rain processes slightly decline in their contribution at the front end of the squall line potentially due to a reduction in cloud water content reaching cold sectors of the storm. Contributions to number exhibit more distinct differences where we can see that K13 has a higher contribution from autoconversion which is substantially buffered by self-collection and evaporation changes. The net result is a similar distribution of rain number rate. The stronger autoconversion rate from K13 exacerbates the clipping problem that we partially addressed with the threshold radius change earlier. Table 2 shows that on the whole precipitation is unchanged when switching to the K13 parameters. There is a very slight increase in warm precipitation at the expense of cold precipitation which is consistent with the fact that we have enhanced autoconversion in K13. This could also be interpreted as a slight increase in the warm rain conversion efficiency as alluded to by comparing the values of . Since K13 has a stronger dependence on the cloud drop number concentration, it stands to reason that it may be more sensitive to changes in that microphysical property. This will be examined as part of the next sub-section.

Drop Number Concentration Sensitivity Experiments
Sensitivity simulations in which the cloud drop number is increased from the default 200 to 800 cm −3 are conducted using both the control and K13 autoconversion parameters. We contrast precipitation rates for the high simulations against their control counterparts in Figures 5 and 6. The effect of increasing the cloud drop number concentration is qualitatively similar for both sets of autoconversion parameters. The net rain mass curve is shifted back slightly potentially due to the delayed onset of raindrop formation. The warm rain contribution decreases throughout with melting of rimed and pristine ice-brought about by enhanced cold processescompensating for the decline. The precipitation peak increases in the higher experiments and can be attributed to the uptick in melting of rimed ice. There is stronger riming and freezing of rain drops which is consistent with 10 of 14 the fact that the higher experiments exhibit more active cold rain processes. At the leading edge of the squall line, the amount of evaporated rain decreases which is consistent with the results of Lebo and Morrison (2014) in which it was shown that a higher leads to larger and less numerous rain drops which in turn decreases rain evaporation within the cold pool. The meridionally ranked distributions ( Figures S5 and S6 in Supporting Information S1) corroborate the zonal narrowing of the convective precipitation band and reveal a slower tapering of convective cells in the stratiform region.
In examining the number rates, we can see a dramatic impact on rain number originating from autoconversion. However, self-collection (and to a lesser degree evaporation) substantially buffers the impact on autoconversion leading to rain number rates that are very similar. Rain number from melting of ice also increases for the higher cloud drop number experiments. By suppressing autoconversion, the higher experiments diminish the problem of excessive rain number clipping.
It should be noted that the increase triggers a slight increase in total condensation and vapor deposition unto ice (Table 2)   can decrease rain evaporation within the cold pool which in turns helps the squall line attain a more optimal balance between the shear and cold pool vortices. A more optimal balance would lead to stronger updrafts and thus higher condensation rates. That said, the condensation increase is small (∼5%) and the decline in the warm rain conversion efficiency helps ensure that the change in precipitation is nil. We can therefore conclude that increasing enhances condensation within the squall line but decreases its ability to convert condensate to precipitation.
Finally, Table 2 confirms that the K13 parameters are indeed more sensitive to the cloud drop number change incurring larger changes in warm and cold precipitation. Increasing leads to a ∼10% drop in while enhancing by 4% and by 1%-3%. We can therefore conclude that suppressing the autoconversion rate does shift this particular idealized convective system toward a stronger reliance on cold processes in keeping with conventional wisdom (e.g., Beydoun & Hoose, 2019;Fan et al., 2016;Khain & Lynn, 2009;Lebo & Morrison, 2014).

The Rain Probability Density Function
Realistic precipitating weather systems do not conform to the approximate meridional symmetry (see Figure 1) that we exploited in our analysis thus far. Furthermore, meridional averaging, while preserving the broad structural context of the squall line, averages out important details about the intensity of local precipitation. A more common approach to evaluating the spectrum of precipitation is the rain probability density function (pdf). Binning grid cells according to precipitation intensity allows quantitatively analyzing the occurrence of precipitation without assumptions about the structure of the system being analyzed. In this section, we demonstrate the tracers' utility for contexts in which creating a physical composite of a precipitating system is not feasible. To create the precipitation pdf, we calculate average precipitation rates (total and tracer contribution) in logarithmically spaced precipitation intensities ranging from 1 to 200 mm/hr and divide that rate by the bin width. For each bin, the precipitation rate for the total precipitation is used. Accordingly, the sum of all tracer precipitation rates is equal to the total precipitation rate in each precipitation bin.
Mass and number precipitation pdfs along with their traced process contributions are shown in Figure 7. The plots on the right-hand side are for the same precipitation pdfs but zoom in on lower precipitation values (<20 mm/hr) that are difficult to discern in the left-hand side plots. We should emphasize that only the total precipitation curve is a precipitation pdf since the values of the process tracers are not used in determining which frequency bin they belong to. The pdf is constructed using the total precipitation rate after which process contributions are computed at the given total precipitation frequency.
Positive contributions to the mass budget for high precipitation rate is very clearly dominated by accretion with a somewhat minor contribution from melting of rimed ice. Melting of pristine ice is not discernible on the left-hand side mass plot but could be seen being relevant for low precipitation rates. Evaporation is the strongest contributor to mass removal for most precipitation bins with higher precipitation rates seeing equally high removal from riming and freezing. Low precipitation rates have a much stronger reliance on cold rain processes. Specifically, cold rain processes dominate their warm counterparts for precipitation bins lower than ∼7 mm/hr. Melting of pristine ice is only higher than melting of rimed ice in the lowest precipitation bin (1 mm/hr) pointing to a surprisingly important role for riming in stratiform precipitation. Interestingly, the highest and lowest precipitation bins have precipitation rates higher than the strongest positive contribution by a particular tracer. This implies that at both the high and low end of precipitation intensity, rain is sourced from more than one process.
The number pdfs tell a similar story: melting of ice explains the number source for low precipitation rates while autoconversion dominates high precipitation rates. As already established in previous analysis, self-collection mirrors autoconversion and substantially reduces the number of raindrops that would otherwise hit the surface. A startling feature of the number pdfs is the fact that while the autoconversion source ceases to increase with increasing precipitation, the resultant number rate continues to rise. The reason for this is the decline of evaporation as a sink of rain number at the highest precipitation rates. Autoconversion overtakes melting as the largest contributor to rain number at around the same precipitation bin accretion overtakes melting of ice mass. Removal of number for cold processes dips through the mid-precipitation bins for reasons that remain unclear.

Conclusions
We introduced process tracers into the WRF model and applied them to an idealized 3D squall line. The tracers confirm long-held intuition that stratiform precipitation is borne from the melting of ice crystals (Houze, 1997).

of 14
Convective precipitation on the other hand was shown to be dominated by accretion and have a minor contribution from cold rain processes. Beyond fundamental understanding of an idealized test case, the tracers were shown to have some utility in microphysics parametrization testing. By tracing process contributions along the entire lifecycle of prognostic rain mass and number, we found that the latter underwent excessive clipping due to copious generation by autoconversion. A simple fix was increasing the threshold radius raindrops are assumed to be formed at. We also revisited the problem of aerosol-cloud interaction in deep convective systems by quantifying how various sources and sinks of rain responded to an increase in cloud drop number.
While physical insights into the sensitivity of a squall line to certain parameter perturbations were revealed, it is important to recognize that the sign and magnitude of the responses we saw here may very well change for different environmental conditions (e.g., shear). However, the tool itself is generalizable to other systems and could therefore be used to explore the robustness (or lack thereof) of a particular response. Future work could thus exploit the new tool introduced here in a rich ensemble of realistic WRF simulations that span various microphysics options as well as various meteorological conditions. An interesting path forward would be to use the tracers to track down process culprits in the mismatch between simulated and observed radar reflectivities. Another potential avenue is to use the tracers as part of an effort to understand precipitation physics in global cloud resolving models. These types of models create immensely high memory output files which severely limit the amount and frequency of traditional diagnostic process rates. By saving information associated with precipitation, process tracers can be written out at the same spatial and temporal resolution as the rain rate without the unnecessary burden of having to reconstruct process pathways offline. (c) Number surface precipitation along with process contributions plotted against frequency of occurrence (d) Same as (c) but zoomed in on lowest precipitation bins. Warm processes dominate at high precipitation intensities while cold processes play a larger role at low precipitation intensities.