Vertical Dependence of Horizontal Variation of Cloud Microphysics: Observations from the ACE-ENA field campaign and implications for warm rain simulation in climate models

1. Physics Department, University of Maryland Baltimore County (UMBC), Baltimore, 21250, USA 2. Joint Center for Earth Systems Technology, UMBC, Baltimore, 21250, USA 10 3. Department of Geography and Atmospheric Science, University of Kansas, Lawrence, 66045, USA 4. Department of Mathematical Sciences, University of Wisconsin — Milwaukee, Milwaukee, 53201, USA 5. Center for Aerosol Science and Engineering, Department of Energy, Environmental and 15 Chemical Engineering, Washington University in St. Louis, St. Louis, 63130, USA 6. Environmental and Climate Science Department, Brookhaven National Laboratory, Upton, 11973, USA 7. Joint Institute for Regional Earth System Science and Engineering, University of California Los Angeles, Los Angeles, 90095, USA 20 8. Jet Propulsion Laboratory, California Institute of Technology, Pasadena, 91011, USA 9. Department of Hydrology and Atmospheric Sciences, University of Arizona, Tucson, 85721, USA 10. Pacific Northwest National Laboratory, Richland, WA 99354, USA 25 To be submitted to the ACP special issue: Marine aerosols, trace gases, and clouds over the North Atlantic


Abstract:
In the current global climate models (GCM), the nonlinearity effect of subgrid cloud variations on the parameterization of warm rain process, e.g., the autoconversion rate, is often treated by multiplying the resolved-scale warm ran process rates by a so-called enhancement factor (EF). In 35 this study, we investigate the subgrid-scale horizontal variations and covariation of cloud water content ( ) and cloud droplet number concentration ( ) in marine boundary layer (MBL) clouds based on the in-situ measurements from a recent field campaign, and study the implications for the autoconversion rate EF in GCMs. Based on a few carefully selected cases from the field campaign, we found that in contrast to the enhancing effect of and variations that tends to make EF>1, 40 the strong positive correlation between and results in a suppressing effect that makes tends to make EF<1. This effect is especially strong at cloud top where the and correlation can be as high as 0.95. We also found that the physically complete EF that accounts for the covariation of and has a robust decreasing trend from cloud base to cloud top. Because the autoconversion process is most important at the cloud top, this vertical dependence of EF should 45 be taken into consideration in the GCM parametrization scheme.

Introduction 49
Marine boundary layer (MBL) clouds cover about 1/5 of Earth's surface and play an important 50 role the climate system (Wood, 2012).  Boucher et al., 2013) and understanding of aerosol-cloud interactions (Carslaw et al.,53 2013; Lohmann and Feichter, 2005). Unfortunately, it turns out to be an extremely challenging 54 task. Among others, an important reason is that many physical processes in MBL clouds occur at 55 the spatial scales much smaller than the typical resolution of GCMs, making the simulation of 56 these processes in GCMs highly challenging. 57 Of particular interest in this study is the warm rain process that play an important role in 58 regulating the lifetime, water budget, and therefore integrated radiative effects of MBL clouds. In 59 the bulk cloud microphysics schemes that are widely used in GCMs ( Morrison and Gettelman,60 2008), continuous cloud particle spectrum is often divided into two modes. Droplets smaller than 61 the "separation size" * are classified into the cloud mode, which is described by two moments of 62 droplet size distribution (DSD), the droplet number concentration (0th moment of DSD) and 63 droplet liquid water content (proportional to the 3rd moment). Droplets larger than * are 64 classified into a precipitation mode (drizzle or rain), with properties denoted by drop concentration 65 and water content ( and ). In a bulk microphysics scheme, the transfer of mass from the cloud 66 to rain modes as a result of the collision-coalescence process is separated into two terms, 67 autoconversion and accretion:( ) = ( ) + ( ) . Autoconversion is defined as the 68 rate of mass transfer from the cloud to rain mode due to the coalescence of two cloud droplets with 69 < * . Accretion is defined as the rate of mass transfer due to the coalescence of a rain drop with 70 > * with a cloud droplet. A number of autoconversion and accretion parameterizations have 71 been developed, formulated either through numerical fitting of droplet spectra obtained from bin 72 microphysics LES or parcel model (Khairoutdinov and Kogan, 2000), or through an analytical 73 simplification of the collection kernel to arrive at expressions that link autoconversion and 74 accretion with the bulk microphysical variables (Liu and Daum, 2004). For example, a widely used 75 scheme developed by Khairoutdinov and Kogan (2000) ("KK scheme" hereafter) relates the 76 autoconversion with and as follows: 77 where and have units of kg kg−1 and cm−3, respectively; the parameter = 1350, and the 78 two exponents = 2.47, = −1.79 are obtained through a nonlinear regression between the 79 variables and and the autoconversion rate derived from large-eddy simulation (LES) with 80 bin-microphysics spectra. 81 Having a highly accurate microphysical parameterizationspecifically, highly accurate 82 local microphysical process ratesis not sufficient for an accurate simulation of warm-rain 83 processes in GCMs. Clouds can have significant structures and variations at the spatial scale much 84 smaller than the typical grid size of GCMs (10 ~ 100 km) (Barker et al., 1996; e.g., Cahalan and 85 Joseph, 1989; Lebsock et al., 2013;Wood and Hartmann, 2006;Zhang et al., 2019). Therefore, 86 GCMs need to account for these subgrid-scale variations in order to correctly calculate grid-mean 87 autoconversion and accretion rates. Pincus and Klein (2000) nicely illustrate this dilemma. Given 88 subgrid-scale variability represented as a distribution ( ) of some variable x, for example the 89 in Eq. (1), a grid-mean process rate is calculated as 〈 ( )〉 = ∫ ( ) ( ) , where ( ) is the 90 formula for the local process rate. For nonlinear process rates such as autoconversion and accretion, 91 the grid-mean process rates calculated from the subgrid-scale variability does not equal the process 92 rate calculated from the grid-mean value of x, i.e., 〈 ( )〉 ≠ (〈 〉) . Therefore, calculating 93 autoconversion and accretion from grid-mean quantities introduces biases arising from subgrid-94 scale variability. To take this effect into account, a parameter is often introduced as part of the 95 parameterization such that 〈 ( )〉 = • (〈 〉). Following the convention of previous studies, is 96 referred to as the "enhancement factor" (EF) here. Given the autoconversion parameterization 97 scheme, the magnitude of EF is primarily determined by cloud horizontal variability within a GCM 98 grid. Unfortunately, because most GCMs do not resolve subgrid cloud variation, the value of EF 99 is often simply assumed to be a constant for the lack of better options. The EF for KK 100 autoconversion scheme due to subgrid variation is assumed to be 3.2 in the two-moment scheme 101 by Morrison and Gettelman (2008), which is employed in the widely used Community Atmosphere 102

Model (CAM). 103
A number of studies have been carried out to better understand the horizontal variations of 104 cloud microphysics in MBL cloud and the implications for warm rain simulations in GCMs. Most 105 of these studies have been focused on the subgrid variation of . Morrison and Gettelman (2008) 106 and several later studies (Boutle et al., 2014;Hill et al., 2015;Lebsock et al., 2013;Zhang et al., 107 2019) showed that the subgrid variability and thereby the EF are dependent on cloud regime 108 and cloud fraction ( ). They are generally smaller over the closed-cell stratocumulus regime with 109 higher and larger over the open-cell cumulus regime that often has a relatively small . The 110 subgrid variance of is also dependent on the horizontal scale ( ) of a GCM grid. Based on the 111 combination of in situ and satellite observations, Boutle et al. (2014) found that the subgrid 112 variance first increases quickly with when is below about 20 km, then increases slow and 113 seems to approach to a asymptotic value for larger . Similar spatial dependence is also reported 114 in Huang et al. (2014), Huang and Liu (2014), Xie and Zhang (2015), and (Grosvenor et al., 2018) ). Furthermore, MODIS cloud retrieval product is known to suffer from 139 several inherent uncertainties, such as the three-dimensional radiative effects(e.g., Zhang and how the EF should be modeled in the GCMs. Second, a good understanding of the vertical 159 dependence of and variation inside of MBL clouds will also help us understand the 160 limitations in the previous studies, such as Z19, that use the column-integrated products for the 161 study of EF. Finally, this investigation may also be useful for modeling other processes, such as 162 aerosol-cloud interactions, in the GCMs. 163 Therefore, our main objectives in this study are to: 1) better understand the horizontal 164 variations of and , as well as their covariation in MBL clouds, in particular their dependence 165 on the vertical height in cloud; 2) elucidate the implications for the EF of the autoconversion 166 parameterization in GCMs. The rest of the paper is organized as follows: we will describe the data 167 and observations used in this study in Section 2 and explain how we select the cases from the ACE-168 ENA campaign for our study in Section 3. We will present cases studies in Section 4 and 5. Finally, 169 the results and findings from this study will be summarized and discussed in Section 6. 170 171

Data and Observations 172
The data and observations used for this study are from two main sources: the in-situ 173 measurements from the ACE-ENA campaign and the ground-based observations from the ARM 174 ENA site. The ENA region is characterized by persistent subtropical MBL clouds that are 175 influenced by different seasonal meteorological conditions and a variety of aerosol sources (Wood 176 et al., 2015). A modeling study by Carslaw et al. (2013) found the ENA to be one of regions over 177 the globe with the largest uncertainty of aerosol indirect effect. As such, the ENA region attracted 178 substantial attention over the past few decades for aerosol-cloud interaction studies. we adopt a * = 20 µ as the threshold to separate cloud droplets from drizzle drops, i.e., drops 209 with < * are considered as cloud droplets. After the separation, the and are derived from 210 the FCDP droplet size distribution measurements. As an evaluation, we compared our FCDP-211 derived results with the direct measurements of from the multi-element water content system 212 (WCM-2000) also flown during the ACE-ENA and found an excellent agreement. We also 213 performed a couple of sensitivity tests in which we perturbed the value of * by 5 µm. The 214 perturbation shows little impact on the results shown in sections 4 and 5. The cloud droplet 215 spectrum from the FCDP is available at a frequency of 10 Hz. Since the typical horizontal speed 216 of the G-1 aircraft during the in-cloud leg is about 100 m s-1, the spatial sampling rate these 217 instruments is on the order of 10 m for the FCDP. 218

Ground observations from ARM ENA site 219
In addition to the in-situ measurements, ground measurements from the ARM ENA site 220 are also used to provide ancillary data for our studies. In particular, we will use the Active Remote 221 Sensing of Cloud Layers product (ARSCL; (Clothiaux et al., 2000;Kollias et al., 2005) which 222 blends radar observations from the Ka-band ARM zenith cloud radar (KAZR), micropulse lidar 223 (MPL), and the ceilometer to provide information on cloud boundaries and the mesoscale structure 224 of cloud and precipitation. The ARSCL product is used to specify the vertical location of the G1 225 aircraft and thereby the in-situ measurements with respect to the cloud boundaries, i.e., cloud base 226 and top (see example in Figure 1). In addition, the radar reflectivity observations from KAZR, 227 alone with in situ measurements, are used to select the precipitating cases for our study. Note that 228 the ARSCL product is from the vertically pointing instruments, which sometimes are not 229 collocated with the in-situ measurements from G1 aircraft. As explained later in the next section, 230 only those cases with a reasonable collocation are selected for our study. 231

ACE-ENA flight pattern 233
The section provides a brief overview of the G1 aircraft flight patterns during the ACE-234 ENA and explains the method for cases selections for our study using the July 18, 2017 RF as an 235 example. As shown in Table 2, a variety of MBL conditions were sampled during the two IOPs of 236 the ACE-ENA campaign, from mostly clear-sky to thin stratus and drizzling stratocumulus. In this 237 study, we are interested in the RFs that encountered the drizzling stratocumulus clouds, since our 238 objective is to understand the implications of subgrid cloud variation for the autoconversion 239 process. The basic flight patterns of G1 aircraft in the ACE-ENA included spirals to obtain vertical 240 profiles of aerosol and clouds, and legs at multiple altitudes, including below cloud, inside cloud, 241 at the cloud top, and in the free troposphere. As an example, Figure 1a shows the horizontal 242 location of the G1 aircraft during the July 18, 2017 RF which is the "golden case" for our study as 243 explained in the next section. The corresponding vertical track of the aircraft is shown in Figure  244 1b overlaid on the reflectivity curtain of ground based KAZR. In this RF, the G1 aircraft repeated 245 multiple times of horizontal level runs in a "V" shape at different vertical levels inside, above and 246 below the MBL (see Figure 1b). The lower tip of the "V" shape is located at the ENA site on 247 Graciosa island. The average wind in the upper MBL (i.e., 900 mb) is approximately Northwest. 248 So, the left side of the V-shape horizontal level runs is along the wind and the right side cross the 249 wind. Note that the horizontal velocity of the G1 aircraft is approximately 100 m s−1. Since the 250 duration of these selected "V" shape hlegs is between 580 s and 700 s, their total horizontal length 251 is roughly 60 km, with each side of the "V" shape ~30 km. These "V" shape horizontal level runs, insights for larger GCM grid sizes. In addition to the hlegs, we also identified the vertical 263 penetration legs in each flight, referred to as the "vlegs", from which we will obtain the vertical 264 structure of the MBL, along with the properties of cloud and aerosol. 265

Case selection 266
As illustrated in Figure 1 a and b for the July 18, 2017 RF, the criterions we used to select 267 the RF cases and the hlegs within the RF can be summarized as follows: 268 • The RF encounters precipitating MBL clouds according to both pilot report and radar 269 reflectivity observations from the ground-based KAZR. 270 • The RF samples multiple continuous in-cloud hlegs at different vertical levels with the 271 horizontal length of at least 10 km and cloud fraction larger than 10% (i.e., the fraction of 272 a hleg with >0.01g m−3 must exceed 10% of the total length of that hleg) 273 • Moreover, the selected hlegs must sample the same region repeatedly in terms of horizontal 274 track but different vertical levels in terms of vertical track. Take the July 18, 2017 case as 275 an example. The hleg 5, 6, 7, 8 follow the same "V" shape horizontal track (see Figure 1a) 276 but sample different vertical levels of the MBL clouds (see Figure 1b). Such hlegs provide 277 us the horizontal sampling needed to study the subgrid horizontal variations of the cloud 278 properties and, at the same time, the chance to study the vertical dependence of the 279 horizontal cloud variations. 280 • Finally, the RF needs to have at least one vleg and the cloud boundary derived from the 281 vleg is largely consistent with that derived from the ground-based measurements. This 282 requirement is to ensure that the vertical locations of the selected hlegs with respect to 283 cloud boundaries can be specified. For example, as shown in Figure 1b according to the 284 ground-based observations, the hlegs 5 and 10 of July 18, 2017 case are close to cloud base, 285 while hlegs 8 and 12 close to cloud top (see also Figure 4). 286 The above requirements together pose a strong constraint on the observation. Fortunately, thanks 287 to the careful planning of the RF which had already taken studies like ours into consideration, we 288 are able to select a total of four RF cases as summarized in Table 3. We will first focus on the 289 "golden case-July 18, 2017 RF and then investigate if the lessons learned from the July 18, 2017 290 RF also apply to the other three cases. Azores high to the south (see Figure 2b), which is a common pattern of large-scale circulation 295 during the summer season in this region (Wood et al., 2015). The Azores is at the southern tip of 296 the cold air sector of a frontal system where the fair-weather low-level stratocumulus clouds are 297 dominant (see satellite image in Figure 2a). The RF on this day started around 8:30 UTC and ended 298 around 12 UTC. As explained in the previous section, we selected 7 hlegs from this RF that 299 horizontally sampled the same region repeatedly in a similar "V" shaped track but vertically at 300 different levels. The radar reflectivity observation from the ground based KAZR during the same 301 period peaks around 10 dBZ indicating the presence of significant drizzle inside the MBL clouds. 302 Among the 7 selected hlegs, the hlegs 5, 6, 7 and 8 are 4 consecutive "V" shape tracks, 303 with hlegs 5 close to cloud base and hleg 8 close to cloud top. The hlegs 10,11 and 12 are another 304 set of consecutive "V" shape tracks with hlegs 10 and 12 close to cloud base and top, respectively 305 (see Figure 1). Using > 0.01 gm−3 as a threshold for cloud, the cloud fraction ( ) of all these 306 hlegs is close to unity (i.e., overcast), except for the two hlegs close to cloud base ( =46% for 307 hleg 5 and =51% for hleg 10). The and derived from the in situ FCDP measurements for 308 these selected hleg are plotted in Figure 3  To obtain a further understanding of the vertical variations of cloud microphysics, we 323 analyzed the cloud microphysics observations from the two green-shaded vlegs 1 and 3 in Figure  324 1b. The vertical profile of the mean and from these two vlegs are shown in Figure 4a and 325 where is either (i.e., = ) or (i.e., = ). 〈 〉 and σ X are the mean value and 335 standard deviation of , respectively. As such the smaller the value the larger the horizontal 336 variation of in comparison with the mean value. As shown in Figure 4c reveals that each of the two modes in these bimodal distributions approximately corresponds to 358 one side of the "V" shape track. As aforementioned, for all the selected 7 "V" shape hlegs, the left 359 side is along the wind and the right side across the wind (see Figure 1). To illustrate this difference, 360 the across-wind side of the hleg is shaded in yellow in Figure 3. It is intriguing to note that the 361 from the across-wind side of the hleg are systematically larger than those from the along-wind side,

Implications for the EF for the autoconversion rate parameterization 381
As explained in the introduction, in GCMs the autoconversion process is usually 382 parameterized as a highly nonlinear function of and , e.g., the KK scheme in Eq. (1). In such 383 parameterization, an EF is needed to account for the bias caused by the nonlinearity effect. If the subgrid variations of and , as well as their covariation, are known, then the EF 389 can be estimated based on its definition as follows 390 where 〈 〉 and 〈 〉 are the grid-mean value, ( , ) is the joint probability density function 391 (PDF) of and . Some previous studies approximate the ( , ) as a bivariate lognormal 392 distribution as follows: 393  autoconversion rate is proportional to 2.47 , the negatively skewed also leads to a negatively 454 skewed in Figure 7b. As a result, the leg-averaged diagnosed from the observation is slightly 455 smaller than that derived based on Eq. (6) by assuming a lognormal distribution. The negative 456 skewness also explains the large error in for hleg 10 seen on Figure 6b. As shown in Figure 7c  457 the observed is also negatively skewed, to a much larger extent in comparison with . Because 458 the autoconversion rate is proportional to −1.79 , the highly negatively skewed results in a 459 highly positively skewed in Figure 7d. As a result, the diagnosed from the observation is 460 much larger than that derived based on Eq. (7) by assuming a lognormal distribution. 461 The and reflect only the individual contributions of subgrid and variations to 462 the EF. The effect of the covariation of and , i.e., the is shown in Figure 6c. 463 Interestingly, the value of is smaller than unity for all the selected hlegs. As explained in Eq. 464 < 1 is a result of a positive correlation between and , as seen in Figure 4d. Finally, it is intriguing to note that the value of = • • in Figure 6d is 487 comparable to Figure 6a, which indicates that the enhancing effect of > 1 in Figure 6b is 488 partially canceled by the suppressing effect of < 1 in Figure 6c. As aforementioned, many 489 previous studies of the EF consider only the effect of but overlook the effect of and . 490 The error in the studies would be quite large if it were not for a fortunate error cancellation. 491

Other Selected Cases 492
In addition to the July 18, 2017 RF, we also found another 3 RFs that meet our criterions 493 as described in Section 3 for case selection. As summarized in Table 3 and shown in Figure 1c  As aforementioned, many previous studies of the EF for the autoconversion rate 534 parameterization consider only the effect of subgrid variation but ignore the effects of subgrid 535 variation, and its covariation with . To understand the potential error, we compared the 536 and both derived based on observations in Figure 10. Apparently, is significantly larger than 537 for most of the selected hlegs, which implies that the considering only subgrid variation 538 would likely lead to an overestimation of EF. This is an interesting result. Note that ≥ 1 by 539 definition and therefore > is possible only when the covariation of and has a 540 suppressing effect, instead of enhancing. Once again, this result demonstrates the importance of 541 understanding the covariation of and for understanding the EF for autoconversion rate 542 parameterization. 543 Having looked at the observation-based EFs, we now check if the EFs derived based on 544 assumed PDFs (e.g., lognormal or bi-variate lognormal distributions) agree with the observation-545 based results. As shown in Figure 11a, the based on Eq. (6) that assumes a lognormal 546 distribution for the subgrid variation of is in an excellent agreement with the observation-based 547 results. In contrast, the comparison is much worse for the in Figure 11b, which is not surprising 548 given the results from the July 18, 2017 case in Figure 6b. As one can see from Figure 5, the 549 marginal PDF of is often broad and sometimes even bimodal. The deviation of the observed 550 PDF from the lognormal distribution is probably the reason for the large difference of in 551 Figure 11b. As shown in Figure 11c, the derived based on Eq. (5) by assuming a bi-variate 552 lognormal function for the joint distribution of and tends to be larger than the observation-553 based results. The reason for this overestimation is because the joint PDF of and is often 554 bimodal as seen in Figure 5. In such case, the small correlation coefficient due to the bimodality 555 is misinterpreted as a rather broad bi-variate lognormal distribution which in turn leads to an 556 overestimated value. • In a few selected "V" shape hlegs, the and follow a bimodal joint distribution 571 which leads to a poor linear correlation between them. The two modes in the bimodal 572 distribution correspond to the along-wind and cross-wind sides of the "V" shape hlegs. 573 • The observation-based physically complete that accounts for the covariation of and 574 has a robust decreasing trend from cloud base to cloud top, which can be explained 575 by the increasing trend of the and correlation from cloud base to cloud top. 576 • The estimated by assuming a monomodal bi-variate lognormal joint distribution 577 between and systematically overestimates the observation-based results, 578 especially for the hlegs with a bimodal and joint distribution. The omission of the 579 variation and its covariation with tends to lead to an overestimation of EF despite 580 the error cancellation. 581 These results provide the following two new understandings of the EF for the autoconversion 582 parameterization that have potentially important implications for GCM. First, our study indicates 583 that the physically complete has a robust decreasing trend from cloud base to cloud top. Because 584 the autoconversion process is most important at the cloud top, this vertical dependence of EF 585 should be taken into consideration in the GCM parametrization scheme. Second, our study 586 indicates that effect of the and correlation plays a critical role in determining the EF. Lately 587 a few novel modeling techniques have been developed to provide the coarse resolution GCMs 588 information of subgrid cloud variation, such as the PDF-based higher-order turbulence closure 589 method-Cloud Layer Unified By Binormals, CLUBB (Golaz et al., 2002;Guo et al., 2015;590 Larson et al., 2002). These models are able to provide parameterized subgrid variance of which 591 can be used in turn to estimate . However, as shown in our study the tends to overestimate 592 the EF.