arXiv CSopen access

Uncertainty quantification via conformal prediction in data assimilation

ORIGINAL ARTICLE arXiv

arXiv:2606.27001v1 [cs.LG] 25 Jun 2026

Uncertainty quantification via conformal prediction in data assimilation Catherine George1 |

| Alireza Javanmardi2

| Tijana Janjić3

Eyke Hüllermeier4

1,3 Mathematisches Institut für Maschinelles

Lernen und Data Science, Katholische Universität Eichstätt-Ingolstadt

Quantifying the evolution of uncertainty is critical to both prob-

2,4 Institute of Informatics, LMU Munich,

abilistic forecasting and data assimilation in numerical weather

Munich Center for Machine Learning (MCML) Correspondence Catherine George, Mathematisches Institut für Maschinelles Lernen und Data Science, Katholische Universität Eichstätt-Ingolstadt, Ingolstadt, Bavaria, 85049, Germany Email: [email protected]

prediction. In this study, we investigate the applicability of conformal prediction (CP), a recent machine learning (ML) method, to quantify uncertainty in a controlled, idealized setting. We use the one dimensional modified shallow water model, designed to mimic the convective process. CP provides a set of possible outcomes with a chosen confidence level. Here, we compare

Funding information Klaus Tschira Stiftung gGmbH Projekt-Nr./FKZ: 00.019.2024: "Uncertainty-aware and physics-informed machine learning for short-range atmospheric forecasts".

and evaluate the average empirical coverage, the average interval length, miss low, miss high and average interval score loss (AISL) for three variants of CP, namely a) Standard CP, b) Normalized CP and c) Conformalized Quantile Regression. We further compare these CP-based uncertainty estimates with traditional ensemble-based measures such as standard deviation intervals and ensemble spread. In addition, we investigate the integration of CP-derived uncertainty within the data assimilation cycle through CP perturbations. Our results highlight the strengths and limitations of each approach, providing insight into the effectiveness of CP to complement common ensemble-based uncertainty quantification in simplified atmospheric models. Keywords – Data Assimilation, Uncertainty Quantification, Ensemble Forecasting, Machine Learning, Conformal Prediction.

1

2

1

C.G EORGE ET AL .

|

INTRODUCTION

Uncertainty is an inseparable part of atmospheric prediction. Due to the chaotic dynamics of the atmosphere, the limitations of numerical models, and the incompleteness of observations, no single forecast can ever be fully certain. Quantifying this uncertainty is crucial for informed decision-making (Leutbecher and Palmer, 2008), especially in high-impact weather events where confidence intervals can be as important as the forecast itself. In operational meteorology, uncertainty estimation plays a key role in risk assessment, early warning systems, and decision support for extreme events. In numerical weather prediction, uncertainty estimates are also required in data assimilation, which combines information from numerical models and observations to approximate the most likely atmospheric state (Kalnay, 2002; Evensen, 2009). Data assimilation methods aim to optimally merge prior model forecasts with observations while accounting for uncertainties in both sources of information. The uncertainty arises from several components, including errors in the background forecast, observational errors, and model imperfections. Accurately quantifying these uncertainties is essential, as they determine the relative weighting of model forecasts and observations during the assimilation process. Among data assimilation approaches, ensemble-based data assimilation methods such as the ensemble Kalman filter (EnKF) represent uncertainty by evolving multiple realizations of the atmospheric state (Evensen, 1994; Burgers et al., 1998). In these methods, uncertainty in the background forecast is characterized through the ensemble spread, which provides a sample-based estimate of the forecast error covariance (Evensen, 2003). Traditionally, ensemble methods have therefore been the cornerstone of uncertainty quantification in numerical weather prediction (Leutbecher and Palmer, 2008). By evolving multiple model realizations from perturbed initial states, ensembles approximate the distribution of possible outcomes (Magnusson et al., 2008). Despite their success, ensemble-based approaches face important challenges. In particular, reliable estimation of forecasting uncertainty requires sufficiently large ensembles, which makes these methods computationally demanding (Chattopadhyay et al., 2023; Janjić et al., 2023). Moreover, ensemble sizes are often limited in operational settings due to computational constraints, which can lead to sampling errors and underestimation of uncertainty in the estimated error covariances (Evensen, 2009). These challenges motivate the exploration of alternative approaches for uncertainty quantification. Motivated by the computational expensiveness and slow convergence of Monte Carlo based ensemble methods in highdimensional systems, uncertainty-aware machine learning (ML) is emerging as a promising alternative or complement to ensembles (Price et al., 2025). ML methods, particularly deep neural networks, are powerful function approximators capable of learning complex nonlinear relationships directly from data (LeCun and Hinton, 2015; Goodfellow et al., 2016; Legler and Janjić, 2022). In recent years, ML has increasingly been explored in meteorology and climate science for tasks such as parameterization and forecasting (Reichstein et al., 2019; Rasp et al., 2018). Beyond improving predictive accuracy, several ML approaches aim to explicitly represent predictive uncertainty. Many ML methods quantify uncertainty through prediction intervals or probabilistic forecasts, which describe a range of plausible future states rather than a single deterministic prediction. For instance, Bayesian neural networks (BNNs) model uncertainty in network weights (Neal, 1996; Goan and Fookes, 2020), while deep ensembles use variability across independently trained models as a proxy for model uncertainty (Lakshminarayanan et al., 2017). Other approaches, such as dropout-based Bayesian approximations (Gal and Ghahramani, 2016) or random forests (Breiman, 2001), can also produce probabilistic predictions. These methods provide valuable insights into predictive uncertainty but typically rely on heuristic assumptions or approximate inference schemes. However, ML-based uncertainty quantification is not without challenges. Many probabilistic predictions are miscalibrated, systematically over or underestimating uncertainty (Guo et al., 2017; Arrieta-Ibarra et al., 2022). As a result, prediction intervals may fail to contain the true state at the desired confidence level. In other words, most ML models lack coverage guarantees, meaning that the probability that the true value lies within the predicted interval is not controlled. In applications such as short-range weather forecasting and data assimilation, where reliable uncertainty estimates are essential, this limitation can lead to misleading confidence in predictions. To overcome this gap, conformal prediction (CP) offers a statistically rigorous and distribution-free framework for uncertainty

C.G EORGE ET AL .

3

quantification (Angelopoulos and Bates, 2022). CP produces prediction intervals with guaranteed finite sample coverage, which means that the true value lies within the set at a user-specified confidence level. Unlike ensemble spread or heuristic ML uncertainty estimates, CP provides formal reliability guarantees under very mild assumptions. Recent work has demonstrated the potential of CP for uncertainty quantification in machine learning applications (Romano et al., 2019; Javanmardi and Hüllermeier, 2023). Because conformal prediction does not assume a particular data distribution, it is particularly attractive for atmospheric applications where distributions are often complex and non-Gaussian. To test and develop new uncertainty quantification methods, researchers often use toy models. In our study, we use the one dimensional modified shallow water (MSW) model of Würsch and Craig (2014), later employed by Haslehner et al. (2016) and Ruckstuhl et al. (2021), for the study of data assimilation methods. The model captures essential convective scale dynamics while remaining computationally affordable. The MSW model reproduces regime switches and rainfall onset delays that resemble convection, making it an ideal benchmark for uncertainty-aware methods. Previous work has demonstrated the importance of enforcing physical constraints in data assimilation (Janjić et al., 2014; Zeng et al., 2017; Ruckstuhl and Janjić, 2018) and explored ML-data assimilation hybrids such as convolution neural network-based corrections to conserve mass (Ruckstuhl et al., 2021). Building on this foundation, the present study introduces CP as a novel post-hoc uncertainty quantification tool for the MSW model. Our contributions are: (i) We adapt CP to quantify uncertainty in the spatio-temporal setting of a convective scale toy model, (ii) we employ different CP algorithms and evaluate and compare their coverage, interval sharpness, miss low, miss high rates and average interval score loss, (iii) similarly we compare the different CP methods with traditional uncertainty quantification methods of standard deviation interval and ensemble spread and (iv) we explore the potential of CP to improve data assimilation by introducing CP perturbations either through the convolutional neural network (CNN) output prior to assimilation or directly within the EnKF analysis. The paper is organized as follows. Section 2 introduces the MSW model, data assimilation, training and architecture of the CNN. Section 3 covers uncertainty quantification by CP and describes CP variants along with its experimental setup. Section 4 presents the approach for incorporating CP in data assimilation to improve it. Section 5 discusses the results and Section 6 concludes the paper.

2

2.1

EXPERIMENTAL SETUP

|

|

Modified Shallow Water Model

In this study, we use the one dimensional MSW model (Würsch and Craig, 2014), adopting the numerical implementation as proposed therein. The model captures the dynamics of horizontal wind velocity ’u’, fluid height ’h’ and rain water content ’r’ and has been used as an idealized framework for studying data assimilation algorithms on convective scale (Ruckstuhl and Janjić, 2018; Janjić et al., 2021). Its design allows for the convection representation through simple threshold based triggering and localized perturbations. It ensures conservation of mass, i.e., the spatial integral of fluid height over the domain is constant in time, and the rain is nonnegative. If rain becomes negative during model integration, negative values are set to zero. We adopt the configuration used in the earlier studies (Ruckstuhl et al., 2021), with a 125 km one-dimensional domain discretized into 250 grid points and integrated forward in time with a step size of 5 seconds. The MSW model equations are given by:

∂u ∂u ∂ (φ + γ 2 r ) ∂ 2u +u + = βu + Du , ∂t ∂x ∂x ∂x 2

(1)

4

C.G EORGE ET AL .

with ( φ=

φc

if h > h c

gh

otherwise,

∂r ∂r ∂ 2r +u = Dr − ηr − ∂t ∂x ∂x 2

(

∂u δ ∂x ,

∂u ifh > h r and ∂x <0

0

otherwise,

∂h ∂ (uh ) ∂ 2h + = Dh . ∂t ∂x ∂x 2

(2)

(3)

(4)

The parameters that control the switch between regimes are fixed and have the following values φc = 899.77 m 2 s −2 , h c = 90.02 m, h r = 90.4 m. In the equation for the evolution of rain, the parameter η = 2.5 · 10 −4 is the removal rate of the rain while δ = 1/300 is the production rate. The diffusion parameters are given by Du = D h = 25000 m 2 s −1 , D r = 200 m 2 s −1 for p u, h and r , respectively. The absolute fluid level is h 0 = 90 m and γ = g h 0 is the gravity wave speed where g is the gravity constant. To initiate convection, stochastic Gaussian forcing βu is added at random locations in each model time step. The βu has a half width of four grid points and an amplitude of 0.002 m/s 2 .

2.2

|

Data assimilation and training the CNN

The experimental setup follows a twin-experiment design.The MSW model run will represent the true state of the system, and we call it the nature run. We generate synthetic observations from a nature run by perturbing the nature run with Gaussian errors for wind velocity and fluid height, and log-normal errors for precipitation. In addition, observations are taken only at grid points where precipitation exceeds 0.005 dBZ, with an additional 10% of grid points containing wind data to mimic supplementary measurements. It should be noted that the nature run is used solely for the generation of observations and the evaluation of results and is not involved in any part of the data assimilation processes. Two data assimilation methods are used to create analysis: the ensemble Kalman filter (EnKF; Evensen, 1994; Burgers et al., 1998; Evensen, 2009) and quadratic programming ensemble (QPEns; Janjić et al., 2014) that preserves mass and non-negativity of rain. The EnKF produces an ensemble of analyses that are propagated forward in time during each assimilation cycle using the MSW model. This process yields a background ensemble at fixed time k designed to represent the uncertainty in the model forecasts. The stochastic EnKF uses this ensemble of predictions (backgrounds) {Zb,i }iN=1 to calculate the sample background error covariance matrix N

Pb =

1 Õ b,i b b [Z − Z̄ ] [Zb,i − Z̄ ]T N −1

(5)

i =1

b

with Z̄ representing the ensemble mean. The EnKF then applies the standard Kalman filter equations to compute the analysis for each ensemble member i = 1, . . . , N using Pb and observation error covariance R. The QPEns does not use Kalman filter equations; instead, it minimizes the cost function with additional physical constraints in each ensemble member and uses Pb and R to weight the background and observations, respectively. It is shown to be more accurate than EnKF for state estimation, due to the explicit additional use of physical information (Janjić et al., 2014; Zeng et al., 2017; Ruckstuhl and Janjić, 2018; Gleiter et al., 2022). However, QPEns is computationally expensive compared to EnKF, which motivated Ruckstuhl et al. (2021) to train a convolutional neural network (CNN) to map EnKF analysis to QPEns analysis, thereby achieving QPEns-like quality

C.G EORGE ET AL .

5

at a fraction of the computational cost. Both QPEns and EnKF produce an ensemble of background and analysis in each data assimilation cycle. Following Ruckstuhl et al. (2021), using an ensemble of N = 10 members initialized from arbitrary states, we start the data assimilation cycling and train a CNN that inputs EnKF analysis ensemble and outputs QPEns analysis ensemble. Also, the CNN is not a standard black-box network; it incorporates a mass conservation constraint to ensure that the physical consistency is preserved. We generate training data using 10 random seeds, consisting of QPEns cycling outputs Q tb ,Q ta : t = 1,2,...,200, where Q denotes the QPEns system, and the superscripts b and a represent the background and analysis states, respectively. Simultaneously, we construct the dataset Xta : t = 1,2,...,200, where Xta denotes the unconstrained solution obtained from the same background Q tb . For each random seed, the data lies in ÒN ×n ×3 , where N is the number of ensemble members, n is the number of grid points (n = 250), and the three channels correspond to the variables u, h, and r, over 200 assimilation cycles. The corresponding output data has the same dimensional structure, with the CNN applied to the EnKF analysis for the subsequent 180 assimilation cycles, starting from timestep t = 20.

The CNN is trained to learn the correction introduced by the constrained QPEns analysis relative to the unconstrained EnKF update. To prevent the influence of large transient increments associated with the initial spin-up phase, the first 20 assimilation cycles are excluded from the training dataset. QPEns cycling is performed initially, and only states obtained from cycle 20 onward, when the system has reached a dynamically balanced state, are used to construct the training pairs. To generate the training data set, QPEns cycling produces background states Q tb and corresponding constrained analyses (Q ta ). At each assimilation timestep, the same QPEns background (Q tb ) is used to compute an unconstrained EnKF analysis (Z ta ) by applying the standard Kalman update without enforcing physical constraints. The CNN is then trained to map the unconstrained analysis (Z ta ) to the constrained QPEns analysis (Q ta ). By using identical background states for both updates, the training procedure isolates the difference between the unconstrained and constrained minimization steps, ensuring that the CNN learns only the effect of constraint enforcement, rather than differences caused by accumulated background errors during cycling. Wind velocity ’u’ and fluid height ’h’ are normalized by subtracting their climatological mean and dividing by the standard deviation, while precipitation ’r ’ is scaled only by its standard deviation to preserve positivity. The overall CNN training pipeline adopted in this study is illustrated schematically in Figure 1. For further details, we refer to Ruckstuhl et al. (2021). Zat EnKF

MSW

Qb t

QPEns

FIGURE 1

CNN Training

Observations

Trained CNN

Qat

Schematic overview of the CNN training workflow within the data assimilation framework.

For the evaluation experiments, QPEns is first integrated independently for 20 assimilation cycles to allow the system to reach a dynamically balanced state. The analysis ensemble obtained at cycle 20 is then used as a common initial condition for the EnKF and CNN corrected EnKF experiments. From cycle 21 onward, three independent data assimilation systems are run: EnKF, QPEns, and EnKF augmented with the trained CNN correction. Each method subsequently performs its own forecast and analysis steps without further interaction, ensuring a consistent and fair comparison of their respective performances.

6

C.G EORGE ET AL .

2.3

Architecture of the CNN

|

To emulate the relationship between the unconstrained and physically constrained analyses, CNN is employed. As in Ruckstuhl et al. (2021), the CNN consists of four convolutional hidden layers, each containing 32 filters with a kernel size of 3. The activation function used in all hidden layers is the scaled exponential linear unit (SELU), which helps maintain normalized activations throughout the network and improves training stability. The output layer is also a convolutional layer, with three filters corresponding to the three model variables - u, h, r and a kernel size of 3. A linear activation is applied to the output channels for u and h, while a ReLU activation ensured non-negative predictions for r. The effective influence radius of the CNN, is five grid points. Training is conducted using the Adam optimizer, and the loss function is defined as the root mean square error (RMSE) between the predicted and target analyses, computed over all grid points and averaged over the three variables u, h, and r. To prevent overfitting and monitor generalization, a validation dataset is used alongside the training data, and the model is trained for 100 epochs with a batch size of 96. The implementation is carried out in Python using Keras (Chollet, 2017) and is based on the implementation of Ruckstuhl et al. (2021).

3

|

UNCERTAINTY QUANTIFICATION VIA CP

In the system described above, the CNN takes the EnKF analysis Z ta,new as input and outputs an estimate Q̂ ta,new of the QPEns analysis Q ta,new 1 . To equip this prediction with a reliable measure of uncertainty, we apply the CP framework to construct a distribution-free prediction interval PI(Z ta,new , α ) ⊂ Ò that contains the sought variable Q ta,new with a high probability 1 − α, where α ∈ (0, 1) is a user-defined error rate, as illustrated in Figure 2. In this context, uncertainty is represented by the length of the PI; a wider interval indicates higher uncertainty. At its core, CP performs hypothesis testing. Given an input Z ta,new , it evaluates every candidate output y ∈ Ò to determine if it could plausibly be the valid label. A candidate is only excluded if the hypothesis is rejected. Hypothesis testing is performed by calculating a nonconformity score, a real-valued function S (Z ta,new , y ) that quantifies how "unusual" the pair (Z ta,new , y ) appears relative to previously observed data. In this study, the QPEns analysis is used as a proxy truth for the evaluation of the CP intervals.

Zbt

MSW

FIGURE 2

EnKF

Zat

Trained CNN

Q̂ at

CP

PIt

CP pipeline for generating prediction intervals from CNN analysis.

There are two main variants of conformal prediction: Full CP and Split CP. Full CP is a transductive method, in which the nonconformity score of each candidate label y is computed by retraining the model on the full dataset augmented with the test point (Z ta,new , y ), and comparing it against the scores of all other data points. While this yields tight prediction intervals, it is computationally prohibitive since it requires a full model retraining for every candidate y . Split CP, by contrast, is an inductive method that trains the model once and uses a separate held-out calibration set to estimate the nonconformity score distribution, making it computationally efficient and well-suited to settings where the model is expensive to train. For this reason, we focus on Split CP in this paper. Suppose we are given a dataset of paired observations { (Z ta,i , Q ta,i ) }iM=1 . Split CP begins by randomly partitioning the data indices [M ] into two disjoint subsets: a training set Itrain and a calibration set Ical . A CNN is 1 With a slight abuse of notation, we treat Z a,i and Q a,i (and similarly Z a,new and Q a,new ) as scalar quantities corresponding to a single variable (among r , u, h ) at

t

t

t

t

assimilation time t . The indexing (e.g., i or “new”) is used generically and may refer to any instance among the N × n instances (e.g., across ensemble members or spatial grid points).

C.G EORGE ET AL .

7

trained exclusively on Itrain to produce point predictions Q̂ ta,i . The nonconformity score is then evaluated on the calibration data, yielding the set of scores {S (Z ta,i , Q ta,i ) : i ∈ Ical }. We denote by q̂ 1−α the (1 − α ) (1 + 1/| Ical | ) quantile of this set of scores. The resulting prediction interval is defined as

PI(Z ta,new , α ) = { y ∈ Ò : S (Z ta,new , y ) ≤ q̂ 1−α }.

(6)

If the data points { (Z ta,i , Q ta,i ) }i ∈Ical ∪ (Z ta,new , Q ta,new ) are exchangeable, meaning their joint distribution is invariant under permutation, then the Split CP procedure described above satisfies the marginal coverage guarantee:
Ð Q ta,new ∈ PI(Z ta,new , α ) ≥ 1 − α .

(7)

This guarantee is distribution-free and holds regardless of the choice of nonconformity score. However, the nonconformity score significantly influences the length of the prediction intervals, a property referred to as efficiency, as well as the shape of the intervals, as we will clarify in the following sections. In this paper, we focus on three popular methods in CP: standard CP (SCP) (Vovk et al., 2022), normalized CP (NCP) (Papadopoulos et al., 2008), and conformalized quantile regression (CQR) (Romano et al., 2019). Figure 2 illustrates the proposed framework of integrating CP to the current setup to predict intervals.

3.1

|

Standard CP

The nonconformity score S (Z ta,i , Q ta,i ) is often defined by the level of disagreement between the point prediction Q̂ ta,i and the true value Q ta,i . Since these are real-valued numbers, a standard choice is the absolute residual: S SCP (Z ta,i , Q ta,i ) = | Q̂ ta,i − Q ta,i |.

(8)

Using this nonconformity score, the prediction interval PI(Z ta,new , α ) can be written as: PISCP (Z ta,new , α ) = [ Q̂ ta,new − q̂ 1−α , Q̂ ta,new + q̂ 1−α ].

(9)

An immediate drawback of this approach is that the prediction intervals have a constant length of 2q̂ 1−α for every point. As mentioned previously, the interval length represents the uncertainty of the prediction; therefore, one would expect the length to adapt to the difficulty of the prediction at each point. This constant-width behaviour fails to account for heteroscedasticity, where some regions of the input space are inherently more uncertain than others.

3.2

|

Normalized CP

To address this, more sophisticated nonconformity scores can be used to achieve adaptivity. One such method is normalized CP, which scales the residual by a heuristic uncertainty estimate σ̂ (Z ta,i ): S NCP (Z ta,i , Q ta,i ) =

| Q̂ ta,i − Q ta,i | σ̂ (Z ta,i )

.

(10)

8

C.G EORGE ET AL .

Accordingly, the resulting prediction interval becomes: PINCP (Z ta,new , α ) = [ Q̂ ta,new − q̂ 1−α σ̂ (Z ta,new ), Q̂ ta,new + q̂ 1−α σ̂ (Z ta,new ) ].

(11)

In this case, the length of the prediction interval at each point Z ta,new is proportional to the uncertainty estimate σ̂ (Z ta,new ) at that point. A common approach is to define σ̂ (Z ta,new ) such that it estimates the local dispersion of the residuals at Z ta,new . For this purpose, an auxiliary regressor can be trained on the training data { (Z ta,i , |Q ta,i − Q̂ ta,i | ) : i ∈ Itrain }.

3.3

|

CQR

Another popular approach that also satisfies the requirement for adaptivity is Conformalized Quantile Regression (CQR) (Romano et al., 2019). Unlike previous methods that rely on a single point estimate Q̂ ta,i , CQR leverages training data Itrain to train two quantile regressors, L̂ (Z ta,i ) and Û (Z ta,i ), which are trained to estimate the lower and upper quantiles of the target distribution (e.g., at levels α/2 and 1 − α/2). For our CNN training step, this means we must replace the RMSE loss function with the pinball loss (also known as quantile loss). The network architecture remains identical to the RMSE based model; only the loss function is modified. In the quantile setting, the output layer is adapted to predict the desired quantile levels, and the model is trained using the pinball loss instead of the mean squared error. The nonconformity score for CQR is defined as S CQR (Z ta,i , Q ta,i ) = max{ L̂ (Z ta,i ) − Q ta,i , Q ta,i − Û (Z ta,i ) }.

(12)

This score is positive if the true value Q ta,i falls outside the estimated bounds and negative if it falls inside. The final prediction interval is constructed by shifting the initial quantile estimates: PICQR (Z ta,new , α ) = [ L̂ (Z ta,new ) − q̂ 1−α , Û (Z ta,new ) + q̂ 1−α ].

(13)

Here, the conformal quantile appears as a correction term applied to the interval produced by the quantile regressor. In particular, a positive calibration quantile indicates that, for a substantial fraction of calibration points, the true labels fall outside the predicted quantile intervals. Consequently, the intervals must be expanded to achieve the desired coverage.

3.4

|

CP Experimental Setup

In our experiments, we generate evaluation data using 10 random seeds, similar to what is explained in subsection 2.2. We then randomly select 5 seeds for calibration and the other 5 seeds for testing. We apply CP independently for each of the three variables (u, h, r ) and for each assimilation cycle. This leads to 5 × N × n calibration instances per variable–assimilation cycle, and 5 × N × n test instances. Furthermore, we consider two different variants for applying conformal prediction: ensemble-based and mean-based approaches. In the ensemble-based approach, we consider the entire 5 × N × n instances and apply CP accordingly. In the mean-based approach, we first take the mean across the ensemble size, and then apply CP to the 5 × n instances. In either case, we treat the test sets similarly. Also, as rainfall is a non-negative variable, the lower bound of the prediction interval is clipped to zero to ensure positivity of rain.

C.G EORGE ET AL .

4

9

INCLUSION OF CP INTERVALS IN DATA ASSIMILATION CYCLING

|

After training the CNN and obtaining the uncertainty estimates, we investigated the possibility of using these estimates during data assimilation with EnKF. The uncertainty of the hybrid of a CNN and the EnKF has been represented through the analysis ensemble during data assimilation. Therefore, thus far, the CP intervals are constructed in a post-processing way, meaning that the data assimilation cycle operates independently of CP, while CP utilizes the resulting analysis to quantify uncertainty as shown in Figure 2. In particular, CP provides an interval where we expect the QPEns analysis, used here as a proxy truth or reference, to lie for every grid point across all assimilation cycles. Due to the small ensemble size, the ensemble spread alone does not adequately represent the full analysis uncertainty, motivating the use of CP to provide a more reliable uncertainty quantification. Furthermore, since the CNN is trained to approximate QPEns and not the nature run truth, the CP intervals quantify uncertainty with respect to the QPEns proxy truth, and not to the nature run. Therefore, we additionally perturb each analysis ensemble member assuming a Gaussian distribution over the CP interval in each grid point independently, creating the CP informed analysis ensemble { z̃a,i }iN=1 , i.e., z̃a,i = za,i + η i

(14)

Here η i is a vector of size of the state, whose each element is a random realization of a normally distributed random variable with mean zero and with standard deviation half of the calculated CP interval for each cycle averaged spatially as mentioned in Equation 15. ηi =

Ct i ξ, 2

ξ i ∼ N (0, 1),

(15)

where Ct is the width of the CP interval at each cycle t and i = 1, . . . , N indexes the members of the ensemble. Note that for the rainfall variable, the conformal scores are computed without clipping the lower bound of the prediction interval to zero to preserve the full uncertainty in the calibration process. Figure 3 provides a schematic overview of the proposed workflow, in which CP intervals are injected as perturbations into the data assimilation cycle. The perturbations are applied either to (a) the CNN analysis or (b) the EnKF analysis, followed by EnKF assimilation. The resulting EnKF analysis is then used by CNN for correcting, which yields updated CNN and EnKF analyzes that incorporate conformal uncertainty information. These two configurations represent different ways of integrating CP within the data assimilation framework, moving beyond the traditional use of CP as a purely post-processing uncertainty quantification method. In addition, investigates how incorporating CP estimates directly within the data assimilation cycle can improve the quality of the analysis, rather than being applied only after the prediction stage. Choose a value randomly from Gaussian assumption

a)

Choose a value randomly from Gaussian assumption

PIt

PIt

b)

MSW

Zbt

EnKF

Zat

Trained CNN

Q̂ at

CP

MSW

Zbt

EnKF

Zat

Trained CNN

Q̂ at

CP

F I G U R E 3 Schematic overview of the integration of CP with the CNN-EnKF data assimilation framework, showing perturbation injection to the (a) CNN output and (b) EnKF analysis stages

10

5

C.G EORGE ET AL .

RESULTS

|

We repeat the random splitting of the evaluation seeds into calibration and test, as explained in subsection 3.4, 10 times, each time performing calibration with the two variants (i.e., ensemble-based and mean-based), for the three algorithms (SCP, NCP, and CQR), for each variable–assimilation cycle, and then evaluating the prediction sets on the test data. For a given test data point (Z ta,new , Q ta,new ) and a prediction interval PI(Z ta,new , α ), we compute the coverage as 1(Q ta,new ∈ PI(Z ta,new , α ) ), the interval size as max PI(Z ta,new , α ) − min PI(Z ta,new , α ), miss low which measures the magnitude by which the true value falls below the lower bound of the prediction interval as min PI(Z ta,new , α ) − Q ta,new · 1{Q ta,new < min PI(Z ta,new , α ) }, miss high which measures the magnitude by which the true value exceeds the upper bound of the prediction interval as Q ta,new − max PI(Z ta,new , α ) · 1{Q ta,new > max PI(Z ta,new , α ) } and the interval score loss (ISL) (Gneiting and Raftery, 2007) as: 

max PI(Z ta,new , α ) − min PI(Z ta,new , α )

2 · min PI(Z ta,new , α ) − Q ta,new · 1{Q ta,new < min PI(Z ta,new , α ) } α 
2 + · Q ta,new − max PI(Z ta,new , α ) · 1{Q ta,new > max PI(Z ta,new , α ) } . α +

(16)

For each random split seed, these values are averaged over all test points to calculate the marginal coverage, average interval size, and average interval score loss (AISL) for the target coverage value 1 − α = 0.90 (α is 0.1 for 90% coverage). Finally, these average values are aggregated across different random split seeds, and the resulting mean and standard deviation are reported accordingly. In the next subsections 5.1-5.3, different CP methods are compared with respect to average coverage and average interval size, the ensemble based approach is considered first.

5.1

|

Uncertainty quantification with SCP

In this approach, the nonconformity score is computed from the absolute error between the CNN prediction and the QPEns analysis, which is considered the reference solution that approximates the truth. As a result, the CP interval directly reflects how far the CNN is expected to deviate from the QPEns. Larger intervals indicate time steps or spatial regions where the CNN predictions become less certain, whereas narrow intervals correspond to high local agreement with QPEns. Figure 4a, shows the average empirical coverage over all assimilation steps for the three prognostic variables u, h, r. Here, average empirical coverage denotes the fraction of grid points whose reference states are contained within the CP intervals, averaged over all ensemble members and random seeds at each time step. Across the variables, the ensemble based SCP consistently achieves coverage close to the nominal 90% level, demonstrating that the calibrated nonconformity distribution yields well-balanced uncertainty intervals. Coverage for the variables u, h, r remains particularly stable but display slightly higher variability around timestep 50 to 75. Nevertheless, no systematic over or under coverage occurs, indicating that SCP provides a robust baseline uncertainty representation across the full assimilation window. The average SCP interval sizes (Figure 4b) reveal how uncertainty evolves over time. The initial jump in the interval size is due to the spin-up time for the data assimilation methods to stabilize, which can be ignored. For u, the interval sizes remain very small (between 4 × 10 −3 and 5.5 × 10 −3 ), indicating increased deviation and less smooth velocity field. The height field (h) displays intervals around 2.5 × 10 −2 and 5 × 10 −2 . The rainfall variable r widths up to roughly 1.5 × 10 −3 and 2.5 × 10 −3 . To illustrate how the SCP intervals influence the physical representation of each state variable, we present spatial snapshots at two representative time steps. The early snapshot (t=25) captures the initial assimilation phase, while the later snapshot (t=160) reflects long-term behaviour after the system has adjusted to repeated assimilation. Together, these snapshots provide a clear

C.G EORGE ET AL .

11

F I G U R E 4 Temporal evolution of the a) average empirical coverage and b) average interval size for SCP method applied to variables u (left), h (middle), r (right) from timestep 20. In a) solid line denote mean ensemble coverage across random seeds, shaded region indicate coverage variability, and the dashed black line represents the nominal 90% target coverage. In b) solid line denote mean ensemble interval size across random seeds and the shaded region indicate interval size variability.

F I G U R E 5 Truth (black dashed) and ensemble mean for QPEns (blue) and NN (orange) at timestep 25 (left) and 160 (right). Gray shading represents uncertainty calculated with SCP Method for the three prognostic variables u (upper), h (middle), r (lower).

visualization of the evolution of uncertainty and physical consistency between variables. The snapshot at timesteps 25 and 160 (Figure 5) shows the ensemble mean of QPEns and CNN as well as the true fields together with the SCP intervals for a randomly selected test sample. For variable r, the interval width is clipped to zero to showcase the positivity of rain. The CP interval for each variable remains essentially constant across all grid points because SCP uses a single global nonconformity threshold and does not adapt locally to variations in spatial scale. As a result, while the interval width is appropriate for most regions, at a few grid points typically near sharper peaks, the interval does not fully enclose the QPEns value, and the truth also lies outside the bounds. This behaviour of truth field not captured fully is consistent with the computation of SCP, which calibrates uncertainty solely based on CNN-QPEns deviations, and also this method cannot adjust its width to local dynamical variability. These limitations motivate the use of other CP methods, which rescales the nonconformity score and allows the interval width to vary across space.

12

5.2

C.G EORGE ET AL .

|

Uncertainty quantification with NCP

The NCP results show how incorporating a scaling term enables the interval widths to adapt locally across both time and space, producing uncertainty bands that better reflect the structure of each variable. Here, the nonconformity score is scaled by a normalization factor ‘σ’, defined as the standard deviation of the CNN ensemble predictions across grid points. The normalization factor σ is estimated using a Random Forest regressor composed of 100 decision trees with a maximum depth of 20. The model is trained separately for each variable to predict the magnitude of the absolute residual between the CNN and QPEns analyses, with five independent random seeds to improve robustness. The resulting σ values are averaged across seeds and used to normalize the residuals in the computation of nonconformity scores, allowing the method to account for heteroscedastic uncertainty. NCP is applied only to the variables u and h, as the rainfall variable r contains many zero valued grid points; using normalization in such regions would lead to unstable or inflated nonconformity scores.

F I G U R E 6 Temporal evolution of the a) average empirical coverage and b) average interval size for NCP method applied to variables u (left), h (middle), r (right) from timestep 20. In a) solid line denote mean ensemble coverage across random seeds, shaded region indicate coverage variability, and the dashed black line represents the nominal 90% target coverage. In b) solid line denote mean ensemble interval size across random seeds and the shaded region indicate interval size variability.

In Figure 6a, the ensemble based NCP maintains average empirical coverage close to the nominal 90% level for all three variables, with noticeably varying intervals compared to the SCP. The average interval sizes (Figure 6b) for the velocity component u is around 2.4 × 10 −3 and 3.2 × 10 −3 , the variable h around 1.5 × 10 −2 and 2.2 × 10 −2 , and the rainfall component near 1.5 × 10 −3 and 2.5 × 10 −3 . These values are smaller than those produced by the SCP because the normalization term compensates for local changes in space. The snapshot of the ensemble mean of a randomly selected test sample at timestep 25 and 160 (Figure 7) further illustrates this behaviour; unlike SCP which imposes a uniform interval width across all grid points, NCP expands the bands in regions with stronger local variability and contracts them in smoother areas, producing spatially responsive uncertainty estimates; although a few sharp peaks remain outside the interval for certain grid points.

5.3

|

Uncertainty quantification with CQR

The CQR framework is constructed by training a quantile regression neural network using the pinball loss at two quantile levels, 1 − α/2= 0.95 and α/2=0.05, enabling the model to learn data driven lower and upper predictive bounds. In Figure 8a, the ensemble based CQR method maintains average coverage close to the prescribed 90% level for all variables but more variability than NCP for variables u and h. The average interval sizes (Figure 8b) follow a pattern like the earlier CP variants: an initial rise during spin-up followed by a rapid contraction after timestep 20. The stabilized widths lie around 3.8 × 10 −3 and 5 × 10 −3

C.G EORGE ET AL .

13

F I G U R E 7 Truth (black dashed) and ensemble mean for QPEns (blue) and NN (orange) at timestep 25 (left) and 160 (right). Gray shading represents uncertainty calculated with NCP Method for the three prognostic variables u (upper), h (middle), r (lower).

for u, approximately 2 × 10 −2 to 3.5 × 10 −2 for h, and roughly 1 × 10 −3 and 2 × 10 −3 for r. The spatial snapshot at timestep

F I G U R E 8 Temporal evolution of the a) average empirical coverage and b) average interval size for CQR method applied to variables u (left), h (middle), r (right) from timestep 20. In a) solid line denote mean ensemble coverage across random seeds, shaded region indicate coverage variability, and the dashed black line represents the nominal 90% target coverage. In b) solid line denote mean ensemble interval size across random seeds and the shaded region indicate interval size variability.

25 and 160 (Figure 9) illustrates how CQR forms locally responsive intervals. The quantile bands widen near regions where the QPEns solution exhibits stronger deviations from the median, while contracting in smoother segments of the domain. The quantile outputs appear in the figure as dashed green (α/2=0.05) and dashed cyan (1 − α/2=0.95) lines, with their midpoint shown in green. For the rainfall variable r, the behaviour is notably different when compared to other methods: because r is highly sporadic and contains many zero valued grid points, the network learns a very narrow lower quantile and a rapidly expanding upper quantile near convective peaks, resulting in visibly asymmetric intervals. This reflects the innate skewness and heavy tailed distribution of the rainfall field rather than instability in the CQR method. The interval for u is noticeably wider than for h, reflecting larger spread in the velocity component, while for h, the quantiles stay relatively close and only widen slightly around the regions of increased variability.

14

C.G EORGE ET AL .

F I G U R E 9 Truth (black dashed) and ensemble mean of QPEns (blue),upper quantile (1 − α/2=0.95,cyan), lower quantile (α/2=0.05,green) and midpoint of quantiles (orange) at timestep 25 (left) and 160 (right). Gray shading represents uncertainty calculated with CQR Method for the three prognostic variables u (upper), h (middle), r (lower).

5.4

Comparison of uncertainty quantification in different methods

|

In this section, the performance of different uncertainty quantification approaches is analyzed through two complementary comparisons as mentioned in subsection 3.4. First, the behaviour of the three CP methods is examined by computing CP intervals using two different strategies: computing CP directly from the ensemble members and computing CP from the ensemble mean. This allows us to evaluate how the method of constructing the CP scores influences the resulting uncertainty estimates and to assess which formulation provides more reliable results. In the second comparison, CP-based uncertainty estimates computed from the ensemble mean are evaluated against traditional uncertainty quantification methods, namely the standard deviation interval and the ensemble mean spread. This comparison provides insight into how CP methods perform relative to conventional approaches and illustrates how CP can be incorporated as an alternative way to quantify predictive uncertainty.

5.4.1

|

Comparison of uncertainty quantification in CP methods

From Table 1 which summarize the AISL values, empirical coverage percentages, interval width and directional error rates (miss-high and miss-low) averaged over the full assimilation window computed against QPEns, the average score loss interval (AISL) matrix (Gneiting and Raftery, 2007) provides a quantitative view of how each CP method behaves across the three prognostic variables. Here, we are focusing on comparing the three CP methods when the conformal scores are computed directly from the ensemble members and from the ensemble mean. Both formulations provide reliable uncertainty estimates with coverage generally close to the nominal 90% target. The comparison is therefore intended not to demonstrate failure of either approach, but rather to identify which formulation offers a more stable and representative implementation of CP-based uncertainty quantification within the ensemble forecasting framework. For SCP, the ensemble-mean formulation generally provides improved uncertainty estimates compared to the ensemble-based implementation. The SCP mean-based approach consistently produces smaller AISL values, such as 0.004±0.001 for u compared to 0.008±0.001 for the ensemble-based formulation, while also maintaining coverage closer to the nominal target (89.93%±4.18% compared to 89.89%±3.02%). The interval widths are also slightly narrower for the mean-based approach (0.002% for u) than for the ensemble-based formulation (0.004%). A similar behaviour is observed for h, where AISL decreases from 0.094±0.027 to 0.055±0.023. For variable r, AISL is lesser for mean based (0.008) than ensemble based (0.012), but interval width is slightly higher by 0.001 with near similar miss rates. Overall, the SCP mean-based formulation provides slightly tighter intervals while maintaining balanced uncertainty

C.G EORGE ET AL .

15

estimates. For NCP, the opposite trend is observed, with the ensemble-based NCP generally providing stronger performance than its mean-based counterpart. For u, the ensemble-based implementation achieves an AISL of 0.003±0.000 with coverage of 89.98%±1.36%, compared to 0.004±0.001 and 90.02%±3.96% for the mean-based approach. A similar pattern is seen for h, where the ensemble-based method produces AISL of 0.021±0.005 compared to 0.052±0.019 for the mean-based method, while maintaining coverage close to the nominal level. The interval widths are also same for both the implementations, and the miss-low and miss-high rates remain relatively symmetric for ensemble based approach and highly asymmetric for mean based approach. But for the variable r, the AISL value for mean based is lower (0.008) but has higher interval width by 0.001 when compared with the ensemble based counterpart. For CQR, the mean-based formulation tends to provide improved performance compared to the ensemble-based version. For instance, for u, the AISL decreases from 0.007±0.001 in the ensemble-based formulation to 0.004±0.001 when using the ensemble mean, while maintaining coverage close to 90%. A similar improvement is observed for h, where the AISL decreases from 0.080±0.026 to 0.051±0.019. For the variable r, mean based have lower AISL values of 0.007 than ensemble based with 0.011 and similar coverage levels. However, the mean-based CQR approach generally produces slightly narrower intervals. To further assess the robustness of the CP-based uncertainty estimates, the same metrics are also evaluated against the true model state rather than QPEns (Table 2). Since the CP intervals are constructed using QPEns as the reference, the empirical coverage with respect to truth is not expected to match the nominal 90% level. Instead, these statistics provide insight into how well the QPEns-based uncertainty representation captures the actual model error. For SCP, the mean-based formulation again shows more consistent performance when evaluated against truth. For u, the mean formulation achieves AISL of 0.010±0.002 but with lower coverage of 64.46%, while the ensemble-based implementation produces a slightly larger AISL (0.011±0.001) and lower coverage (77.99%). The directional miss rates also portray this difference in the ensemble formulation (11.44% miss-low and 10.51% miss-high) compared to the mean-based case (17.99% and 17.55%). For the variable h, similarly mean based have lower AISL values of 0.073 than ensemble based with 0.094, with ensemble-based approach maintaining higher interval width (0.034) than the mean-based formulation (0.024). For the variable r, both approaches produce comparable coverage levels but lower AISL for mean based with slightly higher interval width than ensemble based. For NCP, evaluating against truth reveals larger deviations, particularly for the ensemble-based implementation. For u, the ensemble formulation yields AISL of 0.016±0.001 but only 44.28% coverage, accompanied by large directional miss rates (29.82% miss-low and 25.90% miss-high). The mean-based implementation improves the AISL (0.010±0.001) and increases coverage to 62.83%, although the miss rates remain relatively high. A similar improvement is observed for the variable h, where the ensemble-based formulation achieves 47.41% coverage compared to 80.35% for the mean-based version. For the variable r, however, both implementations perform similarly, with AISL values lower for mean based (0.009), but both approaches with coverage close to 90% and interval width around 0.001-0.002. For CQR, the both implementations demonstrate near similar performances. For u, both formulation produces AISL of 0.011± but ensemble based have coverage of 75.18% coverage, while the mean-based formulation yields lower coverage of 61.70% despite closer AISL values. For h, mean based have lower AISL of 0.068 than ensemble based with 0.078, though the ensemble formulation maintains slightly higher coverage (85.12%) than the mean-based version (84.70%). For the rainfall r, both implementations perform well, with AISL values around 0.009 and coverage around 90%, while maintaining relatively small directional miss rates.

5.4.2

|

Comparison of uncertainty quantification with CP and Traditional methods

Now, we consider traditional ensemble-based methods, namely the standard deviation interval and the ensemble spread, as baseline approaches to compare with ensemble mean based CP methods. The standard deviation interval is scaled to provide a 90% confidence level, ensuring comparability with the nominal coverage of CP methods, and the ensemble spread represents the raw variability of ensemble members around the mean and does not provide an explicit probabilistic coverage guarantee. As a

16

C.G EORGE ET AL .

Method SCP ensemble based

NCP ensemble based

CQR ensemble based

SCP mean based

NCP mean based

CQR mean based

Std Deviation

Ensemble spread

Variable

AISL (mean±std)

Coverage (mean±std)

Interval Width (mean±std)

Miss Low (mean±std)

Miss High (mean±std) 5.17% ± 1.84%

u

0.008 ± 0.001

89.89% ± 3.02%

0.004 ± 0.000

4.94% ± 1.54%

h

0.094 ± 0.027

89.82% ± 3.88%

0.034 ± 0.012

5.22% ± 2.31%

4.95% ± 1.96%

r

0.012 ± 0.005

89.64% ± 3.85%

0.001 ± 0.000

5.71% ± 2.39%

4.65% ± 1.93% 4.83% ± 0.97%

u

0.003 ± 0.000

89.98% ± 1.36%

0.002 ± 0.000

5.19% ± 0.94%

h

0.021 ± 0.005

89.98% ± 1.14%

0.020 ± 0.005

4.73% ± 0.93%

5.29% ± 1.01%

r

0.012 ± 0.005

89.64% ± 3.85%

0.001 ± 0.000

5.71% ± 2.39%

4.65% ± 1.93% 4.29% ± 1.67%

u

0.007 ± 0.001

89.87% ± 2.99%

0.004 ± 0.000

5.84% ± 1.77%

h

0.080 ± 0.026

89.84% ± 3.62%

0.026 ± 0.009

3.49% ± 2.09%

6.66% ± 2.27%

r

0.011 ± 0.000

89.84% ± 3.63%

0.001 ± 0.000

3.02% ± 1.60%

7.14% ± 2.76%

5.34% ± 2.81%

u

0.004 ± 0.001

89.93% ± 4.18%

0.002 ± 0.000

4.73% ± 2.33%

h

0.055 ± 0.023

89.93% ± 4.32%

0.024 ± 0.007

5.36% ± 2.87%

4.71% ± 2.35%

r

0.008 ± 0.004

89.85% ± 4.20%

0.002 ± 0.000

5.49% ± 2.76%

4.66% ± 2.18% 3.87% ± 2.35%

u

0.004 ± 0.001

90.02% ± 3.96%

0.002 ± 0.000

6.11% ± 2.60%

h

0.052 ± 0.019

89.88% ± 4.01%

0.027 ± 0.007

2.06% ± 1.56%

8.06% ± 3.25%

r

0.008 ± 0.004

89.85% ± 4.20%

0.002 ± 0.000

5.49% ± 2.76%

4.66% ± 2.18% 3.64% ± 2.26%

u

0.004 ± 0.001

89.85% ± 4.26%

0.002 ± 0.000

6.50% ± 3.14%

h

0.051 ± 0.021

89.90% ± 3.92%

0.020 ± 0.006

2.20% ± 2.01%

7.89% ± 3.53%

r

0.008 ± 0.003

89.81% ± 4.21%

0.002 ± 0.000

2.56% ± 1.80%

7.63% ± 3.34%

1.65% ± 1.22%

u

0.005 ± 0.000

95.71% ± 2.15%

0.004 ± 0.000

2.64% ± 1.45%

h

0.051 ± 0.019

88.80% ± 3.84%

0.025 ± 0.005

2.33% ± 1.66%

8.88% ± 3.14%

r

0.007 ± 0.003

31.79% ± 8.07%

0.001 ± 0.000

1.83% ± 1.34%

66.38% ± 8.48%

2.33% ± 1.53%

u

-

94.53% ± 2.48%

0.004 ± 0.000

3.12% ± 1.58%

h

-

87.13% ± 4.19%

0.024 ± 0.005

3.74% ± 2.33%

9.13% ± 3.17%

r

-

31.83% ± 7.99%

0.002 ± 0.000

2.55% ± 1.66%

65.62% ± 8.62%

TA B L E 1 Summary of AISL, empirical coverage percentage, interval width and directional error rates percentage (miss-high and miss-low) for the three CP methods (both considering all ensembles and mean based), standard deviation method and ensemble mean spread method averaged over the full assimilation window computed against QPEns.

result, AISL value for ensemble spread method is not able to be computed. From (Table 1), when comparing the mean-based CP methods against the standard deviation interval using QPEns as the reference, all CP approaches demonstrate competitive and, in some cases, improved performance. For u, all three CP methods (SCP, NCP, and CQR) achieve AISL of 0.004±0.001, slightly higher than the performance of the standard deviation interval (0.005±0.000) while maintaining coverage close to the nominal level (89.85-90.02% for CP vs 94.53-95.71% for ensemble spread and standard deviation methods). Notably, CP methods produce slightly narrower intervals (0.002) compared to the standard deviation and ensemble spread approaches (0.004), indicating sharper uncertainty estimates. For h, CP methods again yield comparable AISL values (0.051-0.055) to the standard deviation interval (0.051±0.019), while maintaining similar coverage levels (89.88-89.93% for CP; 88.80% for standard deviation; 87.13% for ensemble spread). For variable r, CP methods achieve AISL values of 0.008, slightly higher than the standard deviation interval (0.007±0.003), but maintain significantly better coverage of 89.81-89.85± than standard deviation and ensemble spread methods with 31% coverage and higher values of miss high (65-66%). The same comparison is evaluated against the true model state (Table 2) and for u, all mean-based CP methods (SCP, NCP,

C.G EORGE ET AL .

Method SCP ensemble based

NCP ensemble based

CQR ensemble based

SCP mean based

NCP mean based

CQR mean based

Std Deviation

Ensemble spread

17

Variable

AISL (mean±std)

Coverage (mean±std)

Interval Width (mean±std)

Miss Low (mean±std)

Miss High (mean±std) 10.51% ± 2.46%

u

0.011 ± 0.001

78.05% ± 4.16%

0.004 ± 0.000

11.44% ± 2.72%

h

0.094 ± 0.026

89.50% ± 3.65%

0.034 ± 0.012

5.60% ± 2.08%

4.89% ± 2.04%

r

0.012 ± 0.004

89.85% ± 3.26%

0.001 ± 0.000

6.67% ± 2.18%

3.48% ± 1.75%

u

0.016 ± 0.001

44.28% ± 4.01%

0.002 ± 0.000

29.82% ± 3.33%

25.90% ± 3.39%

h

0.091 ± 0.014

47.41% ± 4.69%

0.020 ± 0.005

23.37% ± 4.37%

29.22% ± 4.04%

r

0.012 ± 0.004

89.85% ± 3.26%

0.001 ± 0.000

6.67% ± 2.18%

3.48% ± 1.75%

u

0.011 ± 0.001

75.18% ± 4.32%

0.004 ± 0.000

14.08% ± 2.99%

10.74% ± 2.64%

h

0.078 ± 0.024

85.12% ± 4.97%

0.026 ± 0.009

4.49% ± 2.14%

10.39% ± 3.86%

r

0.009 ± 0.004

90.59% ± 3.00%

0.001 ± 0.000

3.94% ± 1.49%

5.47% ± 2.34%

17.55% ± 4.42%

u

0.010 ± 0.002

64.46% ± 6.91%

0.002 ± 0.000

17.99% ± 4.54%

h

0.073 ± 0.024

86.32% ± 4.75%

0.024 ± 0.007

7.99% ± 3.26%

5.69% ± 2.23%

r

0.009 ± 0.004

89.55% ± 3.82%

0.002 ± 0.000

7.76% ± 3.05%

2.69% ± 1.53%

u

0.010 ± 0.001

62.83% ± 6.94%

0.002 ± 0.000

20.18% ± 4.47%

17.00% ± 4.53%

h

0.068 ± 0.020

80.35% ± 6.33%

0.027 ± 0.007

7.22% ± 3.81%

12.44% ± 4.19%

r

0.009 ± 0.004

89.55% ± 3.82%

0.002 ± 0.000

7.76% ± 3.05%

2.69% ± 1.53%

u

0.011 ± 0.002

61.70% ± 6.98%

0.002 ± 0.000

21.49% ± 4.92%

16.81% ± 4.57%

h

0.066 ± 0.022

84.70% ± 5.44%

0.020 ± 0.006

4.51% ± 2.88%

10.79% ± 4.27%

r

0.009 ± 0.003

92.58% ± 2.91%

0.002 ± 0.000

3.97% ± 2.03%

3.45% ± 1.75%

u

0.008 ± 0.001

80.46% ± 4.64%

0.004 ± 0.000

11.17% ± 3.19%

8.38% ± 3.03%

h

0.068 ± 0.020

77.53% ± 5.21%

0.025 ± 0.005

8.48% ± 3.65%

13.99% ± 3.79%

r

0.007 ± 0.003

32.42% ± 8.78%

0.001 ± 0.000

2.56% ± 1.22%

65.02% ± 8.97%

10.89% ± 3.41%

u

-

76.79% ± 4.97%

0.004 ± 0.000

12.32% ± 3.27%

h

-

72.76% ± 5.58%

0.024 ± 0.005

12.43% ± 4.27%

14.81% ± 3.81%

r

-

31.68% ± 8.56%

0.002 ± 0.000

3.66% ± 1.62%

64.65% ± 9.02%

TA B L E 2 Summary of AISL, empirical coverage percentage, interval width and directional error rates percentage (miss-high and miss-low) for the three CP methods (both considering all ensembles and mean based), standard deviation method and ensemble mean spread method averaged over the full assimilation window computed against Truth.

and CQR) produce AISL values around 0.010-0.011, which are slightly higher than the standard deviation interval (0.008). However, CP methods show lower coverage (61-64%) compared to the standard deviation (80.46%) and ensemble spread (76.79%), indicating that the QPEns-based CP interval tends to underestimate the true uncertainty for this variable. This is further reflected in the larger directional miss rates for CP (16-20% miss-high and miss-low) compared to the standard deviation (8-11%) and ensemble spread (10-12%), although CP intervals remain slightly sharper in width (0.002 vs 0.004). For h, the performance of CP methods becomes more comparable to the standard deviation interval and ensemble spread. The AISL values for CP (0.066-0.073) are very similar to the standard deviation (0.068±0.020), with coverage levels of 80-86% for CP compared to 77.53% for the standard deviation and least for ensemble spread with 72.76%. The interval widths are also comparable between methods (0.020-0.027). For variable r, CP methods clearly outperform the standard deviation interval and ensemble spread in terms of coverage. While the standard deviation and ensemble spread approaches shows very poor coverage (31.68-32.42%), all CP methods maintain significantly higher coverage (89-92%). The AISL values are comparable (0.009 for CP vs 0.008 for standard deviation), but the large miss-high rate of the standard deviation and ensemble spread (64-65%) highlights its inability

18

C.G EORGE ET AL .

to capture the asymmetric uncertainty in r. In contrast, CP methods maintain much more balanced and smaller miss rates.

F I G U R E 1 0 Truth (black dashed) and ensemble mean for QPEns (blue) and NN (orange) at timestep 25(left) and 160(right). Gray shading represents uncertainty calculated with mean based SCP Method, std deviation and ensemble spread for the three prognostic variables u (upper), h (middle), r (lower). In Figure 10 for timestep 25 and 160, the SCP mean-based intervals show a consistent balance between sharpness and coverage when compared with traditional methods. For u, the SCP interval closely follows the QPEns mean and remains well aligned with the truth, while maintaining a moderate interval width that is visibly narrower than the standard deviation interval and more stable than the fluctuating ensemble spread. The standard deviation interval tends to be slightly wider and smoother, occasionally overestimating uncertainty, whereas the ensemble spread exhibits larger variability and sporadic spikes. For h, SCP again produces tight and stable intervals that track both QPEns and truth effectively, particularly around localized peak. In contrast, the standard deviation interval appears slightly wider and smoother, missing some sharp variations, while the ensemble spread shows intermittent bursts of large spread, indicating less consistent uncertainty representation. For variable r, where non-Gaussian behaviour is prominent, SCP maintains controlled and adaptive interval widths, capturing key peaks while avoiding excessive widening. The standard deviation interval, although smooth, fails to fully capture these localized extremes, while the ensemble spread shows pronounced spikes and irregular variability.

F I G U R E 1 1 Truth (black dashed) and ensemble mean for QPEns (red) and NN (orange) at timestep 25 (left) and 160 (right). Gray shading represents uncertainty calculated with mean based NCP Method, std deviation and ensemble spread for the three prognostic variables u (upper), h (middle), r (lower).

C.G EORGE ET AL .

19

Now in Figure 11, the NCP (mean-based) intervals exhibit a sharper but more conservative behaviour compared to SCP, particularly in how the interval width adapts relative to QPEns and the truth. For u, NCP produces consistently narrower and smoother intervals than SCP, closely following the QPEns mean but often appearing slightly under covered with respect to the truth, especially in regions with stronger deviations. In contrast, SCP shows slightly wider and more adaptive intervals that better accommodate local variability. The standard deviation interval remains wider than both CP methods, while the ensemble spread shows large fluctuations and instability. For h, NCP maintains tight and structured intervals that capture the main peaks but are generally narrower than SCP, indicating reduced spread around sharp gradients. The standard deviation interval again appears moderately wider, whereas the ensemble spread exhibits high variability. For variable r, since the same SCP-based formulation is used, both methods produce identical interval structures, showing controlled and adaptive interval widths and the standard deviation interval remains relatively rigid and unable to fully adapt to localized extremes, while the ensemble spread shows strong spikes and irregular variability. In Figure 12, the CQR (mean-based) intervals demonstrate a more adaptive and asymmetric

F I G U R E 1 2 Truth (black dashed) and ensemble mean of QPEns (blue),upper quantile (1 − α/2=0.95,cyan), lower quantile (α/2=0.05,green) and midpoint of quantiles (orange) at timestep 25 (left) and 160 (right). Gray shading represents uncertainty calculated with mean based CQR Method, std eviation and ensemble spread for the three prognostic variables u (upper), h (middle), r (lower).

representation of uncertainty compared to SCP and NCP, particularly in how the quantile bounds respond to variations relative to QPEns and the truth. For u, CQR produces intervals that are slightly wider than NCP, allowing better representation of local deviations while still remaining smoother than the highly fluctuating ensemble spread. Compared to SCP, CQR shows comparable alignment with QPEns but with more adaptive upper and lower bounds, capturing asymmetry in deviations from the truth. The standard deviation interval remains more symmetric than ensemble spread which continues to show more variability. For h, CQR effectively captures sharp localized peaks at both timesteps, with the quantile bounds expanding appropriately around strong gradients. The standard deviation interval appears smoother and slightly less flexible to these peaks, while the ensemble spread shows irregular bursts in high variability regions. For variable r, CQR clearly highlights its strength in handling non-Gaussian and skewed behaviour. The quantiles expand asymmetrically around sharp spikes, capturing the structure of the truth more effectively than both SCP and NCP. In contrast, the standard deviation interval remains relatively rigid and fails to adapt to these localized extremes similar to the ensemble spread.

20

5.5

C.G EORGE ET AL .

|

Comparison of EnKF and CNN with and without CP Perturbations

Two configurations are considered for incorporating CP based uncertainty estimates into the data assimilation cycle. In the first configuration, CP perturbations are applied to the CNN analysis, while in the second configuration, CP perturbations are applied to the EnKF analysis (see Figure 3). For each configuration, we evaluate two cases: EnKF and CNN analyses without CP perturbations, and the corresponding EnKF and CNN analyses with CP perturbations applied at the analysis or output stage. For each CP method (SCP, NCP, and CQR) computed based on ensembles, the EnKF and CNN results without CP perturbations are identical since these CP values are not used and any differences between methods therefore arise solely from the inclusion of CP perturbations. RMSE is used to compare the results of each method. RMSE is computed as the mean over 200 experiments for each method. During the first 20 assimilation cycles, all configurations exhibit nearly identical RMSE values, corresponding to the spin-up phase of the data assimilation system.

F I G U R E 1 3 RMSE analysis for the variables u (left), h (middle), r (right) of SCP (a), NCP (b) and CQR (c) methods and comparing EnKF and CNN analysis with and without CNN-CP perturbations.

Figure 13 shows the RMSE evolution for the background and analysis states when CP perturbations are introduced at the CNN output before the EnKF update from timestep 21. In this configuration, the introduction of CP uncertainty consistently improves the EnKF performance compared to the no-CP baseline. This improvement is reflected in both the background and

C.G EORGE ET AL .

21

analysis RMSE, indicating that the perturbations introduced through CP propagate through the assimilation cycle and help maintain a more informative ensemble representation. The background of CP-influenced EnKF and the CP-influenced NN produce are same. This indicates that the background state used by the EnKF directly reflects the CP-perturbed CNN output. For all the CP methods, the reduction in RMSE is consistent across all three variables. In both methods, the EnKF with CP perturbations maintain lower RMSE values compared to the no-CP configuration for both background and analysis estimates. In the CQR method, the EnKF continues to show lower RMSE than the no-CP configuration across all variables, but higher in comparison to SCP and NCP methods especially for u and r. The same bias is observed in the SCP method for u but not in NCP which points us to the lack of potential in this setup to capture the stochastic perturbations added to u at each timestep.

F I G U R E 1 4 RMSE analysis for the variables u (left), h (middle), r (right) of SCP (a), NCP (b) and CQR (c) methods and comparing EnKF and CNN analysis with and without EnKF-CP perturbations.

In the second configuration, CP perturbations are introduced directly into the EnKF analysis rather than the CNN output. Figure 14 shows the corresponding background and analysis RMSE evolution for all variables. In contrast to the first configuration, the CNN output is not used to update the EnKF state when CP perturbations are applied. Instead, the conformal perturbations are added directly to the EnKF analysis, while the CNN output acts only as a post-processing step. As a result, the influence of CP on the EnKF and NN behaviour differs from the previous configuration. For all the CP methods, introducing CP perturbations into the EnKF analysis does not lead to noticeable improvements in RMSE compared to the no-CP EnKF baseline. The background

22

C.G EORGE ET AL .

and analysis RMSE for the EnKF remain very similar to the no-CP configuration for variable h but worse for variable u in SCP and CQR and also for variable r in CQR, indicating that these perturbations do not significantly alter the assimilation performance. However, the NN results in this configuration show slightly higher RMSE compared to the NN without CP. This behaviour arises because the CNN output is not used to influence the EnKF update; instead, it is applied only as a post-processing step. Consequently, the CP perturbations introduced in the EnKF analysis do not improve the NN estimates and can lead to a small degradation in NN performance. For the CQR method, a different behaviour can be observed. In this case, the perturbations lead to a slight reduction in EnKF analysis RMSE across all variables compared to the no-CP baseline. Although the overall improvements remain decent, the reduction in RMSE is more stable and visible than for SCP and CQR, indicating that NCP interacts more effectively with the EnKF update when perturbations are introduced directly into the analysis.

6

|

CONCLUDING REMARKS

In this study, we investigated the integration of CP based uncertainty quantification into a hybrid ML-data assimilation framework. Specifically, we evaluated three variants of the CP-SCP, NCP, and CQR. We analyzed their performance through empirical coverage, interval width, AISL scores, and spatial snapshots. The analysis first focused on the intrinsic behaviour of CP methods, including a comparison between ensemble-based and mean-based formulations, with traditional methods and subsequently on their impact when incorporated into the data assimilation cycle under different perturbation configurations. It is also worth noting that, while CP requires exchangeability of the data, a condition not strictly satisfied by our ensemble dataset, the results confirm that CP remains practically valid in this dynamical system setting, producing coverage consistently close to the nominal level. All CP methods achieve coverage close to the nominal 90% level when evaluated against the reference distribution, confirming their validity for uncertainty quantification. However, clear differences emerge in terms of stability and adaptivity. The comparison between ensemble-based and mean-based formulations reveals that the optimal representation depends on the CP method: SCP and CQR benefit from using the ensemble mean, leading to tighter and more stable intervals, whereas NCP performs better when applied directly to the ensemble members, as its normalization inherently captures local variability. Across methods, NCP provides the best overall balance between sharpness and reliability, while CQR demonstrates clear advantages in handling non-Gaussian and asymmetric variables, particularly for rainfall. And SCP produces more uniform and less adaptive intervals due to its global conformity threshold. When compared with traditional uncertainty quantification approaches, namely the standard deviation interval and ensemble spread, CP methods demonstrate clear advantages. While the standard deviation interval assumes Gaussianity and provides fixed symmetric bounds, CP methods produce adaptive intervals that maintain sharpness while achieving more consistent coverage, especially for non-Gaussian variables. The ensemble spread, lacking a probabilistic guarantee, exhibits high variability and cannot provide calibrated uncertainty estimates. Notably, CP methods significantly outperform traditional approaches for skewed variables such as rainfall, where Gaussian assumptions break down. We further evaluated two perturbation configurations that differ in the placement of CP-based uncertainty within the hybrid framework: perturbations applied at the CNN output and at the EnKF analysis. Introducing CP perturbations at the CNN output consistently improves the EnKF analysis across all methods, with the strongest gains observed for NCP, followed by CQR and SCP. So, applying CP directly at the EnKF analysis stage results in limited or no improvement, indicating that modifying the analysis state after the update is less effective than improving the prior ensemble representation. Also, the CNN itself does not improve under CP perturbations, as it is not retrained to incorporate conformal uncertainty. Overall, the results highlight that both the uncertainty structure (choice of the CP method) and its placement within the assimilation pipeline are critical to enhancing hybrid ML data assimilation systems. Among the approaches tested, the combination of NCP with CNN output perturbation yields the most consistent improvements. Future work will focus on improving the inclusion of CP into data assimilation cycling by moving from the current offline

C.G EORGE ET AL .

23

computation of CP conformity scores to an online framework, where prediction intervals are updated sequentially during data assimilation cycle. Also, we focus on more tightly integrating CP within the data assimilation cycle by retraining the CNN online so that uncertainty estimates evolve consistently with the dynamically updated model state, rather than being applied post hoc. In addition to this, we will focus on analyzing whether incorporating aleatoric and epistemic uncertainty can improve the intervals, producing insights for more reliable and interpretable atmospheric forecasting. By separating inherent data variability from model-related uncertainty, we aim to develop more physically consistent ensemble representations and further improve the performance and reliability of hybrid ML-data assimilation systems.

D ATA A N D C O D E The source code (https://doi.org/10.5281/zenodo.4354602, Ruckstuhl et al., 2020) includes the necessary scripts to produce the training data and CNN and https://github.com/Catherinegeo98/Uncertainty-quantification-via-conformal-prediction-in-dataassimilation contains the necessary codes for CP anaylsis and perturbing CP estimates into data assimilation.

AC K N OW L E D G E M E N T S The authors are grateful for the financial support of the Klaus Tschira Stiftung gGmbH, whose funding made this project possible. The authors also wish to thank Yvonne Ruckstuhl for her valuable insights and for the foundational work established in her prior research, upon which this study builds.

C O N FL I C T O F I N T E R E S T S TAT E M E N T The authors declare that they have no conflict of interest.

REFERENCES Anastasios N. Angelopoulos and Stephen Bates. A gentle introduction to conformal prediction and distribution-free uncertainty quantification. 2022. URL https://arxiv.org/abs/2107.07511. Imanol Arrieta-Ibarra, Paman Gujral, Jonathan Tannen, Mark Tygert, and Cherie Xu. Metrics of calibration for probabilistic predictions. Journal of Machine Learning Research, 23(351):1–54, 2022. URL http://jmlr.org/papers/v23/22- 0658.html. L. Breiman. Random forests. Machine Learning, 45(1):5–32, 2001. Gerrit Burgers, Peter Jan van Leeuwen, and Geir Evensen. Analysis scheme in the ensemble Kalman filter. Mon. Wea. Rev, 126: 1719–1724, 1998. Ashesh Chattopadhyay, Ebrahim Nabizadeh, Eviatar Bach, and Pedram Hassanzadeh. Deep learning-enhanced ensemble-based data assimilation for high-dimensional nonlinear dynamical systems. Journal of Computational Physics, 477:111918, 2023. ISSN 0021-9991. doi: https://doi.org/10.1016/j.jcp.2023.111918. URL https://www.sciencedirect.com/science/article/ pii/S002199912300013X. G. Evensen. Sequential data assimilation with a nonliear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. Journal of Geophysical Research, 99:10143–10162, 1994. G. Evensen. Data Assimilation: The Ensemble Kalman Filter. Springer, 2009.

24

C.G EORGE ET AL .

Geir Evensen. The ensemble kalman filter: theoretical formulation and practical implementation. Ocean Dynamics, 53(4):343–367, 11 2003. ISSN 1616-7228. doi: 10.1007/s10236- 003- 0036- 9. URL https://doi.org/10.1007/s10236- 003- 0036- 9. Y. Gal and Z. Ghahramani. Dropout as a Bayesian approximation: Representing model uncertainty in deep learning. In Proc. ICML, 33rd Int. Conf. on Machine Learning, New York, NY, USA, 2016. Tabea Gleiter, Tijana Janjić, and Nan Chen. Ensemble kalman filter based data assimilation for tropical waves in the MJO skeleton model. Quarterly Journal of the Royal Meteorological Society, 148(743):1035–1056, 2022. doi: https://doi.org/10.1002/qj.4245. URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.4245. Tilmann Gneiting and Adrian E. Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 102:359 – 378, 2007. URL https://api.semanticscholar.org/CorpusID:1878582. Ethan Goan and Clinton Fookes. Bayesian Neural Networks: An Introduction and Survey, page 45–87. Springer International Publishing, 2020. ISBN 9783030425531. doi: 10.1007/978- 3- 030- 42553- 1_3. URL http://dx.doi.org/10.1007/978- 3- 030- 42553- 1_3. Ian Goodfellow, Yoshua Bengio, and Aaron Courville. Deep Learning. MIT Press, 2016. http://www.deeplearningbook.org. Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q. Weinberger. On calibration of modern neural networks. In Doina Precup and Yee Whye Teh, editors, Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 1321–1330. PMLR, 06–11 Aug 2017. URL https://proceedings.mlr.press/v70/guo17a.html. M. Haslehner, Tijana Janjić, and George C. Craig. Testing particle filters on simple convective-scale models. Part 2: A modified shallow-water model. Q J R Meteorol Soc, 142(697):1628–1646, 2016. ISSN 1477-870X. doi: 10.1002/qj.2757. T. Janjić, Y. Ruckstuhl, and P. L. Toint. A data assimilation algorithm for predicting rain. Q J R Meteorol Soc, 147(736):1949–1963, 2021. doi: https://doi.org/10.1002/qj.4004. URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.4004. Tijana Janjić, Dennis McLaughlin, Stephen E. Cohn, and Martin Verlaan. Conservation of mass and preservation of positivity with ensemble-type kalman filter algorithms. Monthly Weather Review, 142(2):755 – 773, 2014. doi: 10.1175/MWR-D- 13- 00056.1. Tijana Janjić, Maria Lukáčová-Medvid’ová, Yvonne Ruckstuhl, and Bettina Wiebe. Comparison of uncertainty quantification methods for cloud simulation. Quarterly Journal of the Royal Meteorological Society, 149(756):2895–2910, 2023. Alireza Javanmardi and Eyke Hüllermeier. Conformal prediction intervals for remaining useful lifetime estimation. International Journal of Prognostics and Health Management, (2), 2023. Eugenia Kalnay. Atmospheric Modeling, Data Assimilation and Predictability. Cambridge University Press, 2002. B. Lakshminarayanan, A. Pritzel, and C. Blundell. Simple and scalable predictive uncertainty estimation using deep ensembles. In Proc. NeurIPS, 31st Conf. on Neural Information Processing Systems, Long Beach, California, USA, 2017. Yann LeCun and Geoffrey Hinton. Deep learning. Nature, 521:436–44, 05 2015. doi: 10.1038/nature14539. S. Legler and T. Janjić. Combining data assimilation and machine learning to estimate parameters of a convective-scale model. Quarterly Journal of the Royal Meteorological Society, 148(743):860–874, 2022. doi: https://doi.org/10.1002/qj.4235. URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.4235. M. Leutbecher and T.N. Palmer. Ensemble forecasting. Journal of Computational Physics, 227(7):3515–3539, 2008. ISSN 00219991. doi: https://doi.org/10.1016/j.jcp.2007.02.014. URL https://www.sciencedirect.com/science/article/pii/ S0021999107000812. Linus Magnusson, Erland Källén, and Jonas Nycander. Initial state perturbations in ensemble forecasting. Nonlinear processes in geophysics, 15(5):751–759, 2008. URL http://www.nonlin-processes-geophys.net/15/751/2008/npg- 15- 751- 2008.html. R.N. Neal. Bayesian Learning for Neural Networks. Springer-Verlag, Berlin, Heidelberg, 1996.

C.G EORGE ET AL .

25

Harris Papadopoulos, Alex Gammerman, and Volodya Vovk. Normalized nonconformity measures for regression conformal prediction. In Proceedings of the IASTED International Conference on Artificial Intelligence and Applications (AIA 2008), 2008. Ilan Price, Alvaro Sanchez-Gonzalez, Ferran Alet, Tom R. Andersson, Andrew El-Kadi, Dominic Masters, Timo Ewalds, Jacklynn Stott, Shakir Mohamed, Peter Battaglia, Remi Lam, and Matthew Willson. Probabilistic weather forecasting with machine learning. Nature, 637(8044):84–90, 2025. ISSN 1476-4687. doi: 10.1038/s41586- 024- 08252- 9. URL https://doi.org/10.1038/s41586- 02408252- 9. Stephan Rasp, Michael S. Pritchard, and Pierre Gentine. Deep learning to represent subgrid processes in climate models. Proceedings of the National Academy of Sciences, 115(39):9684–9689, 2018. doi: 10.1073/pnas.1810286115. URL https://www.pnas.org/ doi/abs/10.1073/pnas.1810286115. Markus Reichstein, Gustau Camps-Valls, Bjorn Stevens, Martin Jung, Joachim Denzler, Nuno Carvalhais, and Mr Prabhat. Deep learning and process understanding for data-driven earth system science. Nature, 566:195, 02 2019. doi: 10.1038/s41586- 019- 0912- 1. Yaniv Romano, Evan Patterson, and Emmanuel Candes. Conformalized quantile regression. Advances in neural information processing systems, 32, 2019. Y. Ruckstuhl and T. Janjić. Parameter and state estimation with ensemble Kalman filter based approaches for convective scale data assimilation. Q J R Meteorol Soc, 144(712):826–841, 2018. doi: https://doi.org/10.1002/qj.3257. Y. Ruckstuhl, T. Janjić, and S. Rasp. Training a convolutional neural network to conserve mass in data assimilation. Nonlinear Processes in Geophysics, 28(1):111–119, 2021. doi: 10.5194/npg- 28- 111- 2021. URL https://npg.copernicus.org/articles/28/111/ 2021/. Vladimir Vovk, Alexander Gammerman, and Glenn Shafer. Algorithmic Learning in a Random World. Springer Nature, 2022. M. Würsch and G. C Craig. A simple dynamical model of cumulus convection for data assimilation research. Meteorologische Zeitschrift, 23:483–490, 2014. Yuefei Zeng, Tijana Janjić, Yvonne Ruckstuhl, and Martin Verlaan. Ensemble-type kalman filter algorithm conserving mass, total energy and enstrophy. Quarterly Journal of the Royal Meteorological Society, 143(708):2902–2914, 2017. doi: https://doi.org/10.1002/ qj.3142. URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.3142.