The moist halo region around shallow cumulus clouds in large eddy simulations

In this study, the moist buffering halo region of shallow maritime cumulus clouds is systematically investigated using large eddy simulations with various grid resolutions and numerical choices. Autocorrelation analyses of cloud liquid water and relative humidity suggest a converged size of 200–300 m for moist patches outside clouds when the model resolution is below 50 m, but may overestimate this size due to noncloudy moist regions. Based on a composite analysis, the structure of the moist halo immediately outside individual clouds is examined. It is found that, regardless of model resolution, the distribution of relative humidity in the halo region does not depend on cloud size, but on the real distance away from the cloud boundary, indicating some size‐independent length scales are responsible for the halo formation. The relative humidity decays with distance more quickly with finer horizontal resolution, which is possibly related to the model resolution dependence of the cloud spectrum. The halo size near the cloud base is larger than that within the cloud layer and this feature is robust across all simulations. Further analyses of backward and forward Lagrangian trajectories originating from the moist halo region reveal the possible role for subcloud coherent structures in cloud‐base halo formation. Possible mechanisms explaining cloud halo sizes and associated length scales are discussed.

environment, while in fact only the near-cloud environment air is mixed into the cloud.The underestimation of the specific humidity of the entraining air leads to smaller entrainment rates being diagnosed compared with direct estimations of entrainment rate using cloud properties in the halo region (Dawe & Austin, 2011;Romps, 2010).Hence better understanding of the moist halo region can help define the correct properties of the entraining air in a plume model of convection parameterization.
Besides dynamical effects, the higher relative humidity in the moist halo region also favours hygroscopic growth of aerosol (Carrico et al., 2003;Feingold & Morley, 2003;Flores et al., 2012;Petters & Kreidenweis, 2007).With higher aerosol concentration, the humidity in the halo region can be increased through mixing of more condensed water into the near-cloud environment and in turn can promote large-scale ascent and stronger convection (Abbott & Cronin, 2021).Aerosol humidification can also lead to a change of optical properties in the near-cloud environment (Altaratz et al., 2013).The gradual decrease of aerosol optical depth from cloud to clear sky in the "twilight zone" (Koren et al., 2007(Koren et al., , 2009)), a transition zone between cloud and cloud-free atmospheres, can have a non-negligible contribution to radiative forcing (Bar-Or et al., 2012;Eytan et al., 2020;Jahani et al., 2020).If such radiative effects of the moist halo region are neglected, remote sensing retrieval algorithms of aerosol properties can be biased toward data far from clouds and can lead to the underestimation of aerosol optical depth and possible uncertainties in radiative forcing associated with aerosol (Koren et al., 2007;Marshak et al., 2021;Mieslinger et al., 2021).Hence, the distribution of relative humidity is critical for estimating the aerosol humidification and the distribution of aerosol optical depth.
Therefore, characterizing the distribution of relative humidity in the halo region and the size of this region, and hence the correct representation of mixing in the halo region, can help advance the development of convection parameterization and improve the accuracy of remote sensing near cloud, shedding light on cloud dynamics, as well as cloud-aerosol-environment interaction.Nevertheless, there are disagreements regarding the moist halo region between theories, observations, and numerical simulations, partly due to different definitions of the cloud halo region.Theoretical studies (Pinsky & Khain, 2019, 2020) simplified the entrainment-mixing process at cloud boundaries using a one-dimensional turbulent diffusion equation and estimated the halo size to be around 100 m.However, observational studies have recorded a large uncertainty in halo size, ranging from less than 100 m to more than 1 km (Laird, 2005;Lu et al., 2003;Perry & Hobbs, 1996;Twohy et al., 2009;Wang & Geerts, 2010).A few high-resolution numerical simulations have been performed to investigate the halo region.Using large eddy simulations, Bar-Or et al. (2012) reported the characteristic scale of exponential decay of relative humidity to be slightly less than 100 m, and Lu et al. (2002) found a dependence of halo size on cloud size, but their horizontal resolutions were rather coarse (100-m grid length).Nair et al. (2021) investigated the interfaces at the edge of cumulus clouds using a direct numerical simulation, but this covered a small region of the cloud edge and could not provide comprehensive information on the halo region.Nair et al. (2021) also performed a high-resolution large eddy simulation with 4.1-m grid length and found that the size of the "invisible shell" is less than 200 m for a shell defined in terms of enstrophy.Heus et al. (2008) performed simulations of shallow cumulus clouds with grid lengths from 12.5-100 m, but they focused mainly on the downdraft shells, which have been found to be wider than the moist halo region (McMichael et al., 2022).The downward mass flux in cloud shells was stronger in finer resolution simulations (Heus et al., 2008) and the integrated mass flux in cloud shells was stronger for larger size clouds (Heus & Jonker, 2008).However, it remains unclear whether the properties of cloud shells can be applied robustly to understand the moist halo region, since we lack a systematic assessment of the sensitivity of moist halo structure to resolution and numerical choices using large eddy simulations.
The present study is designed to investigate the moist halo region around shallow cumulus clouds systematically, including the relative humidity distribution, the halo size, and possible physical processes involved in its formation, using high-resolution large eddy simulations.The rest of the article is organized as follows.Section 2 introduces the large eddy simulations (Section 2.1) and a composite algorithm for determining the relative humidity distribution within the halo region (Section 2.2).Section 3 examines the size of moist patches outside the cloud through autocorrelation analyses.Section 4 investigates general features of the relative humidity distribution within the halo region (Section 4.1), their dependence on model resolution (Section 4.2), and numerical details (Section 4.3).Section 5 reveals connections between the halo regions at different levels, by means of Lagrangian trajectories.A discussion is given in Section 6 and a summary in Section 7.

Large eddy simulations
The Met Office-Natural Environment Research Council (NERC) cloud model (MONC; Brown et al., 2015Brown et al., , 2018) ) is used to perform large eddy simulations of oceanic shallow convection based on the Barbados Oceanographic and Meteorological Experiment (BOMEX).Most of the model configuration follows that of Siebesma et al. (2003), but the grid spacing is changed.The horizontal grid spacings used are 100, 50, 25, and 10 m, in order to investigate the dependence of halo region structure on model resolution.Vertical grid spacings are 40, 25, 25, and 10 m, respectively.All simulations have the same model top at 3 km but the domain sizes are different, with consistent horizontal grids (600 × 600) to save computational resource.The 3D Smagorinsky-Lilly scheme is used for the parameterization of subgrid turbulence (Lilly, 1962;Smagorinsky, 1963).A simple saturation adjustment cloud scheme is used to represent the conversion between water vapour and cloud liquid water.There is no rain formation during our simulation period.In all the simulations, constant surface sensible and latent heat fluxes are prescribed.Rather than interactive radiation, we prescribe the large-scale radiative cooling to represent clear-sky longwave radiation.The radiative cooling is constant (−2 K⋅day −1 ) from the surface to 1.5 km height and decreases linearly to zero at the model top.To close the energy budget, we also prescribe a large-scale subsidence that increases linearly with height up to the inversion at 1500 m, above which it decreases.The subsidence is applied to both moisture and temperature fields.We further prescribe a small moisture tendency in the lowest 500 m to mimic the large-scale horizontal advection.The effects of large-scale pressure gradients are parameterized through imposed geostrophic winds (v g = (−10 + 1.8 × 10 −3 z, 0) m⋅s −1 ) and the Coriolis parameter f = 0.376 × 10 −4 s −1 .Other details of the case specification are available in Siebesma et al. (2003).Our analyses cover a period in the equilibrium state (hour 5-6) of the simulation, with 1 min output frequency.Consistent with the previous intercomparison study of Siebesma et al. (2003), the domain-averaged cloud properties remain steady during this period and thus are suitable for our analyses.

Composite algorithm
We use a spatial composite analysis, the "Onion Algorithm", to examine the distribution of relative humidity in the near environment around each cloud.At each vertical level, all cloudy points are first identified with the cloud liquid water criterion q l > 10 −5 kg⋅kg −1 .Contiguous cloudy points are combined to form an individual cloud object.For each cloud object, we identify its boundary and then investigate the distribution of relative humidity in the near-cloud environment as a function of distance from the cloud edge.Distances away from the edge are measured in terms of the real distance and also the distance normalized by cloud size.For the distributions in terms of real distance, we move outward from the cloud boundary in steps of a single grid box (Figure 1a).For the distributions in terms of normalized distance, at each vertical level we first calculate the effective radius of each cloud object as √ S∕, where S is the area coverage of the cloud object.We then express the radius as a number of grid points.The distribution is evaluated by moving outwards by this number of grid boxes at each step (Figure 1b).Any cloudy points outside the individual cloud in question and that are found during the outward movement are excluded from the composite.Mean properties for a given distance are composited to obtain the distribution in the halo region.Previous studies (Dawe & Austin, 2011;Zhao & Austin, 2005) applied similar ideas to understand the interaction between clouds and environment but were limited to the region adjacent to the cloud edge and thus are not able to cover the whole halo region.

SIZE OF MOIST PATCHES OUTSIDE THE CLOUDS
The size of moist patches outside the clouds is first examined using the spatial autocorrelation functions of relative humidity and cloud liquid water at each vertical level.The spatial autocorrelation function C(R) of a field f is defined as where r is the position vector in the field, R is the displacement position vector, and f * (r) represents the complex conjugate of f (r).The autocorrelation function can be computed with two fast Fourier transforms according to the Wiener-Khinchin theorem.Figure 2 shows the autocorrelation function of relative humidity at different levels.
Physically, the autocorrelation of relative humidity characterizes how the moist patches associated with coherent structures decay with distance.The spatial pattern of large correlation coefficients is found to be elongated along the west-east direction (Figure 2), and takes a more elliptical shape in the subcloud layer (Figure 2a).This is because the morphology of coherent structures is shaped by the east-to-west mean flow, which is largest (10 m⋅s −1 ) in the subcloud layer (Denby et al., 2022).The spatial patterns of the autocorrelation field of cloud liquid water from cloud base and above are closer to a round shape and similar across different vertical levels, consistent with the geometry of the clouds (Figure 3).In addition, the high autocorrelation coefficients of q l are more concentrated near the centre than those of relative humidity, indicating

F I G U R E 3
The same as Figure 2, but for the autocorrelation field of cloud liquid water q l .
that the clouds have more compact structures than the moist region.The autocorrelation of cloud liquid water has similar patterns near and above cloud base, except that the autocorrelation coefficient decays more quickly from the centre than the autocorrelation coefficient of the relative humidity field.Therefore, the sizes of moist patches are larger than the cloud sizes.We define the autocorrelation length scales L RH and L q l as the effective length scales of an enclosed area of the corresponding spatial autocorrelation fields as follows: where A is the area within which the autocorrelation coefficient is larger than e −1 .L RH and L q l can be considered as proxies for the sizes of moist patches and cloud objects, respectively.
Figure 4a shows the time-averaged (5-6 h) vertical profiles of L RH and L q l in the simulations at different resolutions.L RH is clearly larger than L q l at all vertical levels in each simulation.Both L RH and L q l start to converge at 25-m resolution, and the length scales in the 100-m simulation are much larger (about twice) than in the higher resolution simulations.In all simulations, L q l increases quickly with height near cloud base and is then fairly constant throughout the cloud layer.L RH is relatively small near the surface, where the size of turbulent eddies is constrained.It has a local maximum at around 100-150 m height, and decreases through the rest of the subcloud layer and through cloud base to achieve a local minimum at around 1000 m height.Thereafter, it increases again to the cloud top.A slight oscillation of L RH above 1000 m in the 10-m grid length simulation is probably due to a lack of sufficient sampling within a small domain size.Larger L RH in the upper part of the cloud layer might be related to terminal detrainment of moist air out of clouds.Moist patches may be large even if the corresponding clouds have dissipated, since their associated water vapour remains within the vicinity for longer than the cloud lifetime.The difference between L RH and L q l (ΔL = L RH − L q l ) provides a measure of bulk halo size in the autocorrelation field.Since the vertical variation of ΔL is largely controlled by L RH , we can examine how the halo sizes at different vertical levels are connected through a correlation analysis.Figure 4c shows the correlation coefficients between the time series of L RH at different vertical levels during hour 5-6 in the 25-m resolution simulation.The results from other simulations are similar (not shown).As expected, ΔL at a specified level is always highly correlated with that at neighbouring levels.Away from neighbouring levels, high positive correlations are also found at low levels between 250 and 750 m, and at high levels between 1500 and 2000 m.This indicates that the halo region near cloud base may be related to coherent structures in the subcloud layer, and that the halo region in the inversion layer may be associated with overturning structures near the cloud top.It is also found that ΔL at around 1000-1200 m is positively correlated with that in the inversion layer (1500-2000 m).Such a connection between the halo region in the mid-levels of the cloud layer and that at cloud top may indicate a role for downdrafts outside the cloud.Negative correlations between halo sizes at 500-1000 m and those at 1000-1500 m suggest a possible out-of-phase evolution, meaning that an increase of L RH in the mid-levels of the cloud layer is accompanied by a decrease of L RH in the inversion layer and vice versa.We hypothesize that the halo size from cloud top to the mid-levels of the cloud layer is increased due to the enhanced mixing between cloud and environmental dry air.Such mixing results in more negative buoyancy and thus leads to stronger downdrafts that can bring drier air from higher levels downward and decrease the size of the halo region below the mid-level of the cloud layer.

General features
The autocorrelation analyses above might overestimate the actual halo size, because some moist patches are remnants of dissipated clouds without any clouds within them.
To focus directly on the near environment around each cloud, we use the "Onion Algorithm" to assess the distribution of relative humidity away from the cloud edge (Section 2.2). Figure 5 shows the distribution of relative humidity perturbations (relative to the domain mean) outside the cloud in the 25-m grid length simulation at three vertical levels: 600, 1000, and 1500 m, which are representative of the cloud base, cloud layer, and near-cloud top, respectively.Only those cloud objects larger than 100 m are included in the composite analyses.These retained cloud objects are categorized into two groups: large and small, based on the median effective size (220 m near cloud base).The distribution expressed in terms of normalized cloud size shows clear differences between the larger and smaller clouds (Figure 5a,c,e).At all vertical levels, the relative humidity of large clouds decreases much more quickly to match the environment than that of small clouds.In contrast, the distributions expressed as a function of real distance are much more similar for the larger and smaller clouds (Figure 5b,d,f).The same observations can also be made for simulations at other horizontal resolutions (not shown).Hence, the decay of relative humidity within the halo region around shallow cumulus clouds scales better with real distance from cloud edge, indicating that the halo size is determined by some length scale or scales independent of cloud size.Some observational studies previously suggested that the halo size was proportional to the cloud size, but may have lacked sufficient sampling, or they focused on different types of clouds (Lu et al., 2003;Wang & Geerts, 2010).
Although the distributions for larger and smaller clouds are more similar when expressed in terms of real distance from the cloud edge, nonetheless the relative humidity around larger clouds at a given distance is lower than that around smaller clouds.This is consistent with the notion that larger clouds have stronger downdrafts, which in turn lead to a slightly drier halo region (Heus & Jonker, 2008;Rodts et al., 2003;Wang et al., 2009;Gu et al., 2020).This point is more apparent in the simulations with finer resolution and near the cloud top, because the cloud-top downdrafts are much better resolved with higher horizontal resolution.

Dependence on model resolution
As shown by Figure 6, it is important to notice that the distribution of relative humidity in the halo region is affected by the horizontal resolution.The relative humidity decreases more slowly from the cloud edge in the coarser resolution simulations, probably because the full spectrum of eddies responsible for mixing across the edge is less well captured.The decrease of relative humidity in the highest resolution simulation (10-m grid length) resembles an exponential decay, while the shape follows a more quadratic decay at lower resolutions.In other words, the distributions of relative humidity away from the cloud edge have not converged with increasing horizontal resolution, at least above 10-m grid length.Nonetheless, the decay rate of relative humidity is consistently found to be slower near cloud base (Figure 6a,d) than within the cloud layer (Figure 6b,c,e,f), indicating that the formation of the halo region near cloud base and at other vertical levels may be affected by different processes.We discuss this point further in Section 6.
If the outer edge of the halo region is defined as the position where the composited mean relative humidity perturbation approaches zero, then the halo size can be calculated as the distance between the cloud boundary and the outer edge.With this definition, we find that the halo sizes in the 10-, 25-, and 50-m simulations are comparable despite their different decay rates near the cloud edge.In each simulation, the halo size near cloud base is around 200 m and it decreases to around 100 m at higher levels.However, the halo size thus diagnosed is larger in the 100-m simulation at all vertical levels.A robust feature of all simulations is that the halo size is largest near cloud base and smaller within the cloud layer.This is also consistent with the results from autocorrelation analyses, apart from the impact of moist patches left by decaying clouds at levels around the cloud top.Similar vertical variation can also be found for downdraft cloud shells (Jonker et al., 2008).
However, the halo size is sensitive to how we define the outer boundary of the halo region.If a nonzero threshold of the relative humidity perturbation is used, then the halo size is smaller and also dependent on the horizontal resolution.The halo size becomes a monotonic function of horizontal resolution, with finer resolution simulations having smaller halo size due to the more rapid decay of relative humidity.The halo size does not converge within the range of resolutions explored in this study.The explanation for this resolution dependence of halo size may be related to the resolution dependence of cloud number density.Assume we have two large eddy simulations.The model grid lengths are Δx 1 and Δx 2 and Δx 2 < Δx 1 .The mean sizes of cloud objects at a specified vertical level are l c1 and l c2 .The mean sizes of moist regions in the two simulations are l m1 and l m2 .The numbers of clouds across the domain are N 1 and N 2 , respectively.A key result in our simulations, shown by Figure 7a,b, is that the fractional area coverage of cloud and halo regions (defined as the region with relative humidity perturbation larger than one standard deviation outside the clouds) are both independent of model resolution (see the proof in the Appendix).This implies the following equalities: Equation ( 4) can be rewritten as (5) Define L h1 = l m1 − l c1 and L h2 = l m2 − l c2 .L h1 and L h2 can be considered as the sizes of cloud halo regions when the model grid lengths are Δx 1 and Δx 2 , respectively.From Equation ( 5), we can derive the ratio between L h1 and L h2 : Combining Equations ( 3) and (4), we have and therefore Substituting Equation (8) (l c2 = l c1 l m2 ∕l m1 ) into Equation ( 6), the ratio between L h1 and L h2 is Shallow cumulus clouds in our large eddy simulations tend to be smaller and more numerous with increased horizontal resolution (Figure 7c).Similar behaviour can also be found in Brown (1999).Hence, we have N 2 > N 1 .As a result, the ratio L h1 ∕L h2 > 1 from Equation ( 9).This means that the mean size of the moist area around an individual cloud must be smaller in finer resolution simulations.

Sensitivity to numerical choices
It is plausible to speculate that the distribution of relative humidity may be sensitive to the numerical details of the model.The robustness of the composited structure in the halo region is therefore also examined with another large eddy model, the CM1 model (Bryan & Fritsch, 2002).The BOMEX simulations were again performed using horizontal grid lengths of 100, 50, 25, and 10 m, but with a smaller domain size (6.4 km) for computational considerations.Similar features can also be found in these simulations.
The distribution of relative humidity in the halo depends only weakly on the cloud size for a given simulation.Also, the rate of decay of the relative humidity perturbation is larger in the finer resolution simulations and smaller near cloud base (Figure 6d-f).
To test whether the size of the halo region is sensitive to details of the subgrid turbulent schemes (e.g., mixing length scale) or advection schemes, we perform additional sensitivity simulations at 25-m grid spacing.The mixing length scale in the subgrid turbulence scheme in MONC simulations is changed by setting the Smargorinsky constant C s from its default value 0.23 to smaller ones, 0.15 and 0.10.As the MONC model does not have multiple options for advection schemes, we test the sensitivity to the advection scheme using the CM1 model.The advection scheme in the control simulation with CM1 is the third-order WENO scheme (Balsara & Shu, 2000;Jiang & Shu, 1996).We further use the fifth, seventh, and ninth-order WENO schemes for the sensitivity simulations.Figure 8 shows that the general features found in control simulations are not sensitive to the numerical choices.

LAGRANGIAN TRAJECTORIES ANALYSIS
The two independent methods of Sections 3 and 2.2 give some consistent results in terms of the vertical variation of the moist halo region, but they cannot provide a picture of time evolution of air within the halo region.To understand further how the halo regions at different vertical levels are connected, and the physical processes involved, Lagrangian particles are used to trace the air parcels in the halo region (defined as RH ′ >  RH , where  RH is one standard deviation of relative humidity) outside the cloud at all vertical levels and at each model output time during hours 5-6 (1 min interval).The Lagrangian trajectories are calculated following the method of Gheusi and Stein (2002), with some extensions.The positions (coordinates) of model grid boxes are used as Lagrangian labels and are advected with the flow using the same advection scheme as that applied to the scalar fields in the model.The trajectories of labelled particles can then be calculated backward and forward through the advected coordinates.The trajectories for each model output time are calculated both backward and forward for 30 min.We chose the 60-min time window, as it is longer than the entire lifetime of almost all clouds in our simulations.
The particles in the moist halo region at the reference times come from other parts of the domain and thereby are located at different heights before and after the formation of the halo region.Figure 9 shows the distributions of heights of Lagrangian trajectories before (−30 min, −10 min) and after (10 min) the reference times and it can be used to indicate the neighbouring levels that are critical during the formation of a moist halo region.Near cloud base (Figure 9a), 30 min before the reference time, slightly more than 50% of the air parcels in the halo region come from the neighbouring levels (about 250 m below and above).However, the other half of the air parcels originate from the subcloud layer, with most of them being near the surface (Figure 9a).Ten min after the formation of the halo region, about 70% of the air parcels have moved downward and half of them (35% of the total) go back to the subcloud layer.These findings provide clear evidence that the halo region near the cloud base is closely related to coherent structures from the subcloud layer.More than half of the air parcels within the halo region in the middle of the cloud layer (1000 m, Figure 9b) and near the cloud top (1500 m, Figure 9c) come from higher levels, and they descend slowly to form the halo.However, only 10 min after the reference time, more than 65% of the air parcels have already descended to lower levels, suggesting that the formation of the halo region is accompanied by a downdraft (Heus & Jonker, 2008;McMichael et al., 2022).These results provide evidence to support our hypothesis of length scales associated with the moist halo region in the next section.

DISCUSSION
The region with downward motion outside the cloud is usually referred to as a "cloud shell," but it is not necessarily related to higher water vapour (Savre, 2021).Recent studies (McMichael et al., 2022;Savre, 2021) suggested that, from the composited perspective, the region with downward motion outside the cloud is broader than the halo region with higher water vapour.Thus, the moist halo region seems to be a subset of the cloud shell, and it should be emphasized that the moist halo region investigated in this study is not the same as the downdraft cloud shells studied by Jonker et al. (2008); Heus and Jonker (2008); Heus et al. (2008), for example.First of all, the primary formation mechanisms of the moist halo region and the cloud shell are different.Since the large-scale relative humidity and moisture content decrease with height in the simulations, the descending cloud shell alone would result in a drier near-cloud environment outside the cloud, which is not the case.The presence of a moist halo region immediately outside the cloud is thus strong evidence that horizontal mixing occurs near cloud boundaries.The mixing between the detrained cloud condensate and the environmental air leads to evaporation and humidifies the near-cloud environment.Meanwhile, the evaporative cooling starts to drive downward motions and thus the formation of the cloud shell.In this sense, the moist halo region and cloud shell form simultaneously, but the underlying mechanisms are not quite the same.
In addition, the moist halo region always surrounds each cloud object, while the strong downdrafts within the cloud shell are not necessarily present, as shown in Figure 10.The distribution of strong downdrafts outside the cloud also has stronger asymmetry, compared with the moist halo region, probably because of the weak vertical wind shear.Savre (2021) found that, in addition to the buoyancy effect, other mechanical forcings, for example, the pressure gradient force and the horizontal As in (a) but including overlapping downdrafts.The clouds are defined using q l >10 −5 kg⋅kg −1 .The moist halo region is defined as where the relative humidity anomaly is larger than one standard deviation of relative humidity at 600 m.The downdrafts are defined as the region with downward motion stronger than one deviation of vertical velocity at 600 m. advection, may be important for downward motion in the cloud shell.These results indicates that there might be more dynamical processes involved in the formation and maintenance of the cloud shell, which contribute to the asymmetries.Furthermore, in terms of detailed structures, Heus et al. (2008) found that the downward mass flux density was stronger in higher resolution simulations but the size of the downdraft shell was consistent across different grid spacings (their fig.10), which is in contrast with the resolution dependence of the moist halo region.Heus and Jonker (2008) showed that the integrated mass flux in cloud shells depends on cloud size, while our results suggest that the relative humidity distribution in the moist halo region scales with real distance from the cloud edge.These points strongly indicate that the moist halo region is different from the downdraft shell and worthy of in-depth understanding.
The fact that the distribution of relative humidity within the halo region scales better with the real distance away from the cloud edge rather than with cloud sizes indicates some size-independent length scales governing the formation of the halo region.A robust finding from all simulations is that the cloud halo size is largest near cloud base and decreases upwards.In considering this behaviour, assume that the largest overturning structure responsible for the mixing between cloud and environment has a length scale of l 0 .That structure breaks down continuously into smaller scales until the eddy is dissipated.We hypothesize that the halo size should be characterized by the mean size of these continuously breaking eddies.We estimate the mean size using the energy-weighted mean as where E(l) dl = E(k) dk is the energy spectrum at length l or wavenumber k and l K is the Kolmogorov length.Assuming that the energy spectrum follows the "−5/3" power law in the inertial range, we have Here we have used the fact that l K ≪ l 0 .We should keep in mind that the simulations cannot capture the full spectrum across the inertial range, because eddies with sizes smaller than the grid length cannot be resolved.Therefore, the factor proportional to the largest eddy size l 0 will be slightly larger than "2∕5" since fewer small eddies are resolved explicitly.The factor is only used for a rough estimation to have a comparison with our analyses.
As shown in Section 5, backward and forward trajectories of Lagrangian particles reveal a close connection of cloud-base halo formation with subcloud coherent structures.In the subcloud layer, a reasonable first guess for l 0 would be the height of the well-mixed subcloud layer.The mixed-layer height in the BOMEX case is around 500 m and thus we estimate l to be 200 m.This is consistent with both the autocorrelation and composite analyses.In the cloud layer, a reasonable length scale near clouds is the buoyancy length scale (Craig & Dörnbrack, 2008).The buoyancy length scale in our simulations can be estimated as √ e c ∕N, where e c is the turbulent kinetic energy (0.5(u ′ 2 + v ′ 2 + w ′ 2 )) in the cloud and N is the Brunt-Väisälä frequency.The buoyancy length scale describes the maximum vertical displacement that can be induced against the stratification in the environment by buoyancy-driven pressure perturbations and thus the maximum scale of eddies that cross the cloud boundary.The mean value of this buoyancy length scale in the cloud layer is around 150 m and thus results in a mean length scale of 60 m, which is smaller than that near cloud base.
Our large eddy simulations produce converged area fractions of cloud across different resolutions, indicating that the properties of the cloud field are controlled by the large-scale forcing (Brown, 1999;Craig, 1996).The converged area fraction of moist patches across different resolutions is a surprise.Possible reasons for the constancy of halo area fraction might be also related to the prescribed large-scale forcing, as discussed in the Appendix.However, the cloud spectrum changes with model resolution in our simulations, leading to a resolution dependence of the relative humidity distribution away from the cloud edge, as explained in Section 4. Thus, the lack of convergence in the relative humidity distribution in the halo region may be a numerical bias induced by the lack of convergence in cloud number.Whether the distributions converge at even higher resolutions needs further investigation.This may also raise doubt about the fidelity of large eddy models to capture the details of natural clouds realistically, so long as the cloud spectrum depends on resolution, when model grid length is no finer than 10 m.Although previous studies (Siebesma & Jonker, 2000) have shown that large eddy models can reproduce the fractal behaviour of clouds (area-perimeter fractal dimension) reasonably, the distribution of relative humidity changing with horizontal resolution suggests that aspects of the detailed cloud morphology may still be difficult to capture.A recent study found that, in comparison with observations, large eddy models tend to generate more plume-like, rather than bubble-like, clouds (Romps et al., 2021).These results indicate a continuing need for improvement of large eddy models to capture detailed structures associated with cloud geometry better.

SUMMARY
The moist halo region, immediately outside a cloud, is moister than the air further from the cloud and is different from the cloud downdraft shell.It is critical for the interplay between the cloud and the large-scale environment and also has non-negligible impact on radiation.In the present study, we systematically investigated the halo region using large eddy simulations across various model resolutions.Autocorrelation analyses of cloud liquid water and the relative humidity field revealed the converged size of moist patches outside the cloud to be around 200-300 m when the model spacing is below 50 m.This value may overestimate the size of the halo region due to the presence of moist patches left by dissipated clouds.To focus on the structure around individual clouds, we examine the distribution of relative humidity from the cloud edge based on an "Onion algorithm".In contrast to previous studies (Lu et al., 2002;Wang et al., 2009), the distribution of relative humidity in the halo region is independent of cloud size and scales much better with real distance away from the cloud boundary, indicating some size-independent length scales responsible for its formation.However, the distribution of relative humidity depends strongly on model grid spacings, with larger decay rates in higher resolution simulations, leading to smaller halo sizes.This may be related to the inability of the large eddy model to simulate a consistent cloud spectrum across the range of model resolutions explored in this study.Nevertheless, regardless of grid spacings, a robust feature is that the cloud halo size varies vertically, with the largest halo near the cloud base.Lagrangian trajectory analyses suggest that the formation of the halo region at different vertical levels may result from different physical processes.The size of the halo region in the cloud layer is possibly affected by the buoyancy length scale.The halo region near the cloud base is likely related to coherent structures in the subcloud layer and thus is characterized by the depth of the mixed layer.Finally, we want to stress that this study focused only on the halo region outside nonprecipitating shallow cumulus clouds.Whether the conclusions or the physical processes can be applied to understand the halo region of organized convection or deep convection in response to different large-scale forcings, for example, or over different basins or continents, remains unclear.Such studies have larger computational demands and need further investigation.It should also be noted that aerosol impacts were not considered in our simulations, although their role has been discussed in the Introduction.How aerosol-cloud interactions may affect the dynamics near the cloud edge and the stratification through vertical-dependent radiative effects, and thus change the size of the halo region, is also left for future studies.face energy fluxes, together with the prescribed subsidence warming, are in equilibrium with the prescribed radiative cooling so that the whole simulated domain achieves energy balance at the equilibrium period.Because no precipitation occurs in the BOMEX case, there should not be net heating at any vertical level and a steady state can be reached.If simulations at different resolutions achieve a very similar steady state, then we might plausibly expect the evaporative cooling contribution to the energy budget to be consistent with resolution.We know that evaporative cooling occurs predominantly within the moist halo region, where there is mixing between cloud and the environmental air.If we can further assume that the moist halo area fraction controls the total evaporative cooling, then it follows that  h should remain constant when the resolution is changed.

F
I G U R E 1 Schematic diagram of the algorithm to detect the near-cloud environment step-by-step in terms of (a) real distance and (b) normalized distance, outward from the edge of each cloud object.The grey shading represents an example of a cloud object.In (a), cyan, yellow, green, red, blue, magenta, and brown colours represent the environment 1, 2, 3, 4, 5, 6, and 7 grid boxes away from the cloud boundary, respectively.Similarly, in (b), these colours denote the environment 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, and 6.5 times the cloud size (R) away from the cloud boundary, respectively.R is the effective radius of each cloud object R = √ S∕, where S is the area coverage of the cloud object.F I G U R E 2 Autocorrelation field of relative humidity RH in a 25-m grid length simulation at different vertical levels: (a) 250, (b) 600, (c) 1000, and (d) 1500 m.The white contour represents the e-folding line.
Figure 4b shows the vertical profile of ΔL.The halo sizes in the 10-and 25-m simulations are comparable (200-300 m) throughout the cloud layer, while those in the 50-m simulation are somewhat larger, particularly in the upper part F I G U R E 4 Time-averaged (5-6 h) (a) vertical profiles of autocorrelation length scales for relative humidity (L RH , solid lines) and cloud liquid water (L q l , dashed lines) and (b) vertical profiles of halo sizes (L RH − L q l ) in simulations with different horizontal grid lengths: 10 (blue), 25 (red), 50 (green), and 100 m (yellow).(c) Correlation coefficients between the time series (5-6 h) of the autocorrelation length scale of relative humidity at different vertical levels in the simulation with horizontal resolution of 25 m.Due to the symmetry, the lower half of the triangular correlation matrix is not shown.The coefficients are shown within a vertical range of 1000 m from the current level, because the air parcels that form the halo region do not travel more than 1000 m in the vertical, as shown in Figure 9. of the cloud layer.Halo sizes in the 100-m simulation are much larger.

F
I G U R E 5 The composited distributions (perturbations have been interpolated in 10-m intervals before being composited) of relative humidity perturbations as functions of (a, c, e) normalized distance and (b, d, f) real distance outward from the cloud boundary, at (a, b) 600 m, (c, d) 1000 m, and (e, f) 1500 m heights in the 25-m grid length simulation.Large red dots are composites for clouds with radii larger than the median value, while blue small dots are composites for smaller clouds.

F
The composited distributions of relative humidity perturbations as functions of real distance from the cloud boundary, at heights of (a, d) 600 m, (b, e) 1000 m, and (c, f) 1500 m.The left (a, b, c) and right columns (d, e, f) show results from MONC and the CM1 model, respectively.Different horizontal grid lengths are represented with different colours: 10 (blue), 25 (red), 50 (green), and 100 m (yellow).

F
I G U R E 7 (a) Vertical profiles of cloud area fraction in different resolution simulations.(b) Vertical profiles of the area fraction of the halo region outside the clouds in different resolution simulations.The inner boundary of the halo region is defined as the cloud edge and the outer boundary is defined using one standard deviation of the relative humidity perturbation at each vertical level.(c) Vertical profiles of cloud number density ((km 2 ) −1 ) in simulations with different horizontal resolutions.The solid blue, red, green, and yellow lines represent the results from simulations with grid lengths of 10, 25, 50, and 100 m, respectively.

F
The composited distributions of the relative humidity perturbation as functions of real distance from the cloud boundary at (a, d) 600 m, (b, e) 1000 m, and (c, f) 1500 m heights from 25-m grid length simulations.The left column (a-c) shows the results in MONC simulations with different settings of mixing length scale in the subgrid turbulence scheme: C s =0.23 (blue), 0.15 (yellow), and 0.10 (cyan).The right column (d-f) shows the results in CM1 simulations with different orders of WENO advection scheme: third (blue), fifth (red), seventh (green), and ninth (yellow).

F
Probabilitydistributions of the heights of Lagrangian trajectories in the 10-m grid length simulation.Trajectories are calculated for air parcels that form the halo region at the reference times, and different colours represent the distribution at different times relative to the reference time: −30 min (blue), −10 min (red), and 10 min (cyan).The different panels are for the halo regions defined at different vertical levels at reference times of (a) 600, (b) 1000, and (c) 1500 m.The orange dot in (a) denotes the height of the cloud base.

F
I G U R E 10 Snapshot of regions of cloud (red shading), moist halo (blue shading), and downdrafts (yellow shading) at hour 6 and at a height of 600 m.(a) Snapshot showing the cloud and moist halo regions.(b)