Representing surface heterogeneity in land-atmosphere coupling in E3SMv1 single-column model over ARM SGP during summertime

. The Earth’s land surface features spatial and temporal heterogeneity over a wide range of scales below those resolved by current Earth system models. State-of-the-art land and atmosphere models employ parameterizations to represent their subgrid heterogeneity, but the land-atmosphere coupling in ESMs typically operates on the grid scale. Communicating the information of the land surface heterogeneity with the overlying atmospheric boundary layer (ABL) remains a challenge in modeling land-atmosphere interactions. In order to account for the subgrid scale heterogeneity in land-atmosphere coupling, we 5 implement a new coupling scheme in the Energy Exascale Earth System Model version 1 (E3SMv1) that uses adjusted surface variances and covariance of potential temperature and specific water content as the lower boundary condition for the atmosphere model. The new lower boundary condition accounts for both the variability of individual subgrid land surface patches and the inter-patch variability. E3SMv1 single-column model simulations over the Atmospheric Radiation Measurement (ARM) Southern Great Plain (SGP) site were performed to assess the impacts. We find that the new coupling parameterization increases 10 the magnitude and diurnal cycle of the temperature variance and humidity variance in the lower ABL in non-precipitating days. The impacts are primarily attributed to subgrid inter-patch variability rather than variability of individual patches. These effects extend vertically from the surface to several levels in the lower ABL on clear days. We also find that accounting for surface heterogeneity increases low cloud cover and liquid water path. These cloud changes are associated with the change in cloud regime indicated by the skewness of the probability density function (PDF) of the subgrid vertical velocity. In precipitating 15 days, the inter-patch variability reduces significantly, so that the impact of accounting for surface heterogeneity vanishes. These results highlight the importance of accounting for subgrid heterogeneity in land-atmosphere coupling in next generation Earth system models. (8)) contributes to 70% of the increase. We also find that the covariance between potential temperature and specific humidity variance in HET is about 4.8 times larger than that in HOM, and that the inter-patch covariance (i.e., the first term in rhs of Eq. (9)) contributes to 79% of the increase. These results demonstrate that the inter-patch differences of potential temperature and specific humidity is the primary term driving the differences in surface properties shown in Fig. 3. approach considers the subgrid variability in surface states, improving the realism of the heterogeneous land-atmosphere interactions in E3SM. This study provides an initial assessment of this treatment in a state-of-the-art ESM in a SCM configuration. The significant effects of this treatment on turbulence structure in the ABL and on low-level clouds highlight the importance of accounting for the subgrid variability in land-atmosphere coupling. Future studies that assess the impact of the new treatment on the climate system in the fully coupled E3SM will be highly valuable.


Introduction
Land surface heterogeneity can influence the atmospheric boundary layer dynamics (Bou-Zeid et al., 2020), cloud population 20 (Fast et al., 2019;Chen et al., 2020), and rainfall initiation (Shrestha et al., 2014;Gao et al., 2021) on relatively small spatial and temporal scales. Representing the subgrid scale heterogeneity of land surface is essential in Earth system models (ESMs) because their typical horizontal grid spacing of about 100 km is insufficient to resolve the variability of land surface characteristics. A common method for representing this surface heterogeneity in land surface models (LSMs) is to introduce "patches" or "tiles" to delineate different land use types, soil characteristics, topographic characteristics (i.e., elevation, slope and aspect), 25 etc., in a model grid box (Koster et al., 2000;Lawrence et al., 2019). Multiple pre-defined patches (e.g., naturally vegetated land, human managed cropland, urban, and lake) may be present within any given grid cell to collectively summarize the land surface characteristics. Recently, emerging schemes begin to include lateral exchange of water, energy, momentum, and carbon among patches, accounting for the effects of those lateral exchanges on the grid scale features and responses to environmental changes (Chaney et al., 2016(Chaney et al., , 2018de Vrese et al., 2016;Lawrence et al., 2019). 30 Representing the subgrid heterogeneity in ESMs is not limited to the land component, but in the atmosphere component as well. In the U.S. Department of Energy's Energy Exascale Earth System Model version 1 (E3SMv1; Golaz et al., 2019), the E3SM Atmosphere Model version 1 (EAMv1; Rasch et al., 2019) utilizes Cloud Layers Unified By Binormals (CLUBB) Larson et al., 2002;Larson and Golaz, 2005;Bogenschutz et al., 2013) for sub-grid atmospheric heterogeneity in turbulent mixing, shallow convection, and cloud macrophysics. In addition, CLUBB predicts higher-order moments 35 required for the closure of the model equation system. As documented in Rasch et al (2019), EAMv1 produces improved climatology compared to the previous-generation ESMs, even though the model exhibits common and long-standing cloud biases that might be related to deficiencies in the cloud, turbulence, and convection parameterizations (Xie et al., 2018;Brunke et al., 2019;Zhang et al., 2019).
While both the land and atmosphere models employ parameterizations to represent their spatial heterogeneity, the two 40 component models in E3SMv1 communicate at grid scale. Only grid-box averages of land and atmospheric states and fluxes are exchanged between the land and atmosphere models. Neglecting the subgrid information of the land and atmosphere in their coupling potentially leads to simulation biases (Manrique-Suñén et al., 2013;Clark et al., 2015). Previous studies have demonstrated the turbulent mixing process is significantly different when the subgrid variability is considered, impacting the simulated mean atmospheric state (Mahrt, 2000;Molod et al., 2003;de Vrese et al., 2016).

45
The main objectives of this study are to 1) implement a new land-atmosphere coupling scheme in a state-of-the-art ESM to account for the subgrid land-surface heterogeneity in land-atmosphere coupling, and 2) assess its impacts on atmospheric boundary layer characteristics, clouds, and precipitation simulations. This paper is organized as follows. Section 2 introduces the land-atmosphere coupling parameterizations and the single-column model configuration. Section 3 presents the comparison of the original coupling parameterization and the new treatment accounting for the subgrid scale heterogeneity. Section 4 gives 50 a summary of major findings.

Methodology
In this study, we implemented a new land-atmosphere coupling scheme in E3SMv1 that integrates the information on subgridscale heterogeneity to pass on from the E3SM Land Model (ELM) to EAM. The model was configured as a single-column

Land-atmosphere coupling in E3SM
In ELM, spatial land surface heterogeneity is represented as a nested subgrid hierarchy in which the highest level is the grid cell while the lowest level is referred to as patch. Over SGP the grid comprises eight patches as listed in Table 1. The patches are designated to capture the geophysical and biogeochemical differences between broad categories of land cover and plant 60 functional types within the grid cell. The mean states and fluxes to/from the surface are computed at each patch individually and then aggregated to the grid scale because EAM's turbulence parameterization CLUBB operates only on grid cell means.
In E3SM, the exchanges of mass and energy between land and atmosphere are determined using grid cell averages of surface properties to serve as the lower boundary conditions. Therefore, the subgrid information on land heterogeneity is lost in the coupling processes and the state of the atmosphere at the surface is treated as being spatially homogeneous within the grid cell.

65
This default treatment of land-atmosphere coupling is referred to as the "HOM" (homogeneous) method hereafter in this study.
In EAM, CLUBB assumes a multivariate, double Gaussian probability density function (PDF) to describe the variability of subgrid liquid water potential temperature (θ ′ l ), total specific water content (q ′ t ), and vertical velocity (w ′ ), with lower boundary conditions provided by the ELM. Because of the staggered grid configuration in CLUBB, only the second-and fourth-order moments are required as lower boundary conditions of the atmosphere . The fourth-order moment w ′4 is 70 computed as part of the PDF closure scheme at the surface, which will not be discussed in this study. The second-order moments include the turbulent fluxes of momentum (u ′ w ′ , v ′ w ′ ), heat (w ′ θ ′ l ) and moisture (w ′ q ′ t ), as well as the turbulent variances and covariance of liquid water potential temperature and total specific water content (θ ′2 l , q ′2 t , and θ ′ l q ′ t ). These scalar variances and covariance can be adjusted to account for the spatial heterogeneity that emerges over the land surface. Conventional surface layer similarity theory, aka Monin-Obukhov similarity theory (MOST), is used in CLUBB to link the surface scalar variances to the corresponding fluxes (adapted from André et al., 1978): Here the parameters c 1 , c 2 , and c 3 are set to 0.4, 0.4, and 0.2, respectively; w ′ θ ′ l and w ′ q ′ t are area-weighted fluxes of heat and 80 moisture, respectively, derived from ELM at grid level. Similarly, the surface velocity scale (u * ) is approximated based on the area-weighted momentum fluxes as follows: The convective velocity w * is defined as w * = g T0 z i w ′ θ ′ l 1 3 when surface heat flux is positive (i.e., unstable conditions), where g is the gravitational acceleration, T 0 is a reference temperature set to 300 K, z i is the mixed layer height in André et al.

85
(1978) but it is set to 1 m in CLUBB. Under stable or neutral conditions, w * = 0.
It is worth noting that Eq. (4) is insufficient to properly represent the characteristics of a heterogeneous underlying surface.
Since the MOST is applied for each patch in ELM, the aggregated grid cell fluxes neglect the subgrid scale heterogeneity in land-atmosphere coupling. Therefore, the homogeneous treatment of the friction velocity u * might not be representative of the land surface heterogeneity in the context of land-atmosphere interactions, potentially leading to biases in the lower atmosphere.

Implementation of the patch-aggregating approach in E3SM
To account for the surface heterogeneity in land-atmosphere coupling, Machulskaya and Mironov (2018) proposed an approach to aggregating the patch characteristics into equivalent grid-cell scalar variances and covariance. Given the MOST is applied for each patch, this representation thus can retain the subgrid-scale heterogeneity of land surface. In this study, we implemented the Machulskaya and Mironov (2018) scheme in E3SM as briefly described below.

95
The decomposition of a generic variable x into the grid-box mean ⟨x⟩, the fluctuation of patch-mean quantities from the grid-cell mean x ′′ , and the sub-patch fluctuation x s , can be expressed as: The angle brackets (⟨⟩) denote the average quantity over the grid cell of a host model and the overline (¯) denotes a mean state over the patch. The double prime ( ′′ ) denotes a deviation of the patch mean from the grid-cell mean. A schematic diagram 100 illustrating the decomposition is given by Fig. 1 from Machulskaya and Mironov (2018). Based on Eq. (5), we obtain the fluctuation about grid-box mean x ′ = x − ⟨x⟩ = x ′′ + x s . Assume ⟨x ′′ ⟩ = 0 by definition and x s = 0 to simplify the discussion, then we find where x, y represents either θ l or q t . The angle brackets denote the area-weighted average. For reference, we list the formulas 105 for computing the aggregated variances and variance of θ l and q t : Here f i represents the weight of patch i relative to the corresponding grid cell, and the f i s for every single grid cell sum to one.

110
Given the fact that the liquid water content q l is negligible at the surface, we adopt the surface specific humidity q and potential temperature θ for computations at individual patches.
In Eqs. (7-9), the first term in the right-hand side (rhs) represents the contribution of the aggregated grid variance from inter-patch difference, as it is the variance of the patch means with respect to the grid cell mean. The patch means are derived in ELM using MOST based on near-surface atmospheric state. The second term represents the contribution of the aggregated 115 grid variance from the variance of individual patches (i.e., sub-patch variance). MOST is applied for each patch (i.e., applying Eq. (1-3) for each patch) to link that sub-patch turbulent variances with their respective fluxes. This approach is able to retain the patch-level characteristics in the land-atmosphere coupling. A schematic diagram is given in Fig. 1 to illustrate the implementation of the new representation of the heterogeneity in the new land-atmosphere coupling. When aggregating the patch-specific moments into a grid cell, this approach accounts for both the "inter-patch" variability of surface states (i.e., the 120 first term in the rhs of Eq. (7)) and the variability of each patch ("patch" variance denoted by the second term in the rhs of Eq. (7)). Hence, the surface boundary conditions contain the information that characterizes heterogeneous surface fluxes and mean states.

E3SM SCM configurations
We configured the EAMv1 to run at SCM mode with 72 vertical levels over the ARM SGP site. The E3SM SCM has been 125 demonstrated to be a useful tool for model development at a low computational cost (Bogenschutz et al., 2020). We performed two simulations from June 1 to August 31 in 2015. The control simulation uses the default E3SM land-atmosphere coupling, which assumes homogeneous land surface, and is labeled as "HOM". The sensitivity simulation that uses the new land-atmosphere coupling described in Section 2.2 is labeled as "HET". Fig. 2 provides a schematic illustrating how the new land-atmosphere coupling is utilized in the E3SM SCM configuration. Using the spatially heterogeneous states and fluxes 130 provided by ELM at the patch level, we computed the variances and covariance of potential temperature and specific water content, and provided them to EAM as the lower boundary condition in CLUBB. A short-term hindcast approach was applied in order to ensure the model was well constrained by the large-scale conditions (Ma et al., 2015). During the simulation period, the SCM was initiated every day at 00:00 UTC (18:00 local time (LT)) and run for 48 h with prescribed large-scale forcing . The 24 to 48 h forecasts in each simulation were then combined into a continuous time series for analysis. 135 We note that discontinuity between two consecutive days is expected when using the hindcast approach. water potential temperature ⟨θ ′2 l ⟩, (b) total water specific humidity ⟨q ′2 t ⟩, and (c) their covariance ⟨θ ′ l q ′ t ⟩, which are prescribed based on the ELM output and then passed into CLUBB as the surface boundary conditions. The solid lines denote the seasonal mean quantities and correspondingly the filled areas indicate the 25% to 75% quartile over the period from June to August 2015.
3 Results Figure 3 shows the diurnal cycle of the surface properties averaged over the simulation period. The discontinuity caused by the hindcast approach is present at 18:00 LT, as expected. We find that the HET configuration produces significantly larger 140 aggregated variances of surface potential temperature (⟨θ ′2 l ⟩) and specific water content (⟨q ′2 t ⟩) than the HOM configuration ( Fig. 3a-b). This is expected because the HET configuration accounts for the land surface heterogeneity so that the variability of the two scalars increases. Fig. 3c shows that the influence of the HET approach on the temperature-humidity covariance (⟨θ ′ l q ′ t ⟩) depends on the time of the day. During the nighttime and early morning, ⟨θ ′ l q ′ t ⟩ in HET is significantly larger than that in HOM, which is near zero. The covariance continues to decrease in the morning and eventually becomes negative in the early 145 afternoon. Table 2 summarizes the means of the grid-aggregated surface quantities in HOM, HET, and the two components of the HET approach (i.e., the two rhs terms in Eq. (7-9)). We find that the potential temperature variance in HET is 5 times larger than that in HOM. Notably, the inter-patch variance (i.e., the first term in rhs of Eq. (7)) contributes to 82% of the increased gridaggregated variance. Similarly, the specific humidity variance in HET is also 5 times larger than that in HOM. The inter-patch 150 variance (i.e., the first term in rhs of Eq. (8)) contributes to 70% of the increase. We also find that the covariance between potential temperature and specific humidity variance in HET is about 4.8 times larger than that in HOM, and that the interpatch covariance (i.e., the first term in rhs of Eq. (9)) contributes to 79% of the increase. These results demonstrate that the inter-patch differences of potential temperature and specific humidity is the primary term driving the differences in surface properties shown in Fig. 3. Table 2. Patch means of the variances and covariance of temperature and humidity at the surface. This table lists the mean quantities averaged over the JJA season in 2015. The first column (HOM) is derived from the default CLUBB approach using the grid-cell average fluxes provided by ELM. The second column shows the HET-computed surface moments using the subgrid-scale fluxes and states from individual ELM patches. Following Eqs. (7-9), the addition of the last two columns is equivalent to value of the HET:total column.

HOM
HET: total HET: patch HET: inter-patch We next investigate the negative temperature-humidity covariance ⟨θ ′ l q ′ t ⟩ in Fig. 3c. In Fig. 4a that shows the inter-patch variability of surface potential temperature at 14 LT when the covariance is at its most negative, we find that the lake (patch #33) temperature is significantly lower than the grid-cell mean and contributes the most to the total temperature variance even though the lake patch only occupies a small fraction of the grid cell (Fig. 5a). In terms of the humidity variance, the vegetated bare ground (patch #0) turns into the dominant factor affecting the total humidity variance (Figs. 4b & 5b). Not only these two 160 outstanding patches but also other vegetated patches actually exhibit a negative correlation between temperature and moisture likely due to the rapid evaporation in the afternoon, leading to the negative grid-aggregated temperature-humidity covariance under convective boundary layer conditions.

Vertical structure of the atmospheric boundary layer in clear-sky days
In this section, we investigate the effects of surface heterogeneity on the ABL characteristics in clear-sky days, when the 165 maximum daytime (06-18 LT) mean liquid water content below 600 hPa is less than 10 −3 g kg −1 . In the simulation period, there are 17 clear-sky days for this analysis.  l , (d) specific humidity variance q ′2 t , and (e) temperature humidity covariance θ ′ l q ′ t averaged over 14-16 LT for the clear-sky days. For clarity, the overline for atmospheric quantities hereafter denotes grid cell mean.
The average ABL top on clear-sky days is located at around 830 hPa for SGP during the summer season of analysis (Figs. 6a-b). The simulated ABL structures in HOM and HET are very similar except that the HET configuration produces a slightly warmer ABL than HOM. In Figs. 6c-d, we find that the increase in the prescribed surface temperature variance and the surface 170 humidity variance in HET propagates upward through the lower half of the ABL. Fig. 6e shows that the decrease in temperaturehumidity covariance in HET also propagate up towards 920 hPa. The changes in θ ′2 l , q ′2 t and θ ′ l q ′ t can further affect the buoyancy term in CLUBB.
When no cloud is present, the buoyancy production term related to θ l can be written as: Figure 7. Same as Fig. 6, except for (a) the buoyancy production term of θ l , θ ′ l θ ′ v , (b) the turbulent flux of liquid water potential temperature w ′ θ ′ l , (c) the buoyancy flux w ′ θ ′ v , (d) the vertical velocity variance w ′2 , (e) the qt buoyancy term, q ′ t θ ′ v , (f) the turbulent flux of water specific content w ′ q ′ t , (g) the third-order moment of vertical velocity w ′3 , and (h) the skewness of the vertical velocity PDF Skw.
where the coefficient c 0 = 1−ϵ0 ϵ0 θ 0 is approximately 200 K. ϵ 0 = R d /R v , R d is the gas constant of dry air, R v is the gas constant of water vapor, and θ 0 is a reference temperature. According to Eq. (10), the decrease in θ ′ l q ′ t would counteract the increase in θ ′2 l for the buoyancy production of θ l in HET simulations. However, considering the order of magnitude of the θ ′2 l term and the θ ′ l q ′ t term in Fig. 6, we conclude that the buoyancy production term of θ l is dictated by the θ ′2 l term. Thus, the buoyancy term θ ′ l θ ′ v in HET extends vertically in the lower ABL (Fig. 7a). Furthermore, the increase in HET θ ′ l θ ′ v helps generate the turbulent 180 flux of temperature w ′ θ ′ l (Fig. 7b), contributing to the production of the upward buoyancy flux (Fig. 7c) and the turbulent kinetic energy (Fig. 7d) in the lower ABL. The enhanced temperature flux thus leads to enhanced vertical mixing and a slightly warmer ABL (Fig. 6a).
Similarly, the q t -related buoyancy term can be expressed as: 185 By comparing the magnitude of the two rhs terms in Eq. (11), we find that the buoyancy production term of q t is more sensitive to the temperature-humidity covariance θ ′ l q ′ t . As shown in Fig. 7e, the HET surface covariance reduces the buoyancy q ′ t θ ′ v . As a result, the moisture flux, w ′ q ′ t , in the lower atmosphere is decreased (Fig. 6f), which partially offsets the upward buoyancy flux enhanced by the buoyancy term of θ ′ l θ ′ v , as well as the turbulent kinetic energy. Accordingly, there is no significant difference in surface sensible and latent heat fluxes between HOM and HET simulations (not shown). The impact of the HET lower boundary condition does not affect the meteorological states near the surface (e.g., 2-m air temperature, 2-m specific humidity) in this SCM study. Figure 7g shows that the third-order moment of the subgrid vertical velocity w ′3 is very sensitive to the surface heterogeneity.
The w ′3 difference between HOM and HET extends to the entire ABL for the clear-sky days. The budget of w ′3 shows that the source due to buoyancy aloft is increased with heterogeneity in HET simulation, whereas it is partially offset by the decrease 195 in pressure (Fig. 8). In addition, the turbulent advection, which represents vertical transport by updrafts and downdrafts, is also responsible for increasing the HET w ′3 in the upper ABL. This change of w ′3 indicates that the subgrid vertical velocity characteristics have been altered when accounting for surface heterogeneity, as the skewness of the w ′ PDF, Sk w ≡ w ′3 /w ′2 3/2 , increases throughout the ABL (Fig. 7h). The Sk w changes can have implications for cloudy days, which will be discussed in the next section.

200
In summary, with the HET approach applied in E3SM, the diverging patch states play a dominant role in representing the spatial land-surface heterogeneity. The HET approach increases the variances of potential temperature and specific water content, as well as their covariance in the night time and early morning.

Impacts of surface heterogeneity on low clouds
In this section, we discuss the impacts of surface heterogeneity on low clouds in non-precipitating cloudy days. A total of 16 205 days were selected for this analysis when the maximum daytime (06-18 LT) liquid water content below 600 hPa exceeded 10 −3 g kg −1 and the daily precipitation rate (convective + large scale) was less than 0.5 mm day −1 . The daytime and nighttime (18-06 LT) averaged cloud properties and cloud liquid water flux profiles are depicted in Fig. 9. We find that the cloud fraction and cloud liquid water content show consistent vertical profiles (Figs. 9a-b). The HET simulation consistently produces more clouds than the HOM simulation. However, the turbulent liquid water flux are insensitive to surface heterogeneity (Fig. 9c).  Table 3. Daytime (06-18 LT) mean, relative difference, and Student t-test p-value of liquid water path, ice water path, low cloud fraction, total cloud fraction, net radiation flux, sensible heat flux, and latent heat flux, over the cloudy days in HOM and HET simulations.  Table 3 shows that the HET approach produces a 12.6% increase in low cloud fraction and a 8.4% increase in liquid water path, compared to the HOM simulation. These differences are found statistically significant, with p-values of 0.030 and 0.027, respectively. The difference in total cloud fraction and ice water path, however, does not pass the Student's t-test for statistical significance with a significance level of 0.05. This suggests that accounting for surface heterogeneity in landatmosphere coupling only affects low clouds in our E3SMv1 SCM configuration. As a result of more low-level clouds in the 215 HET simulation, surface radiative flux, as well as sensible and latent heat fluxes are reduced ( Table 3).

HOM mean HET mean Relative difference (%) T-test p-value
As the shallow convection mainly occurs in the afternoon over SGP, the vertical profiles averaged over 14-16 LT in the cloudy days are depicted in Fig. 10 with respect to the mean states and variances of θ l and q t as well as their covariance.
The effects of the HET approach in the subcloud layer (below 850 hPa) resemble those in the clear-sky days (Fig. 6). The differences stemmed from the surface extend vertically to the lower ABL. In cloud layers (between 650 hPa and 850 hPa; see Figure 10. Same as Fig. 6, except for non-precipitating cloudy days. Note the y-axis here spans higher than Fig. 6. Figure 11. Same as Fig. 7, except for non-precipitating cloudy days. Fig. 9a), the differences in the scalar variances are attributed to the cloud responses to the surface boundary conditions. Similar to the clear-sky days in Section 3.2, the turbulent heat flux is increased while the moisture flux is decreased comparably in response to the heterogeneous surface moments, thereby leading to slightly larger buoyancy flux and turbulent kinetic energy at the lower subcloud layer. In cloud layers, the differences in turbulent fluxes between HOM and HET configurations are negligible (Fig. 11).
225 Figure 11h shows that Sk w increases slightly below clouds and decreases in cloudy layers when accounting for surface heterogeneity. Since low Sk w indicates symmetric turbulent mixing corresponding to the stratocumulus regime and high Sk w indicates the presence of strong updrafts which corresponds to the shallow cumulus regime (e.g., Golaz et al., 2002Golaz et al., , 2005 13 https://doi.org/10.5194/gmd-2021-421 Preprint. Discussion started: 8 February 2022 c Author(s) 2022. CC BY 4.0 License. Cheng and Xu, 2006;Guo et al., 2014), the differences in Sk w between HOM and HET indicates a change in the turbulence structure and, consequently, the cloud regime. The decrease in Sk w within the cloud layers (Fig. 10d) indicates a shift toward 230 a more symmetric mixing within clouds, which leads to higher cloud fraction and liquid cloud content, consistent with the findings summarized in Table 3.

Impacts of surface heterogeneity on precipitation
In this section, we discuss the impacts on precipitation by accounting for surface heterogeneity in land-atmosphere coupling.
A total of 50 days in the simulation period were selected for this analysis when daily precipitation rate (convective + large 235 scale) exceeded 0.5 mm day −1 in both HOM and HET configurations. The results show that the mean precipitation rate in HOM is 6.1 mm day −1 , slightly lower than the precipitation rate of 6.4 mm day −1 in HET. This difference, however, does not pass the Student's t-test for statistical significance. The diurnal cycle of total precipitation is shown in Fig. 12a. Both HOM and HET simulations produce pronounced diurnal cycle of precipitation, with high precipitation rates in the afternoon mainly contributed by the convective precipitation (Fig. 12b). In the default HOM simulation, the precipitation rate peaks between 240 12 LT and 14 LT. In the HET simulation, the peak time is slightly delayed by about 1 h, but the difference does not pass the statistical significance test. We also find that the differences in precipitation PDF between HOM and HET are very small (< 1%; see Fig. 13).
The fact that the new HET approach does not influence the precipitation characteristics in the E3SM SCM is an interesting feature. We find that it is due to the fact that the inter-patch variability, which is the dominant term in the grid-aggregated variability representing the surface heterogeneity (Eq. (8)), vanishes when precipitation occurs. The sums of the area-weighted values are found to be close to zero in mid-June (see Fig. 5), coinciding with the occurrence of precipitation (Fig. 14). The same relationship between inter-patch variances and rainfall events is found in other time periods as well (not shown). This might be related to the fact that only the grid-box averaged precipitation flux is provided by the atmosphere model to the land model so that every land surface patch receives equal amount of water from the atmosphere model, which reduces the inter-250 patch variability. On the other hand, the changes in cloud cover and solar radiation simulated by the atmosphere model are also uniformly applied to the patches, so the patch differences may be especially small under this energy limited (precipitation) regime. Nonetheless, the mechanism driving this model behavior requires further investigation.

Conclusions
In this study, we implemented the Machulskaya and Mironov (2018) approach that aggregates the patch characteristics into 255 equivalent grid-cell quantities to account for surface heterogeneity in land-atmosphere coupling in E3SM. This new treatment adjusts the surface variances and covariance of potential temperature and specific water content for grid cells that enables a consideration of the spatial subgrid-scale heterogeneity in the land-atmosphere coupling. The grid-cell variances and covariance were computed offline and provided to EAM as part of the lower boundary conditions. We performed E3SM SCM hindcast simulations and compared the results with simulations using the default land-atmosphere coupling that neglects the impacts of    Using the new patch-aggregating approach with patch states and fluxes provided by ELM, we find that the inter-patch variability in potential temperature and specific humidity is the dominant term in representing surface heterogeneity in the grid-aggregated surface boundary condition. In clear-sky days, we find that the effects of surface heterogeneity extends from 265 the surface to the lower ABL. While accounting for surface heterogeneity increases temperature and humidity variances as expected, it reduces the temperature-humidity covariance in the afternoon due to the fact that drier/wetter patches have less/more evaporation and hence higher/lower temperature.
Furthermore, we find that the increase in the surface scalar variances as well as the decrease in the temperature-humidity covariance can directly influence the buoyancy terms and further the turbulent fluxes of temperature and moisture in the lower 270 ABL. We also find that the third-order moment of vertical velocity w ′3 is very sensitive to the new land-atmosphere coupling.
The increase in w ′3 induced by the heterogeneous surface properties is pronounced throughout the entire ABL.
In non-precipitating cloudy days, we find that low-level cloud fraction and liquid water content both increase when surface heterogeneity is considered. These changes are accompanied with the decrease in the skewness of the subgrid vertical velocity PDF, Sk w , indicating the change in the turbulence structure toward a more symmetric turbulent mixing in cloud layers. The 275 changes in low-level clouds also affect the surface radiative flux as well as sensible and latent heat fluxes.
Even though the new land-surface coupling introduces significant changes in clear-sky days and non-precipitating cloudy days, it does not have significant effects on precipitation. We find that the subgrid inter-patch variability vanishes when precipitation occurs, so that the surface moments computed using the new approach are similar to those computed by the original homogeneous treatment in land-surface coupling. As a result, precipitation is not affected by this new treatment.

280
The patch-aggregating approach considers the subgrid variability in surface states, improving the realism of the heterogeneous land-atmosphere interactions in E3SM. This study provides an initial assessment of this treatment in a state-of-the-art ESM in a SCM configuration. The significant effects of this treatment on turbulence structure in the ABL and on low-level clouds highlight the importance of accounting for the subgrid variability in land-atmosphere coupling. Future studies that assess the impact of the new treatment on the climate system in the fully coupled E3SM will be highly valuable.

285
Code and data availability.