How does cloud-radiative heating over the North Atlantic change with grid spacing, convective parameterization, and microphysics scheme in ICON version 2.1.00?

Cloud-radiative heating (CRH) within the atmosphere and its changes with warming affect the large-scale atmospheric winds in a myriad of ways, such that reliable predictions and projections of circulation require reliable calculations of CRH. In order to assess the sensitivities of upper-tropospheric midlatitude CRH to model settings, we perform a series of simulations with the ICOsahedral Nonhydrostatic Model (ICON) over the North Atlantic using six different grid spacings, parameterized and explicit convection, and one- versus two-moment cloud microphysics. While sensitivity to grid spacing is limited, CRH profiles change dramatically with microphysics and convection schemes. These dependencies are interpreted via decomposition into cloud classes and examination of cloud properties and cloud-controlling factors within these different classes. We trace the model dependencies back to differences in the mass mixing ratios and number concentrations of cloud ice and snow, as well as vertical velocities. Which frozen species are radiatively active and the broadening of the vertical velocity distribution with explicit convection turn out to be crucial factors in altering the modeled CRH profiles.

1 Introduction

Clouds have important radiative effects within the atmosphere. They absorb the outgoing infrared radiation that would otherwise escape to space and re-emit it at colder temperatures. They also absorb and reflect incoming solar radiation that would otherwise warm the atmosphere and surface. The relative balance of these warming and cooling effects depends on the cloud phase and altitude. The cooling effect tends to dominate for low-level liquid clouds, whereas the warming effect tends to dominate for high-level ice clouds.

Within the atmosphere, the impact of clouds on atmospheric radiation is generally quantified with cloud-radiative-heating rates, as this heating is what influences circulation. This cloud-radiative heating can be calculated as the difference between all-sky and clear-sky flux divergences. A local heating or cooling rate due to clouds translates to changes in atmospheric temperature and pressure gradients and, hence, the driving forces for winds. The notion that clouds are not only embedded in the circulation but also determine it has become an important theme in recent years within cloud and climate research (e.g., Bony et al., 2015; Voigt and Shaw, 2015; Voigt et al., 2020).

A burgeoning body of work highlights the many ways in which clouds affect circulation via their radiative heating. In the tropics, cloud-radiation interactions cause tightening of the ascent region and expansion of the descent region within the Hadley cell (Albern et al., 2018). Radiative heating from tropical upper-tropospheric clouds also contributes importantly to the eastward extension and strengthening of the North Atlantic jet stream over Europe under global warming (Albern et al., 2019, 2021). Radiative effects of tropical clouds push the midlatitude eddy-driven jet equatorward, while those of extratropical clouds push it poleward (Watt-Meyer and Frierson, 2017). A shift from upper-tropospheric cloud-radiative heating in the tropics to cooling in the midlatitudes also strengthens the meridional temperature gradient and, hence, baroclinicity and static stability (Li et al., 2015; Voigt et al., 2020). Biases in the Southern Hemisphere jet location have also been traced back to underestimated shortwave reflection by clouds there (Ceppi et al., 2012). With regard to internal variability, anomalies in cloud-radiative effects can prolong the North Atlantic Oscillation and intensify or mute the amplitude of the El Niño Southern Oscillation depending on the model framework (Papavasileiou et al., 2020; Rädel et al., 2016; Middlemas et al., 2017). A more exhaustive description of these multifaceted cloud-radiative-circulation couplings is provided by Voigt et al. (2020).

Constraining the cloud-radiative-heating (CRH) profile is essential then to understand current-day circulation, as well as its future changes with increased concentrations of atmospheric greenhouse gases. The vertical distribution of CRH, however, varies dramatically from one model to another and between models and satellite products (Cesana et al., 2019; Voigt et al., 2019). This variability is especially pronounced in the upper troposphere where ice clouds exist and is present even between different reanalysis datasets (Tegtmeier et al., 2022). Our previous work has explored this variability in tropical upper-tropospheric CRH (Sullivan and Voigt, 2021; Sullivan et al., 2022). Structural differences in ice microphysics, such as consistency (or lack thereof) in the treatment of ice crystal size or the initial size at which crystals are nucleated, are important drivers of CRH variability in storm-resolving simulations. High-resolution simulations also indicate that cloud macroproperties like the degree of vertical overlap or the decorrelation length between overlying cloud layers strongly influence radiative properties (Wang et al., 2021).

Wang et al. (2021) targeted tropical and Arctic mixed-phase clouds, and Sullivan and Voigt (2021) and Sullivan et al. (2022) focused on tropical ice clouds because of the large intermodel CRH variability in these regions. Wang et al. (2021) note the influence of the width of the hydrometeor size distribution on CRH errors, while Sullivan and Voigt (2021) pinpoint several ice microphysical factors, such as initial ice crystal size and autoconversion rates, that drive CRH variability. Cesana et al. (2019) have compared heating-rate profiles from several global climate models to CloudSat-CALIPSO data, and Hang et al. (2019) have produced a global climatology of radiative heating decomposed into cloud types from the CloudSat multi-sensor data. But sensitivities of midlatitude, atmospheric CRH to model settings remain relatively unexplored. (Senf et al., 2020) found strong grid-spacing dependence in shortwave top-of-atmosphere fluxes and a reduction in compensating longwave and shortwave biases at the finest grid spacings (∼ 2.5 km) over the North Atlantic. We extend their work on top-of-atmosphere fluxes to examine the in-atmosphere cloud-radiative heating here.

We also build upon recent interest in the grid spacing and microphysics dependence of cloud-radiative heating, looking at how these model settings affect heating rates over the North Atlantic (e.g., Gettelman and Sherwood, 2016; Evans et al., 2017; Vannière et al., 2019; Sullivan et al., 2022). We start by establishing the climatological representativeness of our simulated cloud-radiative heating and present its dependencies on model settings, both in the net and decomposed into longwave and shortwave components. We examine whether these dependencies are due to different frequencies of specific cloud classes or whether the clouds in these classes have different properties. We then trace the changes in cloud class occurrence and condensate back to cloud-controlling factors. We close by identifying three model aspects at the root of the variability in North Atlantic cloud-radiative-heating rates.

Figure 1 The NAWDEX simulation domain covers the entirety of the North Atlantic, as well as the northeastern Canadian seaboard, Greenland, Northern Africa, and Europe. The domain runs from 23 to 80^∘ N latitude and 78^∘ W to 40^∘ E longitude.

2 Methods

2.1 ICON Simulations

Simulations were performed with the ICOsahedral Non-hydrostatic model (ICON) version 2.1.00 of the German Weather Service and Max-Planck Institute for Meteorology over a North Atlantic domain between 23 and 80^∘ N latitudinally and between 78^∘ W and 40^∘ E longitudinally (Fig. 1). We use the same set of simulations as presented in Senf et al. (2020). A brief description of these runs is presented here. After removing the spinup period, the ICON simulations extend over 14 d during the North Atlantic Waveguide and Downstream Impact Experiment (NAWDEX) field campaign, specifically 21–25 and 30 September 2016, 1–5 October 2016, and 14–16 October 2016. NAWDEX was an international multi-aircraft field campaign that took place from 17 September to 22 October 2016 and that was based out of Iceland (Schäfler et al., 2018). NAWDEX studied midlatitude circulations, particularly warm conveyor belts, Rossby waves, and the North Atlantic jet stream, and the physical processes initiating and controlling them.

ICON was run during the NAWDEX period in numerical weather prediction (NWP) mode with the convection scheme of Tiedtke (1989), updated by Bechtold et al. (2008), used at all grid spacings. For the simulations at 2.5 km grid spacing, the deep-convection scheme or both the deep- and shallow-convection schemes are switched off in order to investigate the effect of explicit treatment of convection. The impact of cloud microphysics is explored by switching between the one-moment microphysics of Doms et al. (2005) used in the operational NWP mode and the more sophisticated and computationally expensive two-moment microphysics of Seifert and Beheng (2006), where heterogeneous nucleation is prescribed as in Hande et al. (2015). Although the two-moment microphysics scheme was developed for convection-permitting resolutions, we use it here in combination with parameterized convection also. For either the one- or two-moment scheme, the effective radius of cloud droplets or ice crystals is prescribed from the cloud liquid or ice water content, respectively; this formulation makes microphysics and radiation inconsistent in the two-moment case (Kretzschmar et al., 2020). ICON uses the generalized cloud overlap scheme of Hogan and Illingworth (2000) and a diagnostic cloud cover scheme based upon a probability distribution of vapor mass mixing ratios relative to saturation (Giorgetta et al., 2018). The rapid radiative transfer model (RRTM) evaluates fluxes in our simulations across 16 longwave and 14 shortwave spectral bands using a correlated-k method (Mlawer et al., 1997).

Finally, six different horizontal grid spacings are used to span the range from typical global climate model meshes down to storm-resolving ones: 80, 40, 20, 10, 5, and 2.5 km. Across these grid spacings, the number of grid cells varies by 3 orders of magnitude. In the discussion below, the simulation with a grid spacing of x km is sometimes referred to simply as the x km simulation. Vertical grid spacing is held constant at 75 levels. Lateral boundary conditions with 3-hourly frequency and initial conditions come from the Integrated Forecast System. Surface and aerosol data come from the German Weather Service. We filter out grid points corresponding to land and sea ice from the NAWDEX domain in our results below, focusing only on cloud fields over the ocean to remove differences due to surface albedo, surface temperature, or varying amounts of predicted sea ice.

2.2 Satellite, reanalysis, and AMIP-like data

We compare our heating-rate profiles to those from the 2B-FLXHR-LIDAR data, version P2R04 from CloudSat-CALIPSO data, binned to 2.5^∘ resolution (see Papavasileiou et al., 2020) and remapped to 0.25^∘ resolution, over the North Atlantic domain during September and October between 2006 and 2011. As for the NAWDEX simulation output, we mask the land and sea ice grid points. Ice and liquid effective radii and water contents measured by the CloudSat cloud profiling radar and temperature and humidity profiles from the European Center for Medium-Range Weather Forecast (ECMWF) have been fed to a two-stream radiative transfer model to compute 2B-FLXHR-LIDAR heating rates by L'Ecuyer et al. (2008). We also compare heating rates from the ERA5 reanalysis of the ECMWF to our ICON NAWDEX simulations (Hersbach et al., 2020). The ERA5 reanalysis assimilates radiances from both infrared sounders, such as AIRS and IASI, and geostationary satellites, such as GOES and Meteosat. Heating rates have then been generated within the reanalysis by applying RRTM and assumptions about ice crystal effective size and cloud condensation nuclei concentrations. We download these ERA5 heating rates at 0.25^∘ resolution over our domain from 2012 to 2016 in order to produce a climatologically representative profile.

We also present CRH profiles from other coarse-resolution, AMIP-like simulations with the ECHAM6 atmospheric component of the MPI-ESM model, the LMDz5A atmospheric component of the IPSL-CM5A model, and the ICON atmospheric model version 2.1.00 with a global R2B04 grid, corresponding to a horizontal grid spacing of approximately 160 km. These simulations employ climatological sea surface temperatures from the CMIP5 AMIP protocol and have been analyzed by Voigt et al. (2019). Their CRH profiles are evaluated from over 5 or more years so that we may interpret them as a North Atlantic climatology. In both the ICON NAWDEX and the AMIP-like simulations, cloud-radiative heating is calculated as the difference between all-sky and clear-sky flux divergences.

2.3 Cloud classes

Cloud layering strongly determines CRH, and decomposition of cloud fields into various cloud vertical structure (CVS) classes has proven to be useful in tracing the origins of atmospheric radiative warming and cooling (Oreopoulos et al., 2017; Lee et al., 2020). CVS classes build upon the International Satellite Cloud Climatology Project classification and are defined by cloud fraction thresholds at low (pressure (p) ≥680 hPa), middle (440 hPa $\le p\le \mathrm{680}$ hPa), and high (p≤440 hPa) altitudes. Oreopoulos et al. (2017) define a classification consisting of high, middle, low, high–middle, middle–low, high–middle–low, high–x–middle, high–low, middle–x–low, and high–x–middle–x–low clouds, as well as clear sky; altitude-1–altitude-2 denotes cloudiness at altitudinal range 1 separated by clear sky from cloudiness at altitudinal range 2, whereas altitude-1–x–altitude-2 denotes continuous cloudiness throughout altitudinal ranges 1 and 2.

Figure 2 The cloud vertical structure classification of Oreopoulos et al. (2017) employs cloud fraction in three altitudinal ranges – low, middle, and high – to define 11 classes. We use a subset of these, shown in the red box, and do not distinguish between continuous and discontinuous cloud layers. We also focus on upper-tropospheric CRH influenced mostly by a smaller subset shown in the blue box. Adapted from Fig. 1 of Oreopoulos et al. (2017).

Within the low–middle–high stratification, numerous possibilities exist when looking at the full cloud fraction field, as detailed in the Appendix of Oreopoulos et al. (2017). How many consecutive levels within an altitudinal range must have cloud fractions greater than the threshold for the whole range to qualify as cloudy? Or if 20 % of the cloud exists in the high altitudinal range and 80 % exists in the middle altitudinal range, should it then be classified as isolated middle or high–middle?

We are mostly concerned with a general sensitivity of CRH to isolated versus deeper clouds, so we bypass some of these subtleties by employing a simplified version of the CVS classification with eight classes: isolated high, isolated middle, isolated low, high–x–middle, middle–x–low, high–low, high–x–middle–x–low, and clear sky (Fig. 2). To categorize cloudiness in a given grid cell, thresholds in cloud fraction are verified for the low (p≤680 hPa), middle (440 hPa $\le p\le \mathrm{680}$ hPa), and high (p≤440 hPa) ranges. These two-dimensional low, middle, and high cloud fractions are calculated over the corresponding pressure ranges from the three-dimensional cloud fraction field using the generalized overlap assumption. If, for example, a column of grid cells has more than the threshold cloud fraction in all three ranges, it is classified as high–x–middle–x–low. Or if it has only more than the threshold cloud fraction in the low altitudinal range, it is classified as low. We do not make the distinction between continuous and discontinuous layers of cloudiness. Three sets of thresholds were initially used based upon the following percentiles in the cloud fraction distribution: 60th–60th–25th, 62nd–67th–30rd, and 65th–70th–35th for high, middle, and low altitudinal ranges or cloud classes (Table S1 in the Supplement). The cloud fractions associated with these percentile thresholds change by up to 1 order of magnitude; however, cloud fraction is generally larger than these threshold values when a cloud forms so that the occurrence probability of cloud classes is mostly insensitive to which thresholds are used (Fig. S1 in the Supplement). We show results from the intermediate set of thresholds.

2.4 Hackathon format

The results presented here were generated in a non-traditional hackathon format. Over the course of 2 years, our research group met intermittently for intensive, 3 d periods of data analysis and discussion. Three subgroups focused on the climatological analysis (Sect. 3.1), the cloud class decomposition (Sect. 3.2), and the cloud-controlling factors (Sect. 3.4.2). This format facilitated communication about Python tools to handle the large datasets and a unique, group approach towards performing and organizing analyses.

3 Results

3.1 Climatological cloud-radiative heating in the North Atlantic

Cloud-radiative-heating (CRH) profiles averaged over open ocean in the NAWDEX domain from three global climate model simulations provide a first estimate of variability in North Atlantic climatological CRH (Fig. 3). The most prominent intermodel differences are in the lowermost (p≥800 hPa) and uppermost (p≤300 hPa) troposphere. The atmospheric component of the IPSL-CM5A model predicts by far the largest cloud-radiative cooling in the boundary layer and upper troposphere (maxima of −2.2 and −1.1 K d⁻¹, respectively). These atmospheric coolings are more than 5 times the magnitude of those produced by the MPI-ESM model, while the CRH in ICON falls in between with larger boundary-layer cooling than MPI-ESM but smaller upper-tropospheric cooling. The altitudes of cloud-radiative-cooling maxima also vary by about 80 hPa between the models in both the lower and upper troposphere. CRH profiles averaged over all longitudes between 23 and 80^∘ N mirror those over the NAWDEX domain, meaning that this midlatitude variability is not concentrated only over the North Atlantic. We also note that, on the basis of the ICON simulations, September and October are representative months for the annually averaged North Atlantic CRH (ICON full year versus ICON September–October).

Figure 3 North Atlantic climatological cloud-radiative heating varies 5-fold in coarse-resolution global model simulations. Full (panel a) and upper-tropospheric (panel b) CRH profiles averaged over the NAWDEX domain (23 to 80^∘N and 78^∘ W to 40^∘ E) from the atmospheric components of the MPI-ESM, IPSL-CM5A, and ICON version 2.1.00 models, all with approximately 150 km horizontal grid spacing. The means between 23 and 80^∘N over all longitudes for the three models are shown in the dotted traces denoted as NH mid for Northern Hemisphere mid-latitudes. ICON profiles for both the full year and for only September and October (Sep+Oct in the dashed trace) are shown. The dashed black lines in panel (a) indicate the subset of pressures shown in panel (b).

The circulation effects of the differing CRH in these AMIP-like simulations have been discussed by Voigt et al. (2019); their +4 K simulations show that particularly large CRH differences with warming are concentrated in the upper troposphere. The increase of upper-tropospheric CRH with surface warming results in larger meridional temperature gradients and a poleward expansion of the Hadley cell and extratropical jets. Clear-sky radiative cooling by water vapor provides a strong constraint for upper-tropospheric cloud fraction and cloud top temperature globally (Thompson et al., 2017, 2019). This clear-sky constraint means that, if we can reliably construct the current-day upper-tropospheric CRH, we can also infer what its profile looks like under global warming. We emphasize that radiative cooling from extratropical low-level clouds has non-negligible effects on circulation, for example enhancing baroclinicity (Li et al., 2015). However, given the strong dependence of both current and future circulation on upper-tropospheric CRH, we choose to focus on the model dependencies above 5 km going forward.

We next examine the relative contribution of upper-tropospheric CRH to the total, time mean, spatial mean heating rate within our NAWDEX simulations (Fig. 4). This heating rate climatology for the North Atlantic is constructed from the simulations with the coarsest grid spacing (80 km) and includes the longwave and shortwave cloudy and clear-sky radiative heating rates, as well as dynamic, turbulent, convective, and microphysical heating rates:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\partial T}{\partial t}}& =({\displaystyle \frac{\partial T}{\partial t}}{)}_{\text{CRH}}+({\displaystyle \frac{\partial T}{\partial t}}{)}_{\text{Clr Sky}}+({\displaystyle \frac{\partial T}{\partial t}}{)}_{\text{Dyn}}\\ & +({\displaystyle \frac{\partial T}{\partial t}}{)}_{\text{Turb}}+({\displaystyle \frac{\partial T}{\partial t}}{)}_{\text{Conv}}+({\displaystyle \frac{\partial T}{\partial t}}{)}_{\text{Mphy}},\end{array}\end{array}$$

where the largest component comes from clear-sky longwave radiative cooling (LW Clr Sky) followed by the dynamic heating (Dyn) and clear-sky shortwave radiative heating (SW Clr Sky). Thereafter, from about 9 up to 11 km, the microphysical heating and longwave cloud-radiative cooling are largest, with the latter contributing 14 % to the overall budget. The three smallest components of the budget are convective heating, shortwave cloud-radiative heating, and turbulent heating at these altitudes.

Figure 4 The heating rate budget is dominated by clear-sky radiation and dynamics, but longwave cloud-radiative heating contributes non-negligibly in the upper troposphere. Spatial mean and time mean vertical profiles of heating-rate components at 80 km grid spacing in the one- (panel (a)) and two-moment (panel (b)) microphysics schemes. LW CRH is longwave cloud-radiative heating, SW CRH is shortwave cloud-radiative heating, LW Clr Sky is longwave clear-sky heating, SW Clr Sky is shortwave clear-sky heating, Dyn is dynamics, Turb is turbulence, Conv is convection, and Mphy is latent heating from microphysics and saturation adjustment.

The hierarchy and values of the heating rates are independent of whether we use a one- or two-moment microphysics scheme (Fig. 4a versus b). The longwave cloud-radiative heating profiles do differ qualitatively, however, in whether they exhibit an inflection point. While the longwave cloud component changes from cooling to heating around 7 km in the one-moment setup, it is exclusively cooling at the upper altitudes in the two-moment setup. These heating rates indicate that cloud-radiative heating, especially its longwave component, is non-negligible in the North Atlantic upper troposphere.

Figure 5 Microphysics and convection dependency in the net CRH profile is much stronger than grid-spacing dependency. Upper-tropospheric, time mean, and area mean net cloud-radiative heating from the ICON NAWDEX simulations at grid spacings from 2.5 up to 80 km with a one- (panel (a)) and two-moment (panel (b)) microphysics scheme. The 2.5 km simulations either use only the shallow-convection parameterization (shallow on) or explicitly represent both shallow and deep convection (explicit). The standard deviation and standard error over daily means are depicted as light- and dark-red shades atop the 80 km profile. Profiles from the ERA5 reanalysis in September (dashed black) and October (dotted black), as well as the CloudSat-CALIPSO 2B-FLXHR-LIDAR product (solid black), are also included.

We first construct net CRH profiles from our NAWDEX simulations across six horizontal grid spacings, with shallow convective parameterization only and explicit convection in the 2.5 km simulation and using two different microphysics schemes (Fig. 5). Grid-spacing dependence is subtle. Simulations with coarser grid spacing exhibit larger-magnitude upper-tropospheric CRH, but profiles fall within 1 standard deviation of the 80 km profile over most of the upper troposphere. The CRH changes qualitatively with the microphysics scheme from an S shape in the one-moment scheme (as in the AMIP-like profiles of Fig. 3b) to a uniformly cooling profile in the two-moment scheme.

The most dramatic change occurs in turning off the deep convective parameterization in the two-moment microphysics simulations (Fig. 5b). Omitting the deep convective parameterization in the 2.5 km simulations shifts the upper-tropospheric cooling peak upward by 2 km and narrows its vertical depth relative to the other simulations. The explicit representation of convection also produces prominent heating below 9 km, not present in the other two-moment simulations. Although these results are for the full simulation length in Fig. 5, they are robust for shorter durations down to a single day (Fig. S2 in the Supplement).

Figure 6 Model setting dependency appears in both the longwave and shortwave components. Upper-tropospheric, time mean, area mean shortwave (a, c) and longwave (b, d) cloud-radiative heating with all model settings as in Fig. 5. One- (a, b) and two-moment (c, d) microphysics schemes are shown, as well as profiles from the ERA5 reanalysis in September (dashed black) and October (dotted black) and the CloudSat/CALIPSO 2B-FLXHR-LIDAR product (solid black). Note the different x-axis limits on the (a, c) versus (b, d).

Decomposing the net CRH into its longwave and shortwave components, we find that model dependencies are not isolated within a single component (Fig. 6). Both the longwave and shortwave CRH change more strongly with microphysics and convective scheme than with grid spacing. Interestingly, while the magnitude of longwave cooling increases at coarser grid spacing, that of shortwave heating decreases. Because longwave cooling is about twice as large as shortwave heating, it dominates the net CRH dependence. The larger spread on the longwave profiles also shows that this component drives more of the CRH variability across days.

Atop the simulated CRH values – both net and decomposed into their longwave and shortwave components – we overlay both ERA5 reanalysis values and a CloudSat-CALIPSO climatology over the NAWDEX domain during September and October. ERA5 assimilates observed radiances but still makes cloud microphysical assumptions within its radiative transfer calculations along the lines of a one-moment scheme in which only cloud liquid and ice mass mixing ratios are tracked (e.g., Tiedtke, 1993; Forbes and Tompkins, 2011). The CloudSat-CALIPSO product (2B-FLXHR-LIDAR) incorporates cloud microphysical measurements into its calculation (Sect. 2.2). The ERA5 and CloudSat-CALIPSO profiles differ strongly from one another and from the simulations. The ERA5 profile has a muted version of the S shape from the one-moment simulations, whereas the CloudSat-CALIPSO profile shows uniform upper-tropospheric cooling by clouds as in the two-moment simulations.

Taking CloudSat-CALIPSO as our baseline, simulations with moderate grid spacing (10- or 20 km) and the two-moment microphysics compare most favorably. Using instead the ERA5 reanalysis as our baseline gives an indication of CRH with the cloud environment but not microphysics observationally constrained, and in this case, our simulations with the finest grid spacing (2.5 km) and two-moment microphysics compare most favorably. None of the one-moment profiles mirror the CloudSat-CALIPSO or ERA5 profiles especially well. The messy state of this evaluation highlights the following difficulty: cloud-radiative heating is not directly observed, even from satellites, and associated radiative transfer or microphysical assumptions complicate any model–measurement comparison.

3.2 Cloud class decomposition

We turn next to understanding the strong convective and microphysical scheme dependency in the upper-tropospheric CRH by breaking it down into that associated with various cloud classes. Such a decomposition allows us to determine whether CRH differences are due to variations in heating associated with a particular cloud class or variations in the probability of occurrence associated with a particular cloud class. Stated mathematically, the total CRH is the summation, over all clouds classes i, of the heating associated with a given cloud class weighted by its frequency of occurrence (f_i below):

$$\begin{array}{}\text{(2)}& \text{CRH}=\sum _i{\text{CRH}}_i{f}_i.\end{array}$$

As detailed in Sect. 2.3, eight cloud classes are defined on the basis of cloud cover in three altitudinal ranges. Upper-tropospheric CRH is driven primarily by four of these eight cloud classes: isolated high clouds, continuous high–x–middle clouds, layered high–low clouds, and deep high–x–middle–x–low clouds (blue box in Fig. 2). Physically, isolated high clouds correspond to either dissipating anvil outflow cirrus or cirrus formed in situ. High–x–middle–x–low clouds represent forms of midlatitude deep convection, such as cyclones. The profiles associated with the low, middle, middle–x–low, and clear-sky regions are generally omitted, as these contribute negligibly to the CRH between 5 and 15 km (not shown).

Figure 7 There are no systematic changes in cloud class occurrence with grid spacing. Area-weighted occurrence frequency for eight cloud classes across grid spacings for the simulations with two-moment microphysics. The box shows 25th (Q1), 50th (Q2), and 75th (Q3) percentiles. The whiskers show 1.5 times the interquartile range below the first quartile up to 1.5 times the interquartile range above the third quartile, i.e., [Q1 − 1.5(Q3 − Q1), Q3 + 1.5(Q3 − Q1)]. Diamonds indicate outliers. Figure S3 in the Supplement is the same plot for the one-moment microphysics. Thresholds of the 62nd, 67th, and 30rd percentiles of the cloud fraction distribution are used for high, middle, and low clouds, but mean occurrence is not sensitive to these thresholds (Fig. S1 in the Supplement). The 2.5 km simulation uses neither a deep nor shallow convective parameterization (explicit). The sum of occurrence over all classes equals 1, and the sum over all classes except clear sky equals mean cloud fraction.

Box plots of area-weighted occurrence frequency show negligible grid-spacing dependence for all cloud classes (Figs. 7 and S3 in the Supplement). For the classes including high clouds that are influential for upper-tropospheric CRH, the mean occurrence changes by less than 2 % between the simulations with 80 and 2.5 km grid spacings. Otherwise, these box plots indicate that low clouds are the most frequent, with a mean occurrence around 30 %, followed by deep clouds (H–x–M–x–L) and clear sky, both with mean occurrences of roughly 17 %. Isolated middle clouds are least common, followed by high–x–middle clouds, occurring at an average of 2 % and 3 % of the time, respectively. Isolated high clouds also occur less frequently in this region, with only 6 % coverage on average.

While the occurrence probabilities do not reflect the model dependencies of the net CRH, the cloud-class-filtered CRH does (Fig. 8). The isolated high clouds (high or high–low) uniformly radiatively heat the upper troposphere between 5 and 15 km, whereas deeper clouds (high–x–middle or high–x–middle–x–low) radiatively cool above about 8 km. Isolated high clouds absorb more outgoing longwave radiation (OLR) than clear sky, whereas deep clouds absorb this OLR in the liquid cloud at lower altitudes and re-emit it at colder temperatures from their cloud tops. For the isolated clouds, heating intensifies with finer grid spacing and especially with turning off of the convective parameterization with the two-moment scheme. In contrast, for the deeper clouds, cooling moderates with finer grid spacing. But again, the largest change in the radiative-heating profile comes from turning off the convective parameterization with the two-moment scheme.

Figure 8 Isolated high clouds heat and deep clouds cool the upper troposphere. All cloud classes containing high clouds contribute to the model dependencies in the CRH of the two-moment microphysics simulations. Upper-tropospheric, time mean, area mean net cloud-radiative heating for four of the eight cloud classes with all model settings as in Fig. 5.

Having looked at both f_i and CRH_i from Eq. (2), we conclude that the latter factor drives the overall CRH dependencies. In other words, different model settings do not change the distribution of occurrence of various cloud classes; they only change the CRH profiles associated with these cloud classes. Additionally, these changes are not limited to a single cloud class but rather appear across all of those containing high clouds.

3.3 Cloud properties by class

We have ruled out varying occurrences of different cloud classes and now turn to cloud properties – overall and within the cloud classes – as an explanation for the model dependencies of CRH. An increased magnitude of time mean and area mean cloud-radiative cooling or heating can be due either to a larger amount of condensate in the cloud, a greater coverage of the clouds, or both. We examine cloud liquid water (q_c), cloud fraction, and cloud ice mass mixing ratios (q_i) for the various simulation settings in Fig. 9. q_c increases slightly with finer grid spacing in the two-moment scheme; however, its values are insufficient to drive the model dependencies in CRH (Fig. 9d).

Figure 9 Changes in the cloud ice mass mixing ratio drive the model dependencies of upper-tropospheric CRH in the two-moment simulations. Upper-tropospheric, time mean, and area mean profiles of cloud water mass mixing ratio (panels a and d), cloud fraction (panels b and e), and cloud ice mass mixing ratio (panels c and f) for the one- (a–c) and two-moment (d–f) microphysics simulations, with all model settings as in Fig. 5.

Differences in cloud fraction qualitatively mirror those in CRH for the one-moment scheme (Fig. 9b); specifically, cloud fraction peaks at a lower altitude and has a larger maximum in the simulations with coarser grid spacing, as does the cooling in its net CRH profiles. The correspondence of cloud fraction and net CRH dependence is weaker in the two-moment simulations (Fig. 9e). Cloud fraction is about 2 % larger for the 2.5 km simulations, but otherwise there is no consistent trend with grid spacing or the altitude of maximum cloud fraction. This weak dependence of cloud fraction on model setting appears across the classes with high clouds (Fig. S4 in the Supplement).

The primary driving factor of the large CRH changes with two-moment microphysics and explicit convection is then q_i (Fig. 9f). The amount of cloud ice quadruples from about 5 mg kg⁻¹ in the 80 km simulation to about 19 mg kg⁻¹ in the two 2.5 km simulations (without shallow or any convective parameterization). The one-moment simulations show no such change in q_i with model settings (Fig. 9c). As in Sect. 3.2, we can decompose these q_i differences into those associated with various cloud classes. Figure 10 illustrates that the q_i increases with grid spacing are somewhat larger for the deeper cloud layers – the high–x–middle and high–x–middle–x–low classes – than for the isolated high clouds but occur qualitatively across all the classes with high clouds. Likewise, the lack of grid spacing and convection dependence in q_i for the one-moment schemes is uniform across classes; there are no compensating differences in q_i.

Figure 10 Cloud ice mass mixing ratio increases 4-fold from the coarsest to finest grid-spacing simulations. Diagnostic ice mass mixing ratios from one- (top panels) and two-moment (bottom panels) simulations for the four cloud classes that include high clouds with all model settings, as in Fig. 5.

The model uses only condensate mass to calculate CRH. However, CRH is also physically determined by hydrometeor number, and we examine cloud ice crystal numbers (N_i) from our simulations to understand how their omission may affect CRH. N_i profiles parallel q_i ones for the two-moment microphysics simulations (Fig. 11, top panels). The runs without a deep convective parameterization produce more than 4 times as many ice crystals as those with a convective parameterization. Not only is more ice mass produced in the clouds, it is also distributed over many more hydrometeors. In physically accurate frameworks, larger N_i should promote multiple scattering and eventual absorption of solar radiation, enhancing the shortwave heating peak (Fig. 6c). Distribution of ice mass over many more crystals could also prolong cloud lifetime and enhance CRH. Our simulations permit such a cloud lifetime effect insofar as it is independent of CRH, but the cloud occurrence and cloud fraction results above indicate that it is not dominant.

Along with liquid and ice crystals, upper-tropospheric clouds may also contain snow (q_s) and graupel (q_g). Whereas q_i showed no model dependency for the one-moment simulations, the maximum in q_s changes almost 2-fold from the 80 km simulation down to the 2.5 km one without convective parameterization (Fig. 11, bottom). This monotonic increase in q_s appears for all cloud classes with the largest-magnitude changes from deep clouds in the one-moment scheme. Similarly, the q_g maximum changes by 1 order of magnitude across these model settings between 5 and 15 km (Fig. S5 in the Supplement). It is important to note that snow and graupel do not interact with the radiative transfer scheme in ICON. This exclusion of certain hydrometeors from the radiation scheme is motivated in part by size and in part by lack of a corresponding fractional coverage variable (e.g., Xu and Randall, 1995). Graupel will tend to sediment out more rapidly than the time step used to call the radiation scheme, whereas the fractional coverage of snow, distinct from the liquid or ice cloud fraction, is not a tracked variable. We can therefore conclude that grid spacing dependence for the one-moment microphysics is concentrated in radiatively inactive cloud species.

Figure 11 Strong grid-spacing dependence appears in both the ice crystal numbers from the two-moment scheme and the snow mass mixing ratios from the one-moment scheme. Cloud ice crystal number from the two-moment simulations (top panels) and snow mass mixing ratios from the one-moment simulations (bottom panels) for the four cloud classes that include high clouds with all model settings, as in Fig. 5.

3.4 Understanding cloud property differences

As a final step, we ask why the high, high–x–middle, high–low, and high–x–middle–x–low clouds produce more ice and have slightly higher coverage in the two-moment simulations. We have advocated in our work on tropical cloud-radiative heating for process decomposition as a means of unraveling such differences (Sullivan and Voigt, 2021; Sullivan et al., 2022). This process decomposition can be done in a number of ways. Processes can be classified based upon the temperature range in which they are active to generate an altitudinally stratified recipe for CRH (Sullivan and Voigt, 2021). Processes can also be organized based upon when they occur within the cloud lifecycle to produce a temporally stratified recipe for CRH (Sullivan et al., 2022).

Here, processes are categorized as sources versus sinks of cloud ice. Then q_i variations are understood either in terms of differences in the source–sink formulations or in terms of differences in the inputs to these formulations:

$$\begin{array}{}\text{(3)}& (q_{\mathrm{i}},N_{\mathrm{i}})=\mathit{\varphi }(\text{CPs},\text{CCFs})-\mathit{\psi }(\text{CPs},\text{CCFs}),\end{array}$$

where ϕ and ψ represent microphysical sources and sinks, respectively; CP denotes a cloud parameter like the deposition density of ice crystals; and CCF denotes cloud-controlling factors, a term for the environmental conditions that determine cloud properties (e.g., Stevens and Brenguier, 2009).

Within the two ice microphysics schemes in ICON, ice mass can be consumed by autoconversion, melting, and sedimentation. Because q_c differences are so much smaller than those in q_i, we focus on sink processes that do not involve the liquid phase, namely autoconversion and sedimentation. Ice mass can also be generated by nucleation, droplet freezing, depositional growth, and riming. Somewhat larger cloud water mixing ratios at finer grid spacing in the two-moment simulations may contribute to slightly stronger riming and droplet-freezing tendencies (Fig. 9). However, these processes cannot be the primary driver for the q_i differences of much larger magnitude. We focus instead on nucleation and growth sources.

3.4.1 Cloud ice sources and sinks

Autoconversion is the process of converting between ice and snow, with its rate S_auc represented as follows in the two microphysics schemes:

$$\begin{array}{}\text{(4)}& {\displaystyle }{S}_{\text{auc, 1M}}=\left({\mathrm{10}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}\right)(q_{\mathrm{i}}-q_{\mathrm{i},\mathrm{0}}),\text{(5)}& {\displaystyle }{S}_{\text{auc, 2M}}=E_{\mathrm{ii}}{N}_{\mathrm{i}}{q}_{\mathrm{i}}G({\mathit{\delta }}_{\mathrm{i}},{\mathit{\theta }}_{\mathrm{i}}),\end{array}$$

where q_i,0 is a threshold ice mass mixing ratio before autoconversion initiates, set to 0 in the one-moment scheme; E_ii is the ice–ice collision efficiency; and G is a function of δ_i and θ_i, non-dimensional combinations of gamma distribution parameters representing the ice crystal sizes. The one-moment formulation simply transfers ice to snow over a fixed time constant. This sink is then much stronger than in the two-moment formulation, which incorporates dependence on the crystal numbers and relative sizes.

Snow and ice settle at the following terminal velocities in the one- and two-moment schemes:

$$\begin{array}{}\text{(6)}& {\displaystyle }{v}_{\mathrm{Ts},\mathrm{1}\mathrm{M}}=\left(\mathrm{7.37}\text{m}{\text{s}}^{-\mathrm{1}}{\text{kg}}^{-\mathrm{0.125}}\right)m_{\mathrm{s}}^{\mathrm{0.125}},\text{(7)}& {\displaystyle }{v}_{\mathrm{Ts},\mathrm{2}\mathrm{M}}=\left(\mathrm{8.156}\text{m}{\text{s}}^{-\mathrm{1}}{\text{kg}}^{-\mathrm{0.526}}\right)m_{\mathrm{s}}^{\mathrm{0.526}},\text{(8)}& {\displaystyle }{v}_{\mathrm{Ti},\mathrm{2}\mathrm{M}}=\left(\mathrm{317}\text{m}{\text{s}}^{-\mathrm{1}}{\text{kg}}^{-\mathrm{0.363}}\right)m_{\mathrm{i}}^{\mathrm{0.363}},\end{array}$$

where m_s is the snow crystal mass, and m_i is the ice crystal mass. Ice does not sediment in the one-moment scheme. For a range of hydrometeor masses ∼𝒪( 10⁻¹³ up to 10⁻¹⁰ kg), the terminal settling velocity for snow in the one-moment scheme is much stronger than that for either ice or snow in the two-moment scheme. The sedimentation sink is then also much stronger in the one-moment formulation.

Heterogeneous nucleation occurs on ice-nucleating particles (INPs), represented as follows in the one- and two-moment schemes, respectively:

$$\begin{array}{}\text{(9)}& {\displaystyle }{C}_{\text{INP, 1M}}=(\mathrm{1}\times {\mathrm{10}}^{\mathrm{2}})\text{exp}[-\mathrm{0.2}(T-\mathrm{273}\text{K}\left)\right],\text{(10)}& {\displaystyle }\begin{array}{rl}& C_{\text{INP, 2M}}=\\ & \left\{\begin{array}{l}(\mathrm{4.99}\times {\mathrm{10}}^{\mathrm{4}})\text{exp}[-\mathrm{0.2622}(T-\mathrm{237}\text{K}{)}^{\mathrm{1.2044}}]\\ (\mathrm{7.72}\times {\mathrm{10}}^{\mathrm{4}})\text{exp}[-\mathrm{0.0406}(T-\mathrm{220}\text{K}{)}^{\mathrm{1.4705}}\left]f\right({\text{RH}}_{\text{ice}}),\end{array}\right.\end{array}\end{array}$$

where T is subzero temperature, and RH_ice is the relative humidity with respect to ice. While the one-moment scheme represents only immersion nucleation (Eq. 9), the two-moment scheme represents both a relative humidity-dependent deposition nucleation and immersion nucleation (cases of Eq. 10). Both formulations predict exponential increases in INP as subzero temperature cools but with a much steeper slope in the two-moment than one-moment scheme. Conversely, the absolute INP number from the one-moment scheme is much higher (e.g., Sullivan et al., 2022, their Fig. 10a).

Finally, the rate of depositional growth S_dep is represented with a much more complicated temperature dependence in the two-moment scheme:

$$\begin{array}{}\text{(11)}& {\displaystyle }{S}_{\text{dep, 1M}}=(\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{5}})m_{\mathrm{i}}^{\mathrm{1}/\mathrm{3}}(q_v-q_{\text{sat,i}})\text{(12)}& {\displaystyle }{S}_{\text{dep, 2M}}={\displaystyle \frac{\mathrm{4}\mathit{\pi }{C}_i{D}_i{S}_{\mathrm{i}}f\left(m_{\mathrm{i}}\right)}{[\frac{RT}{p_{\text{sat,i}}\mathcal{D}}+\frac{L_{\mathrm{iv}}}{k_{\mathrm{i}}T}(\frac{L_{\mathrm{iv}}}{RT}-\mathrm{1}\left)\right]}},\end{array}$$

where q_v is the specific humidity, q_sat,i and p_sat,i are the saturation specific humidity and vapor pressure with respect to ice, C_i is the crystal capacitance, f(m_i) represents a mass-dependent ventilation coefficient, k_i is the thermal conductivity of ice, L_iv is the latent heat of sublimation, D_i is the diffusivity of vapor water, S_i is the saturation with respect to ice, and R is the gas constant. Key to both the nucleation and growth sources is the initial mass at which ice crystals are formed. The two-moment scheme initiates its crystals at 10⁻¹⁴ kg, and the one-moment scheme initiates its crystals at a much larger mass of 10⁻¹² kg (e.g., Sullivan et al., 2022, their Table 2). While the two-moment scheme generates fewer smaller crystals, they also stay aloft longer.

3.4.2 Cloud-controlling factors by class

Looking at the cloud ice source and sink formulations above, temperature (T), specific humidity (q_v), and vertical velocity (w) are the most important cloud-controlling factors (CCFs). T and q_v appear explicitly in Eqs. (9)–(12), while the influence of w is felt indirectly by determining saturation with respect to ice (RH_ice or S_i in Eqs. 10 and 12). The strength of w relative to v_Ts also determines whether ice crystals sediment. We examine these inputs across cloud classes and model settings (Fig. 12). Specific humidity differences from the 80 km simulation are quite small (Fig. 12, top row). The simulations with finer grid spacing are drier than the 80 km one below 10 km, but there is not a smooth trend toward lower specific humidity with finer grid spacing.

Figure 12 Simulations without convective parameterization exhibit stronger mean vertical velocities. Differences in the upper-tropospheric, time mean, and area mean specific humidity (top row) and temperature profiles (middle row) from that of the 80 km simulation. Time mean and area mean vertical velocity profiles for all simulation settings (bottom row – note the different x-axis limits in the first and third panels versus the second and fourth panels). Variables associated with the four cloud classes that include high clouds are shown for the simulations with the two-moment scheme only, with all model settings as in Fig. 5.

Profiles of temperature difference from the 80 km simulation mostly indicate a consistent trend of upper-tropospheric temperatures cooling as grid spacing is refined, aside from the 2.5 km simulations (Fig. 12, middle row). Across all classes with high clouds, the 40 km simulation is about 0.5 K cooler than the 80 km one between 5 and 11 km; the 5 km simulation is as much as 1.8 K cooler at these altitudes. These shifts toward colder temperatures below 11 km can help explain the increasing q_i there at finer grid spacings. Colder temperatures will accelerate nucleation of new crystals and depositional growth of existing crystals at warmer subzero temperatures in the two-moment scheme. However, the trend does not hold for the 2.5 km simulations without convective parameterization. Variations in input temperature cannot explain the dramatic increase in q_i with explicit convection.

Vertical velocities increase systematically with refined grid spacing, especially for the deep cloud layers (Fig. 12, bottom row). Deep cloud layers – high–x–middle and high–x–middle–x–low classes – are characterized by ascent throughout, whereas the isolated cirrus – high or high–low classes – have ascending air only above 7 km with descent below. Strengthening ascent will promote nucleation and growth in the same manner as cooling temperature. For the high–x–middle clouds, vertical velocity increases by a factor of 1.8 – from 1.2 to 2.2 m s⁻¹ between 80 and 2.5 km grid spacings. For the high–x—middle–x–low clouds, vertical velocity increases by a factor of 1.4 – from 2.5 to 3.5 m s⁻¹ between 80 and 2.5 km grid spacings.

A subtlety of vertical velocity is that a few instances of strong ascent can drive the majority of ice nucleation (e.g., Donner et al., 2016; Sullivan et al., 2016; Shi and Liu, 2016). The extreme values are more influential than the means depicted in Fig. 12, so we also construct the probability distribution of vertical velocities at 500 hPa from the various simulations (Fig. 13). We note that the ICON model uses no representation of subgrid-scale variability in vertical velocities. The variance of these resolved vertical velocity distributions becomes larger for finer grid spacing and without convective parameterization for both the one- and two-moment microphysics schemes. This distribution broadening indicates that vertical velocities, not only in the mean but also in the extremes, intensify at finer grid spacings.

Figure 13 Instances of the strongest ascent and descent at 500 hPa both become more probable at finer grid spacings. Probability distributions of vertical velocity from the simulations with the one-moment (panel (a)) and two-moment microphysics (panel (b)) with all model settings as in Fig. 5. Because we seek to explain the grid-spacing dependence of q_i, these velocities are not averaged or interpolated to a uniform grid.

A final factor to consider is the separation of convective and grid-scale microphysics with parameterized convection. Within a convecting grid cell, when convection is parameterized, the more sophisticated formulations of Eqs. (4)–(12) are superseded by simpler formulations in the convective microphysics. In particular, liquid condensate is converted to ice using a linear interpolation of temperatures between 273 and 235 K. As a result, the stronger vertical velocities at higher grid spacings have a particularly strong effect in the absence of convective parameterization, as they influence ice formation and growth in all grid cells, not only the non-convecting ones.

This analysis of source and sink processes and the cloud-controlling factors driving them produces a balance in favor of larger ice production within the two-moment scheme, especially with explicit convection. The most important elements in this balance are (1) weaker autoconversion and sedimentation sinks, (2) smaller initial crystal sizes, and (3) more instances of strong vertical velocity in the two-moment setup with explicit convection.

4 Conclusions

Given the importance of cloud-radiative heating – especially its upper-tropospheric values – to large-scale circulation features from the Hadley circulation to the eddy-driven jet, we have explored its dependencies on grid spacing, convective parameterization, and microphysics scheme in a numerical weather prediction model. The combination of parameterized versus explicit representation of convection and a one- versus two-moment microphysics scheme is the most influential model setting for CRH in our simulations. When we use a two-moment microphysics scheme, switching from parameterized to explicit convection has a much more dramatic effect than in the one-moment simulations. We posit that, when convection is parameterized, separation of convective and grid-scale microphysics produces a larger difference in the two-moment case. Sensitivities to grid spacing are more muted than those to the microphysics or convection parameterizations (Fig. 14). This result reflects the increased importance of constraining microphysical uncertainties as we transition toward the higher grid spacings of storm-resolving models.

Figure 14 A schematic overview of our analysis. Results of Sect. 3.1, 3.2, 3.3, and 3.4 are shown in pink, yellow, green, and blue, respectively.

Strong microphysical and convective sensitivity and weaker grid-spacing sensitivity in the CRH profiles do not appear in distributions of cloud class occurrence and appear only weakly in cloud fraction profiles. Instead, it is the cloud ice mass mixing ratio profiles that mirror the CRH dependencies most closely. We can trace these cloud ice mass mixing ratio differences back one additional step to changes in microphysical formulations and cloud-controlling factors (Fig. 14). Radiatively inactive frozen species, like snow and graupel, and the initial ice crystal mass, via its effect on subsequent growth and sedimentation rates, are two influential aspects of the microphysical formulations. Within the cloud-controlling factors, the width of the vertical velocity distribution and the upper-tropospheric temperature vary systematically with model setting.

Importantly, these findings are robust in relation to several factors. The dependencies affect both shortwave and longwave components of the cloud-radiative heating and occur across isolated cirrus, layered cirrus–boundary-layer cumulus, and forms of deep convection (high, high–low, high–x–middle, and high–x–middle–x–low in our decomposition). They are also not dependent on the cloud fraction thresholds used to define these cloud classes (Fig. S1 in the Supplement) or on the simulation duration. The grid-spacing and scheme dependencies already emerge within a single-day simulation (Fig. S2 in the Supplement). The upper-tropospheric CRH variability motivating this work also appears not only across three coarse-resolution global climate models (Fig. 3) but also across four reanalysis datasets (Tegtmeier et al., 2022) and between the ERA5 reanalysis and the CloudSat-CALIPSO 2B-FLXHR-LIDAR data (Fig. 5). While our analysis method could be generalized to other regions or modeling frameworks, the role of q_i and specific microphysical processes or parameters in CRH sensitivity will not necessarily generalize.

This last point highlights a challenge in further constraining atmospheric cloud-radiative heating: even our baseline contains uncertainties or assumptions. The disagreement between the ERA5 and CloudSat-CALIPSO profiles indicates that thermodynamic and wind fields are insufficient to constrain CRH. Both the one- and two-moment microphysics schemes generate quite similar distributions of cloud class occurrence despite drastically different upper-tropospheric CRH profiles (Figs. 7 and S3 in the Supplement). Stated another way, both cloud macrophysical and microphysical properties are needed to predict cloud-radiative heating. This result echoes our previous work on tropical CRH: 4-fold CRH variability can be produced by flipping ice microphysical switches in the model (Sullivan and Voigt, 2021), and cloud ice water content still changes 5-fold when inputs to microphysics schemes are fixed via a Lagrangian piggybacking technique (Sullivan et al., 2022).

A second challenge is that different combinations of model settings may improve model–measurement agreement in top-of-atmosphere or surface radiative fluxes versus atmospheric cloud-radiative heating. Models are often tuned based upon their outgoing longwave radiation, but as we have noted throughout, it is the in-atmosphere heating that feeds back upon circulation. As an example, Senf et al. (2020) assess the top-of-atmosphere cloud-radiative flux using the same set of runs and find the best agreement with the CloudSat-CALIPSO climatology from the simulation with 2.5 km grid spacing and shallow convective parameterization only. Here, it is instead the more moderate grid spacings and parameterized convection that agree best with the CloudSat-CALIPSO CRH (Fig. 5).

One suggestion for progress is to use atmospheric measurements to study atmospheric heating. This study and almost all other existing work use top-of-atmosphere (satellite) measurements to assess and investigate simulated atmospheric CRH. Although these satellite data have much better coverage and provide more robust statistics, it seems natural to use in situ radiative flux and cloud microphysical measurements to further investigate the in-atmosphere link of CRH to small-scale cloud properties. Our future work will adopt this approach in Lagrangian comparisons of simulated and in-situ-retrieved CRH over the Asian monsoon region.