Representing subgrid-scale cloud effects in a radiation parameterization using machine learning: MLe-radiation v1.0

Improvements of Machine Learning (ML)-based radiation emulators remain constrained by the underlying assumptions to represent horizontal and vertical subgrid-scale cloud distributions, which continue to introduce substantial uncertainties. In this study, we introduce a method to represent the impact of subgrid-scale clouds by applying ML to learn processes from high-resolution model output with a horizontal grid spacing of 5 km. In global storm resolving models, clouds begin to be explicitly resolved. Coarse-graining these high-resolution simulations to the resolution of coarser Earth System Models yields radiative heating rates that implicitly include subgrid-scale cloud effects, without assumptions about their horizontal or vertical distributions. We define the cloud radiative impact as the difference between all-sky and clear-sky radiative fluxes, and train the ML component solely on this cloud-induced contribution to heating rates. The clear-sky tendencies remain being computed with a conventional physics-based radiation scheme. This hybrid design enhances generalization, since the machine-learned part addresses only subgrid-scale cloud effects, while the clear-sky component remains responsive to changes in greenhouse gas or aerosol concentrations. Applied to coarse-grained data offline, the ML-enhanced radiation scheme reduces errors by a factor of 4–10 compared with a conventional coarse-scale radiation scheme. This shows the potential of representing subgrid-scale cloud effects in radiation schemes with ML for the next generation of Earth System Models.

1 Introduction

For climate projections, coarse-scale Earth System Models (ESMs) typically have horizontal resolutions of 100–200 km (Chen et al., 2021), in which clouds cannot be resolved explicitly. Therefore, these models require parameterizations of fractional cloudiness, particularly for cloud–radiation interactions at subgrid scales. A widely used approach is the Monte Carlo Independent Column Approach (McICA) (Pincus et al., 2003), where g-points are randomly assigned as cloudy or clear-sky according to the cloud fraction. This stochastic simplification introduces uncertainties in cloud–radiation interactions of up to 100 W m⁻² in surface fluxes, corresponding to relative errors of 10 % or more. However, these errors are unbiased compared to the Independent Column Approach (ICA) (Barker et al., 1999), where entire subcolumns are randomly designed as cloudy or clear, and thus tend to average out in long ESM simulations (Pincus et al., 2003). Vertical cloud overlap is commonly parameterized using the maximum-random overlap assumption (Räisänen et al., 2004), whereby adjacent layers overlap maximally while distant cloud layers overlap randomly. Alternatively, some models employ an all-or-nothing cloud cover scheme, which is a good approximation in high-resolution simulations (Giorgetta et al., 2022; Hohenegger et al., 2023).

Machine learning (ML)-based radiation emulators have been developed for more than two decades (Ukkonen, 2022; Yao et al., 2023; Roh and Song, 2020; Pal et al., 2019; Krasnopolsky et al., 2005; Lagerquist et al., 2023). One potentially appealing aspect of ML-based emulators is their relative speed compared to traditional radiative transfer models, which in principle allows more frequent calls during ESM integration. In practice, however, the speed-up potential has proven limited (Bertoli et al., 2025; Ukkonen and Chantry, 2025; Hafner et al., 2025a), with faster performance achievable through code optimization (Cotronei and Slawig, 2020; Ukkonen and Hogan, 2024). Moreover, improvements from more frequent radiation calls (Hafner et al., 2025a) remain marginal. Nonetheless, ML-based radiation emulators are valuable in specific context, such as differentiable ESMs (Kochkov et al., 2024; Klöwer et al., 2024) – where the derivative can be calculated for a complete model integration step – or GPU-based ESMs, where they are more energy efficient than physics-based schemes (Ukkonen and Chantry, 2025). Despite their high accuracy, emulators still face challenges in cloudy conditions (Hafner et al., 2025b), which may partly explain why they have not yet shown substantial improvements over traditional radiation schemes.

As more high-resolution Global Storm Resolving Model (GSRM) data become available, they offer increasing opportunities to enhance coarse-scale ESMs with machine learning. High-resolution simulations have been applied to nudge coarse-scale models toward fine-scale states (Bretherton et al., 2022), to learn subgrid tendencies directly (Heuer et al., 2024; Busecke et al., 2025), to infer subgrid effects from coarse-scale states (Grundner et al., 2022; Shamekh et al., 2023) and to learn all physics parameterizations (Watt-Meyer et al., 2024). Beyond these applications, high-resolution model simulations enable new strategies for representing fractional cloudiness and cloud overlap in coarse-scale ESMs using ML. Specifically, GSRM output can be coarse-grained to the resolution of the target ESM, and a neural network (NN) can then be trained to learn the subgrid-scale distribution of clouds from the underlying statistics. For instance, Henn et al. (2024) predicted coarse-grained cloud fields to reduce radiative biases; however, the cloud overlap assumption in the radiation scheme remains unchanged, limiting improvements. Leveraging ML to represent subgrid-scale cloud effects on radiation thus provide a promising pathway toward accurate radiation schemes in ESMs for climate projections.

In this work, we demonstrate how the representation of cloud radiative impacts on heating rates can be improved by separating them from the all-sky heating rates and learning the cloud contribution directly from high-resolution simulations. Because these simulations explicitly resolve cloud systems without relying on assumptions about their horizontal and vertical distributions, the ML model can implicitly account for subgrid-scale cloud effects on the radiative heating rates. A similar separation between clear-sky and cloudy radiative fluxes was previously proposed by Chevallier et al. (1998). Moreover, Meyer et al. (2022) focused on learning 3D cloud radiative effects. However, the explicit separation of cloud radiative impacts from all-sky radiation, especially including subgrid-scale effects, remains largely unexplored.

We aim to use ML to encode the true variability of vertical and horizontal subgrid-scale clouds and their radiative impacts from high-resolution simulations, rather than relying on statistical schemes. This raises the central question: “Can the representation of cloud–radiation interactions be improved by training on high-resolution simulation data?” To address this, we compare our ML approach against a state-of-the-art physics-based radiative transfer model, namely RTE+RRTMGP (Pincus et al., 2019) with McICA (Pincus et al., 2003) and maximum-random overlap (Räisänen et al., 2004) to represent subgrid-scale cloudiness. In particular, we evaluate how a physics-based coarse-scale radiation scheme performs when applied to coarse-grained data. This comparison allows us to disentangle differences arising from changes in the input parameters from those due to different resolutions.

This paper is structured as follows. Section 2 describes the method used to learn the cloud radiative impact from high-resolution simulation data. Section 3 introduces the main datasets, generated with the ICOsahedral Non-hydrostatic (ICON) model(Giorgetta et al., 2018, 2022) in both high- and coarse-resolution configurations, and includes a comparison of the input and output variables of the radiation parameterization across these resolutions. In Sect. 4, we present the main results by comparing the physics-based coarse-scale radiation parameterization with the ML-enhanced scheme on coarse-grained test data. Finally, Sect. 5 provides a discussion and conclusion.

Figure 1 Sketch of constructing the cloud radiative impact on heating rates. Radiation schemes calculate fluxes for the same scene once with and once without clouds resulting in all-sky and clear-sky fluxes. The corresponding heating rates can be inferred from the fluxes and the residual yields an approximation of the cloud radiative impact on heating rates for all layers in a column.

2 Learning the Cloud Radiative Impact on Heating Rates

We define the cloud radiative impact on heating rates in a column as the residual between the all-sky and clear-sky heating rates, where clear-sky represents the same atmospheric conditions as all-sky but without clouds (Fig. 1a). In ESMs such as the ICON model (Giorgetta et al., 2018, 2022), radiation parameterizations like RTE+RRTMGP (Pincus et al., 2019) first calculate gas optical properties and then add cloud optical properties. Then, these combined properties are used to calculate all-sky radiative fluxes, and clear-sky fluxes can be obtained omitting the cloud optical properties. Heating rates are then derived from the flux divergence for both all-sky and clear-sky conditions. The residual of these heating rates, which is the cloud radiative impact (CRI) on the heating rates, serves as the training target for the ML-based radiation parameterization:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial T_{\text{CRI}}}{\partial t}}={\displaystyle \frac{\partial T_{\text{all-sky}}}{\partial t}}-{\displaystyle \frac{\partial T_{\text{clear-sky}}}{\partial t}}.\end{array}$$

The heating rates $\frac{\partial T_k}{\partial t}$ in a layer k are calculated from the net flux F_Net at the layer boundaries $k\pm \frac{\mathrm{1}}{\mathrm{2}}$

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial T_k}{\partial t}}={\displaystyle \frac{F_{\text{Net},k+\mathrm{1}/\mathrm{2}}-F_{\text{Net},k-\mathrm{1}/\mathrm{2}}}{c_{\mathrm{v}}{m}_{\text{air}}}},\end{array}$$

where m_air is the mass of moist air per area, and specific heat at constant volume c_v scales with the tracer mixing ratios.

The central idea is to train a neural network (NN) that learns only the cloud impact on heating rates (see Fig. 1). Cloud–radiation interactions are subject to large uncertainties, since the subgrid-scale horizontal and vertical distributions of clouds are not resolved in coarse-scale ESMs. In the hybrid ML-physics radiation parameterization (Fig. 1b), the NN predicts only the cloud radiative impact, while the clear-sky component is still computed by the original radiation parameterization. This design ensures that the ML-enhanced radiation scheme retains sensitivity to changes in GHG and aerosols through the clear-sky part, potentially improving generalization across different climates. The ML-cloud component can respond to GHG and aerosols changes only indirectly, for example through modifications of the cloud distribution. However, this hybrid approach does not capture secondary effects arising from reflected radiation. The validity is further discussed in the Appendix D.

At first glance, a linear decomposition of the clear-sky heating rate and the cloud radiative impacts on heating rates may seem counterintuitive, given the inherently nonlinear interaction among specific humidity, clouds, and trace gases such as ozone (Bony et al., 2015). Nevertheless, we adopt a linear decomposition framework in this work, as illustrated on Fig. 1. Within this framework, the NN is tasked with learning the nonlinear relationships based on the prevailing atmospheric conditions and cloud-related variables (e.g., cloud liquid water and cloud ice).

2.1 Method

We use a bidirectional long short term memory (BiLSTM) based on Hafner et al. (2025b) to learn the cloud radiative impact in an atmospheric column. Bidirectional architectures are particularly well-suited for radiative transfer problems (Ukkonen, 2022; Yao et al., 2023; Ukkonen and Chantry, 2025; Bertoli et al., 2025). Unlike their common usage for temporal sequences, BiLSTMs for radiation scan the vertical dimension in both directions, resembling upward and downward fluxes. The NN consists of one BiLSTM layer with tanh activation and a hidden dimension of 96, with two LSTMs scanning the input vertical profiles in upward and downward directions. The combined output of the LSTMs is then processed by a dense layer that predicts the heating rate at each level. The training is split between shortwave (SW) and longwave (LW) radiation as shortwave temperature tendencies are only calculated for sun-lit areas and longwave temperature tendencies are always computed, totaling in 82k trainable parameters per NN.

The input variables are vertical profiles of mass mixing ratios of specific humidity q_v, cloud liquid q_l, cloud ice q_i, snow q_s and ozone O₃, plus the cloud fraction cl, air density ρ, and temperature T. For SW, we additionally use surface albedo α and incoming solar flux at the top of the atmosphere $F_{\downarrow ,\text{TOA},\text{SW}}$, which is the solar constant weighted by the solar zenith angle and change in Earth-Sun distance to account for daily and seasonal variations. For LW, we additionally use the surface temperature T_surf as input. As mentioned above, the output is the cloud radiative effect on heating rates derived from all-sky and clear-sky heating rates (Fig.  1).

Normalization is important for faster convergence of the training (LeCun et al., 2012) and generalization (Beucler et al., 2024). O₃, ρ, T, and T_surf are normalized using their mean values μ and standard deviation σ, also known as z-score normalization ($x_{\text{norm}}=\frac{x-\mathit{\mu }}{\mathit{\sigma }}$). The normalization factors are computed using all cells of one time step and are one-dimensional such that the vertical structure of the variables remain. The vertical structure is an important aspect that the BilSTM uses to make a vertically correlated prediction. The water related variables q_l, q_i, q_s are normalized by the ambient total (radiatively active) water ($q_{\mathrm{v}}+q_{\mathrm{l}}+q_{\mathrm{i}}+q_{\mathrm{s}}$). Here, radiatively active means used in the radiation scheme. q_v uses z-score normalization providing information that relates to the absolute mass mixing ratios. cl and α are not normalized as they already vary between 0 and 1. $F_{\downarrow ,\text{TOA},\text{SW}}$ is normalized by 1360 W m⁻², which is close to the solar constant. The cloud radiative impact on heating rates is only converted to K d⁻¹, which is on the order of one.

The loss we minimize during training consists of the sum of mean absolute error (MAE) and mean squared error (MSE). We use the Adam optimizer (Kingma and Ba, 2017), and set a learning rate of 10⁻³, which is reduced by a factor of 2 when the validation loss does not decrease by 0.1 % for 5 epochs. To avoid overfitting, we use early stopping, which stops the training if the validation loss does not decrease for 10 epochs.

3 Data

ICON is a weather and climate model permitting simulations across different resolutions, from a few to hundreds of kilometers. For global and long-term applications, the atmospheric component ICON-A (Giorgetta et al., 2018) often runs at 80 km horizontal resolution, and 47 vertical levels covering the altitude range 0–83 km, with parameterizations for radiation, cloud microphysics, turbulence, convection and gravity waves. The ICON model has the option to be run as a GSRM (Stevens et al., 2019; Giorgetta et al., 2022; Hohenegger et al., 2023). The high-resolution simulations used here follow the QUBICC (Quasi-Biennial oscillation in a changing climate) protocol from Giorgetta et al. (2022). The QUBICC simulations have 5 km horizontal resolution, and the vertical dimension spans 83 km on 191 levels. The high-resolution allows the model to run without a convection scheme and gravity wave parameterization, as these processes are starting to be resolved.

We performed QUBICC simulations that cover a total of 40 d evenly distributed across four months: November 2004, January, April, July 2005. The simulations have a physics time step of 40 s and a radiation time step of 6 min. These simulations are run with prescribed sea surface temperatures, sea ice concentrations, greenhouse gas concentrations but no aerosols. The outputs are saved every 192 min. The uneven output interval was chosen to cover a large variability of different solar zenith angles at different locations. The first 6 d of each 10 d period is used for training, the next 2 d for validation and the last 2 d for testing. Using QUBICC data to train our ML model for ICON-A requires coarse-graining the high-resolution QUBICC dataset. All variables are horizontally and vertically coarse-grained from high-resolution simulations as in Grundner et al. (2022). For horizontal remapping, the high-resolution cells are weighted by their cell area that is contained in each coarse cell. High-resolution cells can be contained fully or partially in a coarse cell, which is represented by a smaller cell area. Similarly, the layer thickness was used for vertical coarse-graining, corresponding to a weighted average. This method is valid for all input and output variables used here, i.e., concentrations and fluxes. Cloud-related variables are expressed as mass fraction, ensuring consistent vertical coarse-graining. Absolute mass variables, such as air mass m_air in Eq. (2), are coarse-grained by the (weighted) vertical sum that are partially or fully contained in a coarse layer. We discarded a small number of coarse-grained cells if their surface height showed inconsistencies, which can occur over complex terrain. Specifically, the coarse-grained surface height was computed from the high-resolution grid file and compared to the surface height from the coarse-scale grid file, which is slightly rotated and shifted. Consequently, some high-resolution cells are only partially contained in a coarse cell, which can lead to a small mismatch in surface height. Cells with deviations of more than 0.5 m in surface height were discarded. Then, we randomly sampled 35k and 5k grid points per time step for the training and validation set. For the test set, we randomly selected 75k cells per time step for LW and 35k cells for SW. This yields 2.3 million training samples, 260k validation samples, 1.9 million test samples for shortwave and 4.2 million test samples for longwave radiation.

For comparison with the high-resolution data, we also made coarse-scale ICON-A simulations for the same periods and with a similar configuration at 80 km horizontal resolution to compare differences in various distributions. The ICON-A simulations are based on the version 2.6.4 described in Giorgetta et al. (2018), while the QUBICC simulations are based on the version icon-2024.10 (ICON Partnership et al., 2024). To make the coarse-scale simulation more comparable, we used the same microphysics scheme (Doms et al., 2011) and no aerosols. See Appendix for a comparison with the default microphysics scheme. Both simulations use the radiation scheme RTE+RRTMGP (Pincus et al., 2019). ICON-A uses McICA and maximum-random overlap to represent subgrid-scale cloud-radiative impacts. In the QUBICC simulation, the all-or-nothing scheme is used, assuming horizontal homogeneous distribution of clouds. The physics and radiation time step is 6 min for the ICON-A simulation, matching the radiation time step in QUBICC.

Figure 2 Distributions of water related input variables for ICON-A and coarse-grained QUBICC data. The bold line shows the mean and the shaded area shows 95 % of the spread between the 2.5 % percentile and the 97.5 % percentile. The boxplot is limited by the minimum and maximum values. The box edges are defined at the 25 % and 75 % percentile of the distribution. The black line illustrates the mean of the distribution and the star is the median.

The remaining differences between the coarse-scale ICON-A and high-resolution QUBICC simulation include the horizontal and vertical resolution, the higher temporal resolution of physical process in QUBICC, which is required due to the high spatial resolution. Moreover, QUBICC intends to resolve gravity waves and (deep) convection. Additionally, the radiation scheme RTE+RRTMGP (Pincus et al., 2019) in QUBICC uses snow mixing ratio as an input. Other differences could be due to differences between the code versions and tuning, as we did not retune ICON-A.

3.1 Comparison of input and output variables

If an ML-based scheme is trained on high-resolution simulations like QUBICC and the goal is to couple it to a coarse-scale model like ICON-A, one needs to ensure that the distributions of input and output variables match. Otherwise, the ML-based scheme could be faced with out-of-distribution samples, which can lead to errors that quickly build-up while the model is integrated (Rasp, 2020). Therefore, we analyze systematic differences between the coarse-scale ICON-A simulations and the high-resolution QUBICC simulations. This analysis is conducted by comparing the spatial and temporal means and spread in the input and output used and produced by the radiation parameterization. Specifically, we focus on q_v, q_l, q_i, cl, total cloud cover, q_s, as well as longwave and shortwave heating rates. For the comparison, we use all samples in the test set. The samples of both simulations cover the same time period which is November 2004, January, April, and July 2005. Comparing two simulations with different grids enables us to uncover systematic differences, with a focus on identifying larger variations.

Distributions of water related input variables are shown in Fig. 2. The distributions of specific humidity look similar in the coarse-grained and coarse-scale data set, where the maximum differences between the means is 0.9 g kg⁻¹ (Fig. 2a). The spread in humidity values also overlaps for both simulations. Cloud water has higher values below 3 km for the coarse-scale simulation while cloud water is more evenly distributed throughout the troposphere for the QUBICC simulations (Fig. 2b). Here, the maximum difference of the mean values is 0.02 g kg⁻¹. The distributions of cloud ice have similar shapes, but the coarse-grained distribution peaks higher by 0.002 g kg⁻¹ (Fig. 2c). The mean cloud area fraction is on average larger by 3 % for the coarse-grained data set below 14 km (Fig. 2d). Snow peaks between 5–10 km (Fig. 2f). However, the vertical distribution is wider in the coarse-grained simulation. Nevertheless, snow was not used as input for radiation in ICON-A.

For comparability, the coarse-grained heating rates from QUBICC had to be rescaled by a factor $c_{\mathrm{v}}/c_p$ to account for differences in the two model versions used, where c_p is specific heat at constant pressure.

Figure 3 The distribution of shortwave (top row) and longwave (bottom row) heating rates in coarse-scale and scaled coarse-grained data. The bold line shows the mean and the shaded area shows 95 % of the spread, which is defined as the spread between the 2.5 % percentile and the 97.5 % percentile. The left column shows all-sky heating rates as it is used in the ICON model. The middle column shows clear-sky heating rate computed from clear-sky fluxes, which is a diagnostic output in the ICON model. The right column shows the cloud radiative impact on the heating rate which was computed by subtracting the clear-sky heating rate from the all-sky heating rate.

For SW all-sky heating rates, the mean profiles match very well and the mean difference is 0.18 K d⁻¹ (Fig. 3a). There are only small differences in the spread of heating rates in the troposphere. The heating rate is decomposed into a clear-sky heating rate – which is calculated from the clear-sky fluxes – and the cloud radiative impact on heating rates. Their distributions are also shown in Fig. 3b and c. The SW clear-sky heating rate has a mean difference of 0.11 K d⁻¹  for the mean profiles. For the SW cloud radiative impact, the mean difference is also 0.11 K d⁻¹ . Here, we expect small differences due to mostly resolved clouds in the coarse-grained dataset. However, there is no clear bias between coarse-scale and coarse-grained cloud impacts.

For the LW all-sky heating rates, the mean profiles match well and have a mean difference of 0.18 K d⁻¹ (Fig. 3d). However, the spread in heating rates is slightly different, which co-occurs with the different spread in cloud water at 1 km (Fig. 2b). The LW clear-sky heating rates look very similar in their mean values and their spread where the mean difference is 0.14 K d⁻¹ (Fig. 3e). The LW cloud impact is concentrated in the troposphere (Fig. 3f). The mean impact is almost the same between coarse-scale and coarse-grained simulations with a mean difference of 0.09 K d⁻¹. However, there is a difference in spread, which can be expected, because the coarse-grained data set implicitly includes subgrid-scale cloud effects.

In the unscaled comparison, the clear-sky heating rates show the same bias as the all-sky heating rates for both SW and LW (Fig. B1). However, this bias does not directly translate to the cloud radiative impact on heating rates because adding the cloud impact is a highly non-linear process (Bony et al., 2015). Additionally, this indicates that the mean cloud effect is similar for a (quasi)-hydrostatic and a non-hydrostatic model.

Table 1 All datasets used in this study. Online means that a simulation was run in high or coarse resolution and the output was saved for certain time steps. Offline means that a parameterization (RTE+RRTMGP or MLe-radiation) computed its tendencies based on model output.

Figure 4 Comparison of the baseline and the ML-enhanced radiation scheme on coarse-grained QUBICC data. Results are shown for the shortwave (top rows) and longwave (bottom rows) spectral range. The results are shown separately for clear-sky samples (no clouds, left column), fully cloudy sky samples (second column), and samples with partial cloudiness (third column). Additionally, the results are separated by non-precipitating (fourth column) and precipitating clouds (fifth column). The shown metrics are coefficient of determination R² (green), bias (orange) and MAE (blue) with 95 % of the spread, which is defined as the spread between the 2.5 % percentile and the 97.5 % percentiles. The bias and MAE share the x-axis. The ML-clear-sky panels are gray because the clear-sky fluxes are not calculated by the ML model, but shown as reference.

4 Results

As mentioned above, the simulation outputs of QUBICC and ICON-A can't be used for a sample-by-sample comparison. Therefore, a coarse-scale radiation reference is need that is computed (offline) based on coarse grained QUBICC output. All the dataset used in this study are summarized in Table 1.

Figure 5 As Fig. 4, but for the selected regions shown in Fig. 5 of Bock and Lauer (2024).

For the coarse-scale radiation reference, we use the Python front-end (hereafter pyRTE) (Pincus et al., 2025) of the radiation scheme RTE+RRTMGP (Pincus et al., 2019), which is used in all our simulations. Using pyRTE, we replicate the implementation of radiation in QUBICC as closely as possible. For subgrid-scale variability, we employ McICA (Pincus et al., 2003) together with maximum-random cloud overlap (Räisänen et al., 2004). The procedure is as follows: first we calculate gas optical properties, assign them to the atmospheric state and calculate clear-sky fluxes. Next, we compute cloud optical properties, apply McICA with maximum-random overlap to represent subgrid-scale variability, and add cloud optical properties to the gas optical properties. Because QUBICC also includes snow in its radiation parameterization, we additionally calculate snow optical properties and combine them with the other optical properties, before computing all-sky fluxes. Then, heating rates are obtained applying Eq. (2). We calculate heating rates for all samples in the test dataset. The output of pyRTE is used as baseline as it represents the coarse-scale radiation scheme. This allows a sample-by-sample comparison to the coarse-grained heating rates obtained from QUBICC (first and third row in Figs. 4 and 5). The ML-enhanced radiation scheme predicts the cloud radiative impact, which is added to clear sky heating. The resulting all-sky heating is compared to the coarse-grained QUBICC heating rates (second and third row in Figs. 4 and 5). As evaluation metric, we use mean absolute error (MAE), bias, and the coefficient of determination (R²). The improvement is measured by comparing MLe-radiation to the baseline. The results are presented in Fig. 4 with the column-averaged metrics summarized in Table A1. We restrict the presentation to the lowest 20 km of the atmosphere, as cloud impacts are most relevant in the troposphere (see Fig. 3). Nevertheless, the NN predicts the cloud impact on the entire atmospheric column, and the results for the full column can be found in the Appendix C1.

The first column of Fig. 4 shows clear-sky heating rates, for which the cloud distribution plays no role. Accordingly, only small and statistically insignificant differences are expected, arising mainly from variability in water vapor. While coarse-grained QUBICC heating rates reflect subgrid-scale variability of water vapor, coarse-scale radiation schemes assume a horizontal homogeneous distribution of water vapor. Therefore, the assumption of horizontally homogeneous input parameters introduces a small error of 0.367 K d⁻¹ (SW) and 0.571 K d⁻¹ (LW) for the baseline. For comparison with other studies, Table A1 also reports the root mean squared error (RMSE). For clear-sky heating rates of the coarse-scale baseline, the RMSE is 0.443 K d⁻¹ (SW) and 0.688 K d⁻¹ (LW) compared with the coarse-grained QUBICC heating rates. Hogan and Matricardi (2022) developed a fast tool for computing gas-optical properties and reports an RMSE of 0.18 K d⁻¹. Although the error source differs (gas optical properties vs. spatial resolution), the magnitudes are comparable.

For completeness, we also evaluated the ML-model on clear-sky samples. However, it is not intended for clear-sky scenes as the cloud impact is zero. Therefore, the corresponding panels are grayed out. The MAE is 0.049 K d⁻¹ for SW and 0.028 K d⁻¹ for LW.

The second column of Fig. 4 shows results for fully cloudy samples (total cloud cover of 100 %). For the baseline, the MAE peaks near 10 km, exceeding 0.5 K d⁻¹ for SW and 1 K d⁻¹ for LW. The corresponding R² are low, with averaged values of 0.83 (SW) and 0.66 (LW), compared to 0.98 for the ML-enhanced scheme. The average R² values are computed by averaging over the vertical levels. R² is weighted by variability, and an R² of zero indicates that the error is as large as the variability itself. Since, McICA within a coarse-scale radiation parameterization is supposed to produce unbiased noise, the R² is therefore a less informative metric. In contrast, the bias reveals that the coarse-scale radiation scheme systematically struggles to represent the cloud impact near 10 km, particularly for SW. The ML-enhanced scheme produces nearly unbiased heating rates, with MAEs of 0.106 K d⁻¹ (SW) and 0.127 K d⁻¹ (LW), representing errors 4–6 times smaller than those of the baseline.

The third column of Fig. 4 shows results for partially cloudy scenes, defined here as total cloud cover between 10 %–90 %. In these cases, both schemes exhibit smaller errors than in fully cloudy scenes, consistent with the weaker overall cloud radiative impact. For the baseline, the bias is substantially smaller than fully cloudy conditions but still show a pronounced peak near 1 km for LW and a double peak at 1 and 10 km for SW. In contrast, the ML-enhanced samples produces nearly unbiased heating rates, with an MAE of 0.082 K d⁻¹ (SW) and 0.068 K d⁻¹ (LW), representing errors that are 5–10 times smaller than those of the baseline.

To interpret the double-peak error observed in partially cloudy samples, we divided the samples into precipitating and non-precipitating clouds, as a rough proxy for deep and shallow convection. Non-precipitating (warm) clouds were identified using thresholds of maximum 0.01 mm h⁻¹, total cloud cover larger than 10 %, and vertically integrated ice water path of less than 1 g m⁻². For reference, drizzle has a precipitation rate of 0.2 mm d⁻¹ (Wood, 2012). Precipitating clouds were identified by total cloud cover > 10 % and precipitation more than 3 mm h⁻¹. For reference, Zhao et al. (2024) reports mean precipitation of 3.5 mm h⁻¹ for deep convective cores over the tropical ocean. For non-precipitating clouds, both the coarse-scale radiation scheme and the ML-enhanced scheme show an MAE peak near 1 km for SW and LW. However, the ML-enhanced scheme achieves substantially smaller errors of 0.080 K d⁻¹ for SW and 0.069 K d⁻¹ for LW, which are 6–10 times lower than those of the baseline. For the precipitating clouds, the baseline exhibits an MAE peak at 10–12 km, while the ML-enhanced scheme shows enhanced error in the upper troposphere but without distinct peak. Instead, the ML-enhanced scheme shows a broader peak in MAE between 12–14 km. On average, however, the ML-enhanced error remain about 4–5 times smaller than those of the baseline.

We further evaluated the performance of the ML-enhanced scheme across five selected regions characterized by different predominant cloud regimes, following the classification of Bock and Lauer (2024). The results are shown in Fig. 5 and summarized in Table A2. In the arctic region (70–90° N), errors remain confined below 10 km consistent with the lower tropopause height in this region. For the baseline, the SW MAE is 0.215 K d⁻¹, even smaller than in the clear-sky conditions, although the spread in MAE is slightly larger. For LW, the MAE is 0.614 K d⁻¹, exceeding the clear-sky values. In contrast, the ML-enhanced scheme achieves errors that are 4–8 times smaller than those of the baseline.

In the Southern Ocean (30–65° S), the baseline exhibits large errors of 0.417 K d⁻¹ for SW and 0.760 K d⁻¹ for LW, with a characteristic double peak in the MAE at 1–2 and 10 km. The ML-enhanced scheme also reproduces this double peak structure in the MAE spread but reduces the MAE by a factor of 4–7 relative to the baseline. Over the tropical ocean (30° N–30° S), the baseline shows large errors in the upper troposphere, likely associated with deep convection. In this region, the ML-enhanced scheme again reduces the MAE by a factor of 5–9. A subregion within the tropical region, the Pacific ITCZ region (0–12° N, 135° E–85° W), shows similar behavior but with even larger errors at higher altitudes.

The stratocumulus region is represented by three sub-regions: the Southeast Pacific (10–30° N, 75–97° W), Southeast Atlantic (10–30° S, 10° W–10° E), and Northeast Pacific (15–35° N, 120–140° W). In these regions, both models show two peaks in the MAE: one in the lower troposphere at 1–2 km and another in the upper troposphere at 10–12 km. Notably, the upper tropospheric peak is larger in both SW and LW for both models. Nevertheless, the ML-enhanced scheme achieves errors 5–9 times smaller than those of the baseline.

5 Conclusions

ESMs struggle to represent subgrid-scale cloudiness and commonly rely on statistical schemes, such as McICA, to account for subgrid-scale cloud radiative effects (Pincus et al., 2003). Although random, unbiased errors can be mitigated by large-scale atmospheric mixing (Pincus et al., 2003), the column-by-column error can be large. ML algorithms trained on high-resolution, global storm-resolving simulations now provide an opportunity to represent fractional cloudiness in radiative transfer more accurately. To bridge scales, the high-resolution model output is coarse-grained to the target resolution, such that the coarse-grained variables implicitly retain subgrid-scale effects. Then, the coarse-grained variables implicitly include subgrid-scale effects.

We developed a hybrid physics-ML radiation parameterization, where the physics-based component computes clear-sky fluxes, while the ML component predicts cloud impact, implicitly accounting for subgrid-scale variability. This ML-enhanced framework offers a more robust and generalizable radiation scheme: the physics-based parameterization retains its responsiveness to changes in GHGs and aerosols, thereby mitigating potential out-of-distribution issues in climate projections. The ML component is implemented as a BiLSTM neural network, which has previously demonstrated strong performance in radiation applications (Ukkonen, 2022; Yao et al., 2023; Ukkonen and Chantry, 2025; Bertoli et al., 2025; Hafner et al., 2025b). For training, we use data from high-resolution QUBICC simulations with a horizontal resolution of ≈ 5 km and 191 vertical layers expanding up to 83 km. These fields are coarse-grained to ≈ 80 km and 47 vertical layers, matching the target resolution for a coarse-scale ESM. For comparison and to assess systematic differences between high-resolution and coarse-scale models, we additionally perform a coarse-scale ICON-A simulation. The distributions of the relevant input variables are found to be comparable between the coarse-scale and coarse-grained simulations.

We find that a coarse-scale radiation scheme (baseline) performs well for clear-sky samples, but exhibits large errors in cloudy conditions, reflecting its inability to represent subgrid-scale distributions. In contrast, the ML-enhanced radiation scheme consistently outperforms the baseline, reducing errors from unresolved clouds in the radiative transfer calculations by a factor of 4–10. Although the ML-enhanced radiation scheme does not explicitly resolve subgrid-scale distributions, it learns how specific combinations of grid-scale mean states map to heating rates that implicitly include subgrid effects. THe baseline showed substantial biases of 1–3 K d⁻¹ in the upper troposphere at 10–15 km for precipitating clouds, highlighting the strong influence of subgrid-scale cloud ice on heating rates. In general, both coarse-scale radiation scheme and the ML-enhanced scheme produce larger errors in cloudier conditions, but the ML-enhanced scheme consistently yields smaller errors. These results emphasize the need for a more explicit treatment of subgrid-scale clouds, particularly in the upper troposphere.

Therefore, we conclude that high-resolution model data combined with ML can improve the representation of cloud–radiation interactions in coarse-scale radiation parameterizations. Nevertheless, the presented approach has caveats. High-resolution simulations at 5 km horizontal resolution cannot resolve shallow convection directly, leaving associated cloud radiative effects on heating rates – particularly within the planetary boundary layer – unresolved (Stevens et al., 2019). Using finer horizontal resolutions could help reduce these uncertainties. In addition, aerosols and heterogeneous GHG concentrations are typically absent from current high-resolution models. If future simulations include substantial variability in GHG and aerosol concentrations, these could be incorporated as additional NN inputs to capture secondary effects of reflected radiation.

One of the next steps is the online implementation of the ML-enhanced radiation scheme in a coarse-resolution model such as ICON-A. Although the online stability remains to be tested, the comparison with the coarse-scale model and results from previous stable hybrid simulations (Hafner et al., 2025a) are promising, suggesting potential improvement for climate projections. An additional advantage of the presented scheme is that clear-sky fluxes can be computed less frequently, while the cloud radiative impact can be updated every time step. This provides a pathway to both reducing computational costs and improving the representation of cloud–radiation interactions.

Appendix A: Default microphysics scheme

One major difference between the ICON-A and QUBICC versions is the microphysics scheme. ICON-A uses a modified version of the Lohmann and Roeckner (1996) scheme and the QUBICC simulation uses the graupel scheme described in Doms et al. (2011). While both schemes are single-moment schemes, the latter treats precipitating tracers like snow, rain and graupel as prognostic variables while the former only diagnoses snow and rain. Moreover, it is known that cloud ice is too large in the upper troposphere in ICON-A (Doktorowski et al., 2024), which we also see when comparing cloud ice in ICON-A and QUBICC (Fig. A1).

Figure A1 As Fig. 2 but for the default microphysics scheme based on Lohmann and Roeckner (1996).

Table A1 Bulk statistics for heating rate results of the coarse-scale ML-based radiation emulator on coarse-grained QUBICC data. MAE is mean absolute error and R² is coefficient of determination. RMSE is root mean squared error. The percentage values in parentheses denote the relative values of MAE, bias and RMSE.

Table A2 Bulk statistics for results of cloud radiative effect on heating rates. MAE is the mean absolute error and R² is the coefficient of determination. RMSE is the root mean squared error. The percentage values in parentheses denote the relative values of MAE, bias and RMSE.

Appendix B: Calculation of heating rates

When comparing the unscaled heating rates between coarse-scale and coarse-grained simulations, we find huge biases of up to 10 K d⁻¹ between the mean heating rates, especially in the stratosphere. We found that there is a difference in heating rate calculation between the code versions. Usually, the heating rate is calculated from flux divergence, pressure difference and constants:

$$\begin{array}{}\text{(B1)}& {\displaystyle \frac{\partial T_k}{\partial t}}={\displaystyle \frac{g}{c_p}}{\displaystyle \frac{F_{\text{Net},k+\mathrm{1}/\mathrm{2}}-F_{\text{Net},k-\mathrm{1}/\mathrm{2}}}{P_{k+\mathrm{1}/\mathrm{2}}-P_{k-\mathrm{1}/\mathrm{2}}}},\end{array}$$

where g is gravitational acceleration, c_p specific heat at constant pressure, F_Net is the difference between downward and upward flux, P is pressure. k is defined at the center of a layer (also full levels) while $k\pm \frac{\mathrm{1}}{\mathrm{2}}$ is defined at the layer boundaries (also half levels). This form of converting fluxes to heating rates is usually found in hydrostatic models but does not work in ICON because pressure is a diagnostic variable. Instead, the density is kept constant and specific heat at constant volume c_v needs to be used for the conversion (Zängl et al., 2014). This transforms Eq. (B1) to Eq. (2). In the ICON-A version, they use c_p which is valid for quasi-hydrostatic models because then the hydrostatic pressure holds $\mathrm{d}P/g=\mathit{\rho }\mathrm{d}z$. Additionally, the specific heat is scaled only by water vapor, while all tracers are included in the code version used for QUBICC. Therefore, we scaled the coarse-grained heating rates in Fig. 3 with the ratio $c_{\mathrm{v}}/c_p$, which is on average 0.7. The unscaled heating rates are shown in Fig. B1.

Figure B1 As Fig. 3 but here the coarse-grained heating rates are not scaled.

Appendix C: Results for full vertical column

Figure C1 As Fig. 4 but for the full column. The large errors in the upper stratosphere for the baseline are related to rounding errors increasing with height.

Appendix D: Linear decomposition under changing CO₂

The linear decomposition assumption of clear-sky heating and cloud radiative impact may be questionable for different greenhouse gas concentrations. To provide some validity, we estimate the error induced by this assumption for 4×CO₂. For a direct estimation, we select 4000 random samples from the test set and calculate the cloud radiative impact for the reference climate (CO₂ concentration of 2004) offline using pyRTE. Next, we increase CO₂ by a factor of 4 and repeat the calculation. The mean absolute difference between cloud_effect_reference and cloud_effect_4xCO2 gives an estimate of the error that is induced by the linear decomposition assumption for 4×CO₂. The vertical resolved error and bias are shown in Fig. D1. The x-range is the same as in Figs. 4, 5 and C1 for comparison with the errors induced by subgrid-scale clouds and the remaining errors of MLe-radiation. In general, the error is smaller in the stratosphere than in the troposphere. For SW, the error is around 2.5 % and for LW it is around 4 % in the troposphere. This around 10 times smaller than the error from subgrid-scale clouds and around the same magnitude as the errors from MLe-radiation.

Figure D1 Mean absolute error (MAE) and bias induced by the linear decomposition assumption for 4×CO₂. For comparison, the x-range is the same as in Figs. 4, 5 and C1.