2.1. Split conformal prediction without retraining (SCP)

More formally, given $ n $ training samples $ {\left\{\left({X}_i,{Y}_i\right)\right\}}_{i=1}^n $ , the main objective is to construct the prediction interval $ \mathcal{C}\left({X}_{n+1}\right)\subseteq \mathrm{\mathbb{R}} $ for a test point $ {X}_{n+1} $ . The interval $ \mathcal{C} $ contains the unknown value $ {Y}_{n+1} $ with the miscoverage rate $ \alpha \in \left(0,1\right) $ , satisfying the marginal coverage guarantee

(2.1) $$ \mathrm{\mathbb{P}}\left\{{Y}_{n+1}\in \mathcal{C}\left({X}_{n+1}\right)\right\}\ge 1-\alpha, $$

for any joint distribution of the feature vectors $ X\in {\mathrm{\mathbb{R}}}^p $ of dimension $ p\ge 1 $ and labels $ Y\in \mathrm{\mathbb{R}} $ and any sample size $ n $ with the assumption that all samples $ {\left\{\left({X}_i,{Y}_i\right)\right\}}_{i=1}^{n+1} $ are drawn exchangeably from the joint distribution.

The procedure for constructing the prediction intervals is straightforward. First, a regression model is trained on a proper training subset indexed by $ {\mathcal{I}}_1\subset \left\{1,\dots, n\right\} $ to predict $ \hat{\mu}(x) $ . Then, a disjoint calibration subset indexed by $ {\mathcal{I}}_2\subset \left\{1,\dots, n\right\} $ is used to compute the nonconformity scores

(2.2) $$ {R}_i=\mid {Y}_i-\hat{\mu}\left({X}_i\right)\mid i\in {\mathcal{I}}_2, $$

i.e., the residuals. Computing

(2.3) $$ {Q}_{1-\alpha}\left(R,{\mathcal{I}}_2\right):= \left(1-\alpha \right)\left(1+\frac{1}{\mid {\mathcal{I}}_2\mid}\right)\hbox{-} \mathrm{th}\hskip0.42em \mathrm{empirical}\ \mathrm{quantile}\ \mathrm{of}\;{\left\{{R}_i\right\}}_{i\in {\mathcal{I}}_2}, $$

with $ \mid {\mathcal{I}}_2\mid $ the cardinality of the calibration set, yields the prediction interval

(2.4) $$ \mathcal{C}\left({X}_{n+1}\right)=\left[\hat{\mu}\left({X}_{n+1}\right)-{Q}_{1-\alpha}\Big(R,{\mathcal{I}}_2\Big),\hat{\mu}\left({X}_{n+1}\right)+{Q}_{1-\alpha}\Big(R,{\mathcal{I}}_2\Big)\right], $$

which satisfies (2.1). Setting $ \alpha =0.1 $ corresponds to 90% prediction intervals.

The coverage property (2.1) is marginal, i.e., averaged over all samples $ {\left\{\left({X}_i,{Y}_i\right)\right\}}_{i=1}^{n+1} $ (Romano et al., Reference Romano, Patterson and Candes2019; Feldman et al., Reference Feldman, Bates and Romano2021). The coverage guarantee holds for the distribution of the test point $ {X}_{n+1} $ , but may not hold for the conditional distribution of $ {Y}_{n+1} $ given $ {X}_{n+1} $ ,

(2.5) $$ \mathrm{\mathbb{P}}\left\{{Y}_{n+1}\in \mathcal{C}\left({X}_{n+1}\right)\hskip0.4em |\hskip0.7em {X}_{n+1}\right\}\ge 1-\alpha, $$

which is referred to as conditional coverage and, in general, cannot be achieved (Lei and Wasserman, Reference Lei and Wasserman2012). As a result, the coverage probability can fall below $ 100\left(1-\alpha \right)\% $ for subsets of the data, such as certain regimes or ranges of the target variable.

2.2. Conformalized quantile regression

While SCP allows for constructing marginally valid prediction intervals in a straightforward way at almost no additional computational cost, Eq. (2.4) reveals that the length of the prediction interval $ \mathcal{C}\left({X}_{n+1}\right) $ , is fixed to $ 2{Q}_{1-\alpha}\left(R,{\mathcal{I}}_2\right) $ , independent of the input $ {X}_{n+1} $ . SCP is based on the implicit assumption that the spread of the residuals is constant for all inputs $ {X}_{n+1} $ , i.e., homoscedasticity, which is in practice often not the case. To address this issue, Romano et al. (Reference Romano, Patterson and Candes2019) proposed CQR. In CQR, a QR model is used to derive prediction intervals that can handle heteroscedasticity and thus are potentially more informative.

QR (Koenker and Bassett, Reference Koenker and Bassett1978) is a statistical method to estimate conditional quantiles of $ {Y}_{n+1} $ given $ {X}_{n+1} $ , instead of the conditional mean. Here, a QR model is trained on the proper training dataset to predict the lower and upper quantile $ {\hat{Q}}_{\alpha_{\mathrm{lo}}}(x) $ and $ {\hat{Q}}_{\alpha_{\mathrm{hi}}}(x) $ , which constitute initial estimates of the lower and upper bound of the prediction interval, $ \hat{\mathcal{C}}(x)=\left[{\hat{Q}}_{\alpha_{\mathrm{lo}}}(x),{\hat{Q}}_{\alpha_{\mathrm{hi}}}(x)\right] $ , with $ {\alpha}_{\mathrm{lo}}=\alpha /2 $ and $ {\alpha}_{\mathrm{hi}}=1-\alpha /2 $ , e.g., $ \hat{\mathcal{C}}(x)=\left[{\hat{Q}}_{0.05},{\hat{Q}}_{0.95}\right] $ . To quantify the error of this ad hoc prediction interval, nonconformity scores

(2.6) $$ {E}_i:= \max \left\{{\hat{Q}}_{\alpha_{\mathrm{lo}}}\left({X}_i\right)-{Y}_i,{Y}_i-{\hat{Q}}_{\alpha_{\mathrm{hi}}}\left({X}_i\right)\right\}, $$

are computed on the calibration dataset for each $ i\in {\mathcal{I}}_2 $ . For a new input $ {X}_{n+1} $ , the prediction interval is conformalized by computing $ {Q}_{1-\alpha}\left(E,{\mathcal{I}}_2\right) $ analogously to SCP (2.3). This yields the conformalized prediction interval.

(2.7) $$ \mathcal{C}\left({X}_{n+1}\right)=\left[{\hat{Q}}_{\alpha_{\mathrm{lo}}}\left({X}_{n+1}\right)-{Q}_{1-\alpha}\Big(E,{\mathcal{I}}_2\Big),{\hat{Q}}_{\alpha_{\mathrm{hi}}}\left({X}_{n+1}\right)+{Q}_{1-\alpha}\Big(E,{\mathcal{I}}_2\Big)\right]. $$

3.1. Microphysical process rates

We aim to reconstruct MPRs from output of the ICON model, which were not included in the model output at the time of simulation, using a two-moment bulk microphysics scheme.

Bulk parameterization schemes describe cloud microphysical properties within each model grid volume with statistical bulk quantities, i.e., moments of the size distribution functions of distinct hydrometeor categories. The two-moment microphysics scheme (Seifert and Beheng, Reference Seifert and Beheng2006) parameterizes cloud microphysics with prognostic number concentrations and mass mixing ratios of cloud droplets, raindrops, cloud ice, snow, graupel, and hail, corresponding to the first two moments of the respective size distribution functions. MPRs describe a change in the prognostic variables through phase transitions and interaction processes (Figure 2). (In this work, we only consider process rates related to a change in the mass mixing ratio.) Here, we consider two warm-rain processes, which both describe an increase in the rain mass mixing ratio through the formation of raindrops from coalescing cloud droplets (autoconversion) and cloud droplets and raindrops (accretion). Evaporation of rain causes a decrease in rainwater mass. Rain mass also decreases by freezing of rain to snow, graupel, and hail at temperatures $ T<0{}^{\circ}\mathrm{C} $ . At $ T>0{}^{\circ}\mathrm{C} $ , all frozen hydrometeors can melt to rain, leading to an increase in the rain mass mixing ratio. Furthermore, we consider riming, whereby liquid particles freeze upon contact with a frozen hydrometeor. The riming rate comprises eight individual processes between each category of hydrometeors in the liquid and solid phase. The melting, freezing, and riming rates represent the sum of all individual processes. As inputs to the ML models, we use mass mixing ratios and number concentrations together with temperature, pressure, and density. (These are listed in the Supplementary Material A.2. A short discussion of exchangeability can be found in Supplementary Material A.3.)

Figure 2.

Schematic of selected microphysical processes (arrows) between the six hydrometeor categories (boxes) and water vapor in the two-moment microphysics scheme of Seifert and Beheng (Reference Seifert and Beheng2006). Horizontal (green) arrows represent interaction processes, vertical (blue) arrows represent phase transitions. For simplicity, riming is shown separately. If a process occurs more than once, arrows represent contributions to the total process rate. From an ML perspective, boxes represent input features (mass mixing ratios and number concentrations) and arrows represent targets (process rates).

Figure 2. Long description

The flowchart is organized into three vertical tiers.

At the top, the Gas phase contains a box for Water vapor.

In the middle, the Liquid phase contains boxes for Cloud droplets and Raindrops. Two horizontal green arrows point from Cloud droplets to Raindrops, labeled Accretion and Autoconversion. A curved green arrow loops from Raindrops back to itself. A blue arrow labeled Evaporation points from Raindrops up to Water vapor.

At the bottom, the Solid phase contains four boxes: Cloud ice, Snow, Graupel, and Hail. Blue arrows represent phase transitions between these and the liquid phase:

* A blue arrow labeled Melting points from Cloud ice to Raindrops.

* Two blue arrows connect Snow and Raindrops, labeled Freezing (upward) and Melting (downward).

* Two blue arrows connect Graupel and Raindrops, labeled Freezing (upward) and Melting (downward).

* Two blue arrows connect Hail and Raindrops, labeled Freezing (upward) and Melting (downward).

On the far right, a separate simplified diagram shows a box for Liquid hydrometeors above a box for Frozen hydrometeors. A blue arrow labeled Total riming points downward from liquid to frozen, and a curved green arrow loops from the Frozen hydrometeors box back to itself.

A legend in the top right corner identifies blue vertical arrows as Phase transitions and green horizontal arrows as Interactions.

3.2. ICON model simulations and sampling procedure

For training, calibration, and evaluation, we use a dataset compiled from simulation output of the ICON model (version 2.6.6) (Zängl et al., Reference Zängl, Reinert, Rípodas and Baldauf2015) in a limited-area configuration with $ \overline{\Delta x}\approx 2 $ km effective grid spacing. The simulations are performed with the two-moment microphysics scheme (Seifert and Beheng, Reference Seifert and Beheng2006) with a fast-physics timestep $ {\mathrm{t}}_{\mathrm{fast}}=20\hskip0.1em \mathrm{s} $ . The modeling domain extends from $ -0.4{}^{\circ} $ W to $ 17.7{}^{\circ} $ E and from $ 43.7{}^{\circ} $ N to $ 57.26{}^{\circ} $ N, corresponding to the operational ICON-D2 configuration (Reinert et al., Reference Reinert, Prill, Frank, Denhard, Baldauf, Schraff, Gebhardt, Marsigli, Förstner, Zängl, Schlemmer, Blahak and Welzbacher2025). Each vertical column contains 65 levels up to 22 km above ground. The dataset comprises 19 one-day simulations (24 hours each) from January 2022 to July 2023 with an output time step of 10 minutes. We perform simulations for 1 day for each month. For training purposes, the MPRs are added as additional output variables in our setup. Further details on the model configuration are given in Supplementary Material A.1. The set of input features for each target variable is aligned with the computation of the process rate in the ICON subroutine and listed in Supplementary Material B.1. The size of the initial dataset is reduced by random sampling of model grid points. Furthermore, the output from high altitudes is discarded under the assumption that they are cloud-free. We use the simulation output for 2022 for training (even-numbered months) and validation (odd-numbered months); the simulated days for 2023 are used for the test dataset. The training dataset is randomly split into the proper training and the calibration set. The datasets for each target MPR are retrieved from these initial datasets by filtering out all grid points where the value of the MPR fulfills $ \left|\mathrm{MPR}\right|\ge {10}^{-12}\hskip0.1em \mathrm{kg}\;{\mathrm{kg}}^{-1}\hskip0.1em {{\mathrm{t}}_{\mathrm{fast}}}^{-1} $ and $ {\sum}_k{q}_k\ge {10}^{-12}\hskip0.1em \mathrm{kg}\;{\mathrm{kg}}^{-1} $ over all mass mixing ratios $ {q}_k $ in the set of input features, which corresponds to the threshold $ {q}_{\mathrm{crit}}={10}^{-12}\hskip0.1em \mathrm{kg}\;{\mathrm{kg}}^{-1} $ in the two-moment microphysics scheme. (In the ICON model, $ {q}_{\mathrm{crit}} $ is a critical mass mixing ratio threshold used in the parameterization of certain processes, such as autoconversion.) This allows us to train the regression models to predict logarithmically transformed targets, as the range of the target variables spans multiple orders of magnitude. Thus, the proper training dataset contains $ 1.12\times {10}^7 $ samples, and each validation, calibration, and evaluation dataset contains $ 1.6\times {10}^6 $ samples, corresponding to a 70%/10%/10%/10% split. The input features are scaled to the range $ \left[0,1\right] $ with min-max scaling.

3.3. Training and calibration

A separate model is trained for each process rate. In order to obtain prediction intervals with SCP, we train random forest (RF), gradient boosting (XGB), and neural network (NN) models to obtain point predictions for the process rates. For this, we use the scikit-learn library for the RF models (Pedregosa et al., Reference Pedregosa, Varoquaux, Gramfort, Michel, Thirion, Grisel, Blondel, Prettenhofer, Weiss, Dubourg, Vanderplas, Passos, Cournapeau, Brucher, Perrot and Duchesnay2011), the XGBoost library (Chen and Guestrin, Reference Chen and Guestrin2016) for the XGB models, and PyTorch (Paszke et al., Reference Paszke, Gross, Massa, Lerer, Bradbury, Chanan, Killeen, Lin, Gimelshein, Antiga, Desmaison, Köpf, Yang, DeVito, Raison, Tejani, Chilamkurthy, Steiner, Fang, Bai and Chintala2019) for the NN models. Calibration is performed with the absolute error residual (2.2) and $ \alpha =0.1 $ . For CQR, we employ three QR architectures: quantile RF (QRF) models using quantile_forest’s RandomForestQuantileRegressor (Meinshausen, Reference Meinshausen2006; Johnson, Reference Johnson2024), quantile gradient boosting (QXGB), and quantile feed-forward neural networks (QNN). We use the Adam optimizer (Kingma and Ba, Reference Kingma and Ba2017) for training the QNN with a learning rate scheduler with a minimum learning rate of $ {10}^{-6} $ . The maximum number of epochs is 150. The QNN is trained with the quantile loss (Koenker and Bassett, Reference Koenker and Bassett1978), frequently also referred to as the pinball loss,

(3.1) $$ {L}_{\alpha}\left(y,\hat{y}\right)=\alpha \left(y-\hat{y}\right)\unicode{x1D7D9}\left\{y>\hat{y}\right\}+\left(1-\alpha \right)\left(\hat{y}-y\right)\unicode{x1D7D9}\left\{y\le \hat{y}\right\}, $$

where $ \hat{y} $ is the model’s prediction of the “true” $ y $ . All three QR models are trained to predict the $ 5\% $ - and $ 95\% $ -quantiles. More details on the training and the model hyperparameters are given in Supplementary Material B.2. The predicted upper and lower quantiles are calibrated with the nonconformity score function (2.6) and $ \alpha =0.1 $ .