Numerical framework and performance of the new multiple-phase cloud microphysics scheme in RegCM4.5: precipitation, cloud microphysics, and cloud radiative effects

We implement and evaluate a new parameterization scheme for stratiform cloud microphysics and precipitation within regional climate model RegCM4. This new parameterization is based on a multiple-phase one-moment cloud microphysics scheme built upon the implicit numerical framework recently developed and implemented in the ECMWF operational forecasting model. The parameterization solves five prognostic equations for water vapour, cloud liquid water, rain, cloud ice, and snow mixing ratios. Compared to the pre-existing scheme, it allows a proper treatment of mixed-phase clouds and a more physically realistic representation of cloud microphysics and precipitation. Various fields from a 10-year long integration of RegCM4 run in tropical band mode with the new scheme are compared with their counterparts using the previous cloud scheme and are evaluated against satellite observations. In addition, an assessment using the Cloud Feedback Model Intercomparison Project (CFMIP) Observational Simulator Package (COSP) for a 1-year sub-period provides additional information for evaluating the cloud optical properties against satellite data. The new microphysics parameterization yields an improved simulation of cloud fields, and in particular it removes the overestimation of upper level cloud characteristics of the previous scheme, increasing the agreement with observations and leading to an amelioration of a long-standing problem in the RegCM system. The vertical cloud profile produced by the new scheme leads to a considerably improvement of the representation of the longwave and shortwave components of the cloud radiative forcing.

through the explicit prognostic simulation of one or more hydrometeors.
Simpler microphysics schemes treat the cloud water prognostically and diagnose the precipitating water (e.g. Rotstayn, 1997;Pal et al., 2000). Observational data show that between -23 • C and 0 • C 25 the occurrence of supercooled water is not negligible (Matveev, 1984). Often cloud schemes diagnose the fraction of cloud water that is in the ice phase according to the temperature (e.g. DelGenio et al., 1996). The diagnostic of cloud water into liquid and ice components assumes implicitly that processes within the cloud are fast compared to the model time step, implying that the cloud variables are always in equilibrium. Therefore, a diagnostic representation is unable to describe the temporal 30 variability and evolution of mixed-phase clouds and a prognostic treatment of ice and water is necessary to represent the respective microphysical processes, including their contrasting sedimentation rates, in mixed phased clouds. More complex microphysics schemes have been therefore introduced to treat separately the cold and warm cloud microphysics by solving prognostic equations for cloud liquid water and ice (e.g. Fowler et al., 1996;Lohmann and Roeckner, 1996). These issues are es-35 pecially relevant as climate models approach high resolutions at which cloud physics processes, including convection, need to be explicitly described without the use of parameterization schemes (e.g. Prein et al. 2015).
The Regional Climate Model RegCM version 4 (or RegCM4) of the International Centre for Theoretical Physics (ICTP) is a widely used system that has been applied to local and regional seasonal 40 forecasting and climate change problems for all regions of the globe (e.g. Sylla et al., 2010;Diro et al., 2012a, b;Nogherotto et al., 2013;Coppola et al., 2014;Fuentes-Franco et al., 2014). The model has a wide choice of physical parameterizations for processes such as deep convection, but, to date, uses a simple diagnostic stratiform cloud scheme with a single prognostic cloud variable (Pal et al., 2000). There is a need not only to improve the representation of the cloud processes in the 45 RegCM modelling system, but also to conduct a comprehensive evaluation of the simulated clouds in RegCM integrations, which have received limited attention relative to the surface climate of the model.
In this paper we first present a description of the revised numerics and microphysics of the new 5phase prognostic parameterization scheme for stratiform clouds. The scheme is then tested in a series 50 of experiments with the RegCM4 run using the tropical band configuration of Coppola et al. (2012), which allows an analysis of the scheme's performance in different climatic settings. The cloud variables are compared to the existing RegCM4 SUBEX scheme, and the new parameterization is also assessed using the recently available COSP simulator package, which allows for direct comparison with a range of cloud-relevant satellite products, using model variables in a forward radiative trans-55 fer calculation to avoid uncertainties in retrieval techniques (Bodas-Salcedo et al., 2011). The final section summarizes the findings and makes suggestions for future developments of the scheme.
2 Methodology 2.1 Regional climate model The new cloud microphysics parameterization is introduced into the International Centre for Theo-60 retical Physics (ICTP) Regional Climate Model RegCM version 4. RegCM4 is a three-dimensional compressible, hydrostatic, primitive equation atmospheric model based on the dynamics of the NCAR mesoscale model Version 5 (MM5; Grell et al. 1994) and described in Giorgi et al. (2012).
In the current version of RegCM4 the resolved scale cloud microphysics is treated by the Subgrid Explicit Moisture Scheme (SUBEX, Pal et al. 2000), which calculates fractional cloud cover as 65 a function of grid point average relative humidity and includes only one prognostic equation for cloud water. Rain is calculated diagnostically and it forms when the in-cloud liquid water exceeds a temperature-dependent threshold (autoconversion). Rain is assumed to fall instantaneously within the model's time step and to grow by accretion of cloud droplets. SUBEX does not treat cold cloud microphysics and the fraction of ice is diagnosed as a function of temperature in the RegCM4 ra-70 diation scheme from radiative transfer calculations . The diagnostic split of ice and liquid water assumes that below -30 • C clouds consist of ice and for temperatures above -10 • C clouds are liquid only. This representation is an augmentation of an earlier scheme (Giorgi et al. 1993 which was in turn a simplified version of the scheme of Hsie and Anthes 1984).

75
The new cloud microphysics scheme considers cloud ice as a separate prognostic variable and also solves prognostic equations for rain and snow, accounting for the major microphysical pathways between these categories (Fig. 1). The model includes four hydrometeors in total: cloud liquid water and ice, rain and snow. Each variable is expressed in terms of the grid-mean mixing ratio q x (kg kg −1 ) and the governing equations for the mass mixing ratios of water vapour q v , cloud water q c , 80 cloud ice q i , rain q r and snow q s take the form: where S i includes the microphysical source and sink terms for each hydrometeor, representing the conversion of water substance between microphysical categories (see Figure 1). The second term on the right hand side represents the source of the variable q x from the layer above due to grav-  Forbes et al., 2011). The scheme has the advantage of being conservative, numerically economical, stable at all timesteps, and employs a numerical solution framework that is trivially expandable to a larger numbers of microphysical variables, facilitating the future representation of hail and graupel categories, or various ice crystal size bins. However, Tompkins (2005b) highlights that the scheme is highly diffusive for fast falling 100 species. Following Tompkins (2005b), the equations are solved using the upstream approach, which utilises the forward difference quotient in time and the backward difference quotient in space. For the time step n, dropping the large-scale advection and diabatic contributions as these terms are handled elsewhere in the model outside the microphysics scheme, the discretized equations are: (2) 105 It is seen that the microphysical pathways have been divided between two terms A and B, according to the timescale of the process they describe. Processes that are considered to be fast relative to to the variable q x by the explicit processes. We note that there is not a definitive justification for how microphysical processes are allocated to each solution category. As sedimentation is in the 115 downwards direction and there is no transport within the cloud scheme in the upward direction, the equations can be simply integrated one layer at a time from the top to the bottom layer of the model, making the solution numerically efficient as in each layer the solution of a n × n matrix equation is required, where n is the species number. An n = 3 category system at model level k is given by: where the index k − 1 represents the layer lying above the solution layer. Unlike implicit terms, explicit terms can possibly reduce a cloud variable to zero or negative values. In order to avoid this, and therefore to ensure that all variables remain positive definite at the end of the time step, the initial vector A containing the explicit source and sink terms is generalised using an anti symmetric matrix A, whose elements A xy > 0 represent a source for the variable q x and a sink for q y :  All the terms in the diagonal, A xx , represent microphysical source that are considered "external" to 120 the scheme, such as the cloud water detrainment from the (mass-flux) shallow and deep convection schemes. For each time step, before calling the solvers, the sum of all sinks of each variable is scaled to avoid negative values, a method that avoids negative values while guaranteeing total water conservation. For each microphysical pathway the change of phase is associated with a release or absorption of latent heat. Regarding the enthalpy budget, rather than summing the microphysics 125 pathways (as in the schemes of Tiedtke, 1993;Swann, 1994, for example), which can easily give rise in coding errors and resulting non-conservation when modifying microphysical parameterizations in operational and/or evolving models, the source/sink is calculated using the explicit conservation of the liquid water temperature T L defined as: 130 Since dTL dt = 0, the rate of change of the temperature is given by the equation: where L(x) is the latent heat (of fusion or evaporation depending on the processes considered), D qx is the convective detrainment and the third term in the brackets is the sedimentation term. Unlike the ECMWF IFS, the RegCM4 cloud fraction is not prognostic, but rather uses a diagnostic approach which has the advantage of simplifying the implementation and numerical cost, but has a number of disadvantages. The fractional cloud cover C is calculated following the semiempirical cloudiness parameterization developed by Xu and Randall (1996), which uses the large-scale relative 140 humidity RH and average condensate (cloud water and cloud ice) mixing ratiosq l = q l + q i to give implicit information concerning the subgrid-scale total water distribution (see review in Tompkins, 2002) and the resulting cloud cover: In theory, such a scheme also incorporates the impact of sub-grid temperature variability on cloud 145 fraction, since temperature fluctuations are implicitly incorporated into the statistics of the cloud resolving model simulations to which the scheme is fitted, however temperature fluctuations are likely underestimated in the small 2D domains used in Xu and Randall (1996), although Tompkins (2005a) showed that temperature variability is in general far less important relative to total water variability above the boundary layer. One key disadvantage of using a diagnostic cloud fraction approach is that 150 the treatment of ice supersaturation in the clear part of the model grid box at temperatures below -38C, such as in the scheme of Tompkins (2007) is not permitted. This is because standard RH based schemes (Sundqvist et al., 1989;Xu and Randall, 1996, e.g.) diagnose overcast conditions when the gridbox is saturated. Modifying the diagnostic relation to introduce a higher threshold for nucleation at cold temperatures (Koop et al., 2000) would not be able to represent the hysteresis between pre 155 and post ice nucleation, in other words, a separate memory is required of where in the grid box nucleation has occurred. The Tompkins (2007) scheme was able to use the prognostic cloud fraction to accomplish this by assuming the nucleation and subsequent ice crystal diffusive growth timescales was fast compared compared to the model timestep, thus assuming precisely ice saturated conditions in the cloudy portion of the grid box.

b) Condensation and evaporation
The formation of stratiform clouds associated with large-scale lifting of moist air or with radiative cooling is treated as a function of the variation in time of the saturation mixing ration, following Tiedtke (1993). In fact if the saturation mixing ratio decreases, condensation occurs while as it 165 increases evaporation takes place. The variation in time of the saturation mixing ratio can be written as: This equation shows that the rate of change of the saturation mixing ratio is linked to diabatic cooling (dT /dt) diab and to the vertical motion with a grid mean vertical velocity ω, where (dq sat /dp) ma is 170 the variation of q sat along a moist adiabat.
Condensation occurs when: The condensation rate C 1 is proportional to the amount of cloud and is equal to: and all the increase of cloud is a source of cloud water unless the process occurs within cold clouds, in which case condensation is a source of ice as homogeneous freezing takes place.
The scheme treats two processes that induce evaporation: the large scale descent and the diabatic heating, giving rise to E 1 , and the turbulent mixing of cloud air with drier environmental air, producing E 2 , so that the total evaporation E is given by: As opposed to condensation, evaporation is proportional to the increase of the saturation mixing ratio and to the amount of cloud following: It is reasonable to assume that the cloud water content within clouds is homogeneously distributed 185 in the horizontal direction, therefore the evaporation is not changing the cloud cover until it reduces to zero.
As assumed by Tiedtke (1991), the evaporation by turbulent mixing is considered proportional to the subsaturation of the environment so that: 190 where k = 3 · 10 −6 s −1 is the diffusion coefficient per unit area.

c) Condensation from detrainment
As an input from the convection scheme the microphysics scheme receives the detrained mass flux D that is assumed to condensate into cloud water or into ice diagnostically using a coefficient α, 195 function of temperature. This process is applied for all types of convection, namely deep, shallow and mid-level and represents an important extension of the model's cumulus parameterization.
The source of water/ice cloud content is given by: Autoconversion is the mechanism by which rain or snow droplets form from the aggregation of cloud water or ice particles. This process plays a crucial role in the development of precipitation. For this reason we have implemented four different parameterizations of the process, all following the form: where P is the autoconversion rate, P 0 the autoconversion rate once the autoconversion has started, and T ≤ 1 is a function that describes the threshold behaviour of this process (Liu and Daum, 2004).
The four parameterizations of autoconversion in the scheme employ different threshold functions: 210 an "all-or-nothing" approach, described in Kessler (1969) and three exponential approaches using smooth threshold functions.
The first following Sundqvist (1978): 215 the second following Beheng (1994): and the third following Khairoutdinov and Kogan (2000): The autoconversion of cloud droplets distinguishes between maritime and continental clouds by assumed to freeze instantaneously. For temperature above this threshold supercooled water and ice 230 are allowed to coexist, they are assumed to be well mixed and are distributed uniformly through the cloud. At temperatures below this threshold the liquid water is assumed to freeze instantaneously and the process is a source of cloud ice. The ice crystal is then assumed to grow at the expense of the water droplets through the Wegener-Bergeron-Findeisen process following Rotstayn et al. (2000).
The melting of ice and snow is parameterized taking into account also the cooling due to the evapo-235 ration of liquid water during the melting process. Therefore, the wet-bulb temperature is used instead of the dry-bulb one. Melting occurs if the wet-bulb temperature is greater than 0 • C. The part of the box containing precipitation is allowed to cool to T melt =0 • C over a time scale τ . The wet-bulb temperature T w is parameterized through a numerical approximation suggested by Wilson and Ballard (1999). All rain freezes in a time step if the temperature is lower than 0 • C. This process represents 240 a sink for rain and a source for snow. Since freezing would lead to an increase of temperature due to the latent heat release the scheme ensures that the temperature does no exceed the 0 • C threshold. For a more detailed description of the parameterization of microphysical processes we refer the reader to the IFS Documentation, Cy40r1, Part IV: Physical Processes.  with focus on the two extreme seasons, December-January-February (DJF) and June-July-August (JJA).

Precipitation
An unambiguous assessment of the effect of the new scheme on precipitation performance is extremely difficult. On the one hand, the simulation of precipitation is sensitive to the use of different  ison with the TRMM data, this tends to yield an improved agreement with observations over land and a deterioration over oceans. As already mentioned, this conclusion likely depends on the model configuration, however it is clear from Figure 2 that the new microphysics produces a realistic sim-295 ulation of precipitation, particularly over land, throughout the tropics and sub-tropics. It should also be mentioned that the MIC scheme itself is sensitive to different parameters affecting the production of precipitation, and in particular the ice and snow fall speed and the choice of the autoconversion threshold (Nogherotto, 2015).

Cloud fractions 300
In this section we present an analysis of the cloud fractional cover. This is accomplished by applying the COSP postprocessing tool to the model output to produce cloud variables comparable to those available in the International Satellite Cloud Climatology Project (ISCCP, Rossow et al. 1996) ISCCP. As already mentioned, this post-processing was carried out only for the seasons December in the ACCESS model by Franklin et al. (2013). Total cloud fractions are calculated by the model using the approach of Xu and Randall (1996) and the random overlap assumption, which tends to maximize total cloud cover. The evaluation of total cloud cover is carried out using the ISCCP D1 data set (Rossow et al., 1996) Figure 3 shows the total cloud cover in the SUB and MIC simu-310 lations for the selected seasons, postprocessed with COSP's ISCCP simulator. These are compared with the corresponding observed ISCCP total cloud amounts for the same seasons ( Fig. 3e and Fig.   3f). The ISCCP's observed total cloud fraction averaged over the domain is 66.07% and 64.66% for DJF and JJA, respectively. These values are 68.44% in DJF and 65.35% in JJA for the SUB run, and 61.52% in DJF and 60.04 % in JJA for the MIC (Table 2). Therefore, the SUB scheme produces 315 generally larger cloud fractions than the MIC, and the observations lie within the two model configuration data. In general, both schemes capture the horizontal distribution of clouds over the band    (Table 3), and 335 slightly lower than the CALIPSO product. The largest differences between the two schemes occurs in the simulation of high and medium level clouds. Compared to the SUB scheme, the MIC produces much lower values of high clouds ( 25% vs. 64% for the domain average) and greater values of mid-level clouds ( 11% vs. 7%), in both cases considerably increasing the agreement with the CALIPSO data. A possible explanation could be related to the different approach in treating the con-340 vective detrainment: while in MIC the detrainment produced by the convection scheme is given as an input to the microphysics scheme and is therefore subjected to microphysical processes, in SUB the detrainment is a source of cloud liquid water and is not involved in the formation of rain until the following time step. Another possibility is that the SUB scheme does not include ice physics, which would be dominant at high altitudes. For example, ice crystals tend to aggregate faster than 345 cloud droplets and thus precipitate more efficiently to lower levels. Note that the difference of the results between the assessments with the ISCCP (Figure 3) and CALIPSO (Figure 4) data suggests that the SUB scheme tends to overestimate optically thin clouds not detected by ISCCP. In fact IS-CCP is able to detect clouds with optical depths greater than 0.15-0.25 (over ocean and land), while CALIPSO can measure optically thinner clouds with depths greater than 0.03 Rossow et al. (1996).

350
An even more accurate analysis of cloud vertical distribution can be carried out with the use of the Multi-angle Imaging SpectroRadiometer MISR  data. MISR uses nine cameras spanning much of the range of angles over which cloud reflectivity varies, thereby leading to a more accurate retrieval of albedo than the use of a single camera. The MISR retrievals can be processed to produce joint histograms of cloud top height (CTH) and optical depth, although Naud et al. (2002) 355 found that in the case of multi-layered clouds MISR often "sees" through the thin upper level clouds and mostly refers to low level clouds layers. To compare with the MISR retrievals, we postprocessed 13 Geosci. Model Dev. Discuss., doi:10.5194/gmd-2016-31, 2016 Manuscript under review for journal Geosci. Model Dev.   the RegCM4 data with the MISR simulator described in Marchand et al. (2010). Figure 5   In this section we assess the cloud influence on the model radiation budget via an analysis of the CRF (Ramanathan et al., 1989), defined for the shortwave (SW) and longwave (LW) spectra as: where F is the net downward (i.e. downward minus upward) shortwave (SW) or longwave (LW) CERES observations (Wielicki et al., 1996). Figure 6 shows   Table 4).     effect (leading to closer agreement with observations) on the partitioning of the CRF in its shortwave and longwave components, although the total cloud forcing is similar to that of the old scheme due 400 to cancellation of biases. This is mostly attributed to the reduction of high level clouds found in the previous section.
18 Geosci. Model Dev. Discuss., doi:10.5194/gmd-2016-31, 2016 Manuscript under review for journal Geosci. Model Dev. 2. Conversely, the new scheme had a strong effect on the simulation of cloudiness, and in particular it produced to a decrease in simulated upper level thin stratocumulus clouds, which increased agreement with observations and lead to an amelioration of a long-standing problem in the RegCM system (e.g. Giorgi et al. 1999). In general, the new scheme improved the 420 vertical cloud profile in the model.

3.
Despite having a small effect on the total CRF, the new scheme considerably improved its partitioning into longwave and shortwave components. This is mostly because of the reduction of the upper level cloud bias in the original scheme noted above.
The preliminary tests described here of the new microphysics scheme introduced in RegCM4 provide