Blanka Balogh
email: blanka.balogh@meteo.fr

D Saint-Martin

A Ribes

How to Calibrate a Dynamical System With Neural Network Based Physics?

. Recently, interest has grown in using machine learning (ML) to develop parameterizations, noted 𝐴𝐴 f𝑓 , promising more precise yet affordable parameterizations for climate models (e.g.,

Introduction

The ordinary differential equations describing a climate model, can be written ẋ𝐱 = (𝐱𝐱) + 𝜙𝜙(𝐱𝐱), (1) BALOGH ET AL.

10.1029/2022GL097872

2 of 9 evaluates the error between the estimate function (i.e., f ϕ or 𝐴𝐴 f𝑓 ) and a target function f, the best accurate approximation of ϕ-as even high-resolution model data is only an approximation of the ground truth subgrid-scale processes. This approach allows us to fit the models on a point by point basis. However, a climate model's performance lies in accurate predictions of long-term statistics, measured by online metrics. Physical parameterizations typically require an additional calibration step, once it has been plugged-in to the dynamical model. In contrast, there is usually no tunable parameter included in ML parameterizations. Nevertheless, tuning is very important for the development of climate models, to ensure model stability, to calibrate the value of long-term statistics and to reduce error due to potentially missed interactions between the dynamics and the physics. This online calibration can be achieved regarding an online metric, m. Without online calibration, even though a parameterization is efficient regarding the offline metric ℓ, there is no guarantee that the climate model will be accurate regarding the online evaluation metric [START_REF] Brenowitz | Machine learning climate model dynamics: Offline versus online performance[END_REF].

The development of current ML parameterizations is still hampered by issues addressed by online calibration in the case of physical parameterizations, such as numerical instabilities (e.g., [START_REF] Brenowitz | Interpreting and stabilizing machine-learning parametrizations of convection[END_REF][START_REF] Rasp | Coupled online learning as a way to tackle instabilities and biases in neural network parameterizations: General algorithms and Lorenz 96 case study (v1.0)[END_REF]. These issues are currently handled by enforcing physical conservation laws [START_REF] Beucler | Enforcing analytic constraints in neural networks emulating physical systems[END_REF] or using ML to replace part of a physical parameterization [START_REF] Yuval | Use of neural networks for stable, accurate and physically consistent parameterization of subgrid atmospheric processes with good performance at reduced precision[END_REF]. Additionally, the target function f can also be imperfect. In this case, even though the fit of the ML model is excellent, errors in the imperfect training data set will also be learned and result in a high online error m, also called long-term biases. This is a well known issue with physical parameterizations. The "art of tuning" model parameters consists in finding a compromise between offline and online metrics, going back and forth in optimizing the parameters with respect to ℓ and m [START_REF] Couvreux | Process-based climate model development harnessing machine learning: I. A calibration tool for parameterization improvement[END_REF][START_REF] Hourdin | The art and science of climate model tuning[END_REF][START_REF] Schmidt | Practice and philosophy of climate model tuning across six US modeling centers[END_REF].

To resolve these issues, we propose to include some tunable parameters into the ML parameterization inputs. This allows the calibration of new ML-based parameterization with respect to an online metric in an efficient way. Our approach to find the best value of tunable parameters is less empirical than the methods used to tune physical parameterizations. It relies on the minimization of the online metric m on long-term statistics, akin to the methodology described in [START_REF] Schneider | Earth system modeling 2.0: A blueprint for models that learn from observations and targeted high-resolution simulations: Earth system modeling 2.0[END_REF] and in [START_REF] Cleary | Calibrate, emulate, sample[END_REF]. The method will be demonstrated using Lorenz'63 and Lorenz'96 toy models [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF][START_REF] Lorenz | Predictability: A problem partly solved[END_REF]. These toy models offer a simple framework to perform proof-of-concept atmospheric modeling experiments using ML tools (e.g., [START_REF] Chattopadhyay | Data-driven predictions of a multiscale Lorenz 96 chaotic system using machine-learning methods: Reservoir computing, artificial neural network, and long short-term memory network[END_REF][START_REF] Scher | Generalization properties of feed-forward neural networks trained on Lorenz systems[END_REF]. Lorenz'96 model has also the advantage of implementing a "parameterization." The large scale variable is indeed interacting with the non-linear subgrid-scale variable.

The paper is organized as follows. Methodology is described in Section 2. Section 3 demonstrates our calibration method on two basic examples using Lorenz'63 and Lorenz'96 models. Conclusions are drawn in Section 4, discussing more broadly the results we have obtained.

Methodology

Step 1: Building a Tunable Neural Network Parameterization

The first step toward designing and optimizing tunable neural network (NN) parameterizations is to take into consideration some uncertain but tunable parameters, when fitting the NN. The target function f is often fitted to outputs from high-resolution simulations. This target function depends on uncertain parameters, noted 𝐴𝐴 𝜽𝜽 ∈ ℝ 𝑝𝑝 , used to generate the high-resolution simulations: f ≡ f(x; θ). In traditional approaches, the value of model parameters are fixed to the "best estimate" value of these parameters, θ 0 . Thus, f ≡ f(x; θ 0 ) does not depend on any parameter and θ is not included in the NN model input (e.g., [START_REF] Gentine | Could machine learning break the convection parameterization deadlock[END_REF]Yuval & O'Gorman, 2020). The novelty of our approach is to keep part of the uncertain parameters θ among the input variables of the NN, and thus to retain not only the dependency of f from x but also from some model parameters θ.

The NN learning sample, of size N, is 𝐴𝐴 {(𝐱𝐱, 𝜽𝜽)𝑖𝑖, 𝑓𝑓 (𝐱𝐱𝑖𝑖; 𝜽𝜽𝑖𝑖)} 1≤𝑖𝑖≤𝑁𝑁 . The NN is fitted by optimizing an offline metric ℓ, also called "loss function." Typically, this loss function is Mean Squared Error (MSE), computed for each predicted time step individually:

( ; ), ̂ ( ; ) = 1 ∑ =1 ‖ ‖ ‖ ‖ ( ; ) -̂ ( ; ) ‖ ‖ ‖ ‖ 2 . ( 2
)
BALOGH ET AL.

10.1029/2022GL097872

3 of 9

The parameterization obtained after training the NN is noted, 𝐴𝐴 f𝑓 (𝐱𝐱; 𝜽𝜽) . It can be incorporated into the dynamical model so as to replace ϕ in Equation 1. The resulting dynamical system can be used to generate a validation time series of the dynamical model, noted 𝐴𝐴 [𝐱𝐱] f𝑓 .

Step 2: Optimizing the Tunable Neural Network Parameterization

The goal of this second step is to tune the optimal value of the parameters θ. Whereas the NN model was trained offline on a learning sample, calibration relies on an online validation metric, m, given by: Theoretically, the statistics 𝐴𝐴  are computed over a time series of infinite length. In reality, we only compute an estimate of the long-term metric over a finite length simulation. The length of these simulation is chosen so that the online metric m no longer depends on the x initial condition. To simplify notations, the estimate of the longterm statistics over validation time series 𝐴𝐴 [𝐱𝐱] is noted, 𝐴𝐴  . Even though this time series is sufficiently long, the

𝑚𝑚(𝜽𝜽) = 1 𝑀𝑀 𝑀𝑀 ∑ 𝑘𝑘=1 ‖ ref 𝑘𝑘 -𝑘𝑘(𝜽𝜽)‖ 2 , ( 3
)
resulting metric ̂ ( ) = 1 ∑ =1 ‖ ‖ ‖  ref - ( ) ‖ ‖ ‖ 2
will be noised, which can lead the optimizer to a local minimum of 𝐴𝐴 𝐴 𝐴𝐴 , instead of the global minimum. To address this issue, minimization will be performed over a smoothed version of 𝐴𝐴 𝐴 𝐴𝐴 , obtained by kriging (or Gaussian Process Regression, [START_REF] Cressie | Statistics for spatial data[END_REF]. The kriging metamodel will be fitted to a sample of

𝐴𝐴 { 𝜽𝜽𝑖𝑖, m𝑚 (𝜽𝜽𝑖𝑖) } .
The obtained kriging metamodel is noted 𝐴𝐴 𝐴 𝐴𝐴 . The optimal value of model parameters, θ*, is obtained when 𝐴𝐴 𝐴 𝐴𝐴 reaches its (absolute) minimum value.

In summary, we train a NN parameterization depending on some tunable parameters θ, using an offline metric, ℓ. Since the main purpose is to reduce long-term prediction errors of the dynamical model, the value of θ is subsequently optimized regarding an online metric, m. In practice, optimization is done over a kriging metamodel emulating m as a function of θ. In the following, our method will be demonstrated using the Lorenz'63 (hereafter L63, Lorenz, 1963) and the Lorenz'96 (hereafter L96, Lorenz, 1996) models.

The Lorenz'63 Model

The L63 system consists of a set of ordinary differential equations that can be expressed:

̇ 1 = ( 2 -1) , ̇ 2 = 1 ( -3) -2, ̇ 3 = 1 2 -3. (4)
Temporal evolution of the L63 state variable, x = (x 1 , x 2 , x 3 ), is governed by Equation 4, which involves a set of three model parameters, (σ, ρ, β). Equation 4 admits a chaotic solution in the vicinity of (σ 0 , ρ 0 , β 0 ) = (10, 28, 8/3).

The online minimization metric m (Equation 3) depends strongly on the identification of variables quantifying the long-term statistical behavior of the L63 system. The first candidate of such a long-term variable is the mean value for each one of the three L63 state variables. However, the average value of x 1 and x 2 is independent of L63 parameters (see also Figure S1 in Supporting Information S1). Thus, only the average value of x 3 , noted μ (3) , can be tuned. This quantity is useful to retrieve the optimal value of at most one L63 parameter. In addition to μ (3) , standard deviations over the time series can also be computed. Standard deviations will be noted σ (i) for each one of the L63 parameters, with i ∈ {1, 2, 3}. However, as highlighted by Figure S1 in Supporting Information S1, a strong correlation appears between the three standard deviations. Hence, it is possible to retrieve the optimal value of at most two L63 parameters by optimization. Thus, the following L63 idealized case study will use 𝐴𝐴  = ( 𝜇𝜇 (3) , 𝜎𝜎 (1) , 𝜎𝜎 (2) , 𝜎𝜎 (3) ) .

BALOGH ET AL.

10.1029/2022GL097872 4 of 9

The Lorenz'96 Model

The 2-level L96 dynamical system depending on parameters (h, c, F, b) is given by:

d𝑥𝑥𝑘𝑘 d𝑡𝑡 = -𝑥𝑥𝑘𝑘-1 (𝑥𝑥𝑘𝑘-2 -𝑥𝑥𝑘𝑘+1) -𝑥𝑥𝑘𝑘 + 𝐹𝐹 - ℎ𝑐𝑐 𝑏𝑏 𝐽𝐽 ∑ 𝑗𝑗=1 𝑦𝑦𝑘𝑘𝑘𝑗𝑗𝑘 (5) 1 𝑐𝑐 d𝑦𝑦𝑘𝑘𝑘𝑘𝑘 d𝑡𝑡 = -𝑏𝑏𝑦𝑦𝑘𝑘𝑘𝑘𝑘+1 (𝑦𝑦𝑘𝑘𝑘𝑘𝑘+2 -𝑦𝑦𝑘𝑘𝑘𝑘𝑘-1) -𝑦𝑦𝑘𝑘𝑘𝑘𝑘 + ℎ𝑐𝑐 𝑏𝑏 𝑥𝑥𝑘𝑘𝑘 (6)
where

𝐴𝐴 {𝑥𝑥𝑘𝑘} 1≤𝑘𝑘≤𝐾𝐾 are K large-scale variables coupled to J small-scale variables 𝐴𝐴 { 𝑦𝑦𝑘𝑘 0 ,𝑗𝑗 } 1≤𝑗𝑗≤𝐽𝐽 for each 1 ≤ k 0 ≤ K. The coupling term 𝐴𝐴 𝐴𝐴𝑘𝑘 = -ℎ𝑐𝑐 𝑏𝑏 ∑ 𝐾𝐾 𝑘𝑘=1
𝑦𝑦𝑘𝑘𝑘𝑘𝑘 can be seen as a "subgrid-scale parameterization." Thus, Equation 5 can be analogous to a (very) simple climate model, as described by Equation 1. In this comparison, 𝐴𝐴  (𝑥𝑥𝑘𝑘) = -𝑥𝑥𝑘𝑘-1 (𝑥𝑥𝑘𝑘-2 -𝑥𝑥𝑘𝑘+1) -𝑥𝑥𝑘𝑘 + 𝐹𝐹 and ϕ(x k ) = B k . B k can be approximated as a function of the large-scale state variables, x k , and is often modeled with polynoms (e.g., [START_REF] Arnold | Stochastic parametrizations and model uncertainty in the Lorenz '96 system[END_REF]. This subgrid-scale parameterization will be the target variable of the NN model. Given the symmetry of the L96 model, it is common to fit the NN to predict B k as a function of data from only one spatial variable (e.g., [START_REF] Gagne | Machine learning for stochastic parameterization: Generative adversarial networks in the Lorenz '96 model[END_REF][START_REF] Rasp | Coupled online learning as a way to tackle instabilities and biases in neural network parameterizations: General algorithms and Lorenz 96 case study (v1.0)[END_REF][START_REF] Watson | Applying machine learning to improve simulations of a chaotic dynamical system using empirical error correction[END_REF]. In our case, the NN learns to predict a single B k as a function of (x k , c). This implies that the NN input is of size 2 and the output of size 1. To simplify notations, in the following, the NN will be noted, 𝐴𝐴 B𝐵(𝑥𝑥𝑥 𝑥𝑥) .

As in the L63 case, the mean value of the state variables does not depend on input parameters (h, c, F, b), and is therefore not used in the online metric m. Thus, the metric is based on standard deviations only, estimated on time series obtained with the NN parameterization.

Case Studies

Perfect Model Calibration

The Lorenz'63 Model

The objective is to fit an NN model, noted 𝐴𝐴 f𝑓 , to approximate the L63 time derivative as a function of the state variable x and parameters θ = (ρ, β):

x𝐱 = f𝑓 (𝐱𝐱, 𝜽𝜽). ( 7
)
The learning sample is generated by an optimal sampling method, used to select relevant 𝐴𝐴 (𝐱𝐱, 𝜽𝜽) ∈ ℝ 3 × ℝ 2 values to build the NN learning sample. Latin Hypercube Sampling (hereafter LHS, [START_REF] Mckay | Latin Hypercube Sampling as a tool in uncertainty analysis of computer models[END_REF]) is such a sampling method. The interest of using a specific sampling method to train stable and accurate NN based parameterization has been shown in [START_REF] Balogh | A toy model to investigate stability of Ai-based dynamical systems[END_REF]. The boundaries of the LHS are set around plausible values for both x and θ. The resulting learning sample of size

N LHS is 𝐴𝐴 [(𝐱𝐱, 𝜽𝜽)] LHS = {(𝐱𝐱, 𝛽𝛽)𝑖𝑖, 𝑓𝑓 (𝐱𝐱𝑖𝑖; 𝜽𝜽𝑖𝑖)} 1≤𝑖𝑖≤𝑁𝑁 LHS .
The target variable 𝐴𝐴 ẋ𝐱 is computed with L63 model equations, parameterized with σ = σ 0 . The learning sample is of size N LHS = 10 7 . Values for θ = (ρ, β) are sampled in 𝐴𝐴 [26.5, 32] × [1.5, 3.2] . The NN model consists of n l = 7 hidden layers of type "Dense." More specific details about the NN architecture are available in the Table S1 in Supporting Information S1. R 2 score over an independent subset of 20% of the learning sample is monitored during training. The best weights regarding the R 2 score are loaded after 30 epochs. The final model has R 2 = 1.00, which is not surprising given the low complexity of L63 model, and is consistent with case studies focusing on emulating L63 dynamics with NNs [START_REF] Rasp | Coupled online learning as a way to tackle instabilities and biases in neural network parameterizations: General algorithms and Lorenz 96 case study (v1.0)[END_REF].

The fitted NN model is then used to generate time derivatives, which are integrated with a Runge-Kutta 4 time stepping scheme with a temporal increment Δt = 0.05. The resulting validation time series or "orbit" 𝐴𝐴

[𝐱𝐱] f𝑓 is of length 1000 Model Time Units (hereafter, MTU, where 1 MTU = 20Δt) not including a spin-up of 200 MTU.

The validation orbits are used to compute long-term metrics and thus to assess the online performance of the model, measured by 𝐴𝐴 𝐴 𝐴𝐴 :

m𝑚(𝜽𝜽) = ( μ𝜇 (3) (𝜽𝜽) -𝜇𝜇 (3) ref ) 2 + 3 ∑ 𝑛𝑛=1 ( σ𝜎 (𝑛𝑛) (𝜽𝜽) -𝜎𝜎 (𝑛𝑛) ref ) 2 , ( 8
)
where f𝑓 ). In this example, the "reference" data set consist in a long time-series of length 3000 MTU (considered as "infinite" lenght), generated by integrating L63 equations (Equation 4) with parameters (σ 0 , ρ 0 , β 0 ).

A kriging metamodel learns to approximate 𝐴𝐴 𝐴 𝐴𝐴 as a function of θ, over a learning sample of size N k = 750. To build the kriging learning data set, θ values are sampled in the interval 𝐴𝐴 [26.5, 32] × [1.5, 3] . The sampling interval for θ has been slightly reduced compared with those used to generate the NN learning sample to avoid potential out-of-sample issues. In practice, we use a LHS to generate a sample of (x 0 , θ) to compute the validation orbits of length 1000 MTU on which the long-term metrics are estimated. Including the initial state variable x 0 in the LHS sample reduces the noise related to the (finite) length of the validation orbits from the kriging learning sample. For the reasons explained in the Methodology section, dependency of

𝐴𝐴 𝐴 𝐴𝐴 on x 0 is ignored. The kriging model 𝐴𝐴 𝐴 𝐴𝐴 is fitted to 𝐴𝐴 { 𝜽𝜽𝑖𝑖, m𝑚 (𝜽𝜽𝑖𝑖) } 1≤𝑖𝑖≤𝑁𝑁 𝑘𝑘
. The minimum of 𝐴𝐴 𝐴 𝐴𝐴 is retrieved with BFGS optimizer [START_REF] Fletcher | Practical methods of optimization[END_REF]. The optimal value found for the tunable parameters is θ* = (27.9, 2.64) (Figure 1). The reference orbit was generated using (ρ 0 , β 0 ) = (28, 8/3): θ* is close to the values of the parameters used to generate the reference orbit.

The Lorenz'96 Model

The NN model is trained to predict 𝐴𝐴 𝐴𝐴= -ℎ𝑐𝑐 𝑏𝑏 ∑ 𝐽𝐽 𝑗𝑗=1 𝑦𝑦𝑗𝑗 as a function of x and θ = c. The learning sample is built by sequentially integrating the L96 equations (see Equations 5 and 6, with K = 8 and J = 32), using a fourth order Runge-Kutta time stepping scheme with an increment Δt = 0.005. The length of training integrations is 3,5 MTU (where 1 MTU = 200 Δt), not including 1,5 MTU of model spin-up. We perform N i = 500 integrations, the initial conditions and θ values (x, y, c) of which are sampled using LHS, where x = (x 1 , x 2 , …, x K ) and y = (y 1,1 , y 1,2 , …, y K,J ). The sampling interval for c values is 𝐴𝐴 [6,14] . The NN is two layers deep, and has 32 nodes on each. More specific details about the NN architecture are available in Table S2 in Supporting Information S1. The model is trained over 30 epochs using the MSE loss function (Equation 2). The validation data set is made of 15% randomly chosen samples from the learning data set. R 2 score is monitored over this data set during fitting. Best weights regarding the validation R 2 score are saved and loaded after 30 epochs. The final model has R 2 = 0.89 and is noted, 𝐴𝐴 B𝐵(𝑥𝑥𝑥 𝑥𝑥) .

The fitted NN is then used to generate validation time series 𝐴𝐴 [𝐱𝐱] B𝐵 of length 15 MTU, using 𝐴𝐴 B𝐵 instead of B in the L96 equations (Equation 5). Long-term metric 𝐴𝐴 𝐴 𝐴𝐴 is computed over 𝐴𝐴 [𝐱𝐱] B𝐵 , using standard deviation values only:

m𝑚(𝑐𝑐) = ( σ𝜎 (𝑐𝑐) -𝜎𝜎ref ) 2 , ( 9
)
with σ ref (resp. 𝐴𝐴 𝐴 𝐴𝐴 (𝑐𝑐) ) the standard deviation computed over the reference time series (resp. the validation time series). A time series of length 15 MTU, generated by integrating L96 equations (Equations 5 and 6) with (h 0 , F 0 , b 0 , c 0 ) = (1,10,10,10), is considered as the "reference" data set. [7., 13.] (Figure 2). The value of c minimizing the online metric is retrieved on 𝐴𝐴 𝐴 𝐴𝐴 , by using BFGS minimization from SciPy python package. The optimal value of c is c* = 9.922. This value is very close to the value used to generate the reference data set, that is, c 0 = 10.

Imperfect Model Calibration: The Lorenz'96 Model

We now investigate the case where one of the L96 model parameters is carrying biases. To reproduce this situation, the L96 NN parameterization has been trained on output from the L96 model using the reference value of the model parameters, that is, (h, F, b) = (h 0 , F 0 , b 0 ). However, the NN model will be implemented in an L96 system where one of the model parameters is set to a value different from its reference value. We will show that optimizing the value of the tunable parameter included in the NN still allows us to obtain the wished long-term statistics.

Namely, we set the value of F to a range of biased values, F b ≠ F 0 . For each F b , we generate validation time serie using the NN parameterization, 𝐴𝐴 B𝐵(𝑥𝑥; 𝑐𝑐) , and compute the corresponding long-term statistics. To underline the dependence of the long-term statistics on F b , the online metric will be noted:

m𝑚𝐹𝐹 𝑏𝑏 (𝑐𝑐) = ( σ𝜎𝐹𝐹 𝑏𝑏 (𝑐𝑐) -𝜎𝜎ref ) 2 , ( 10
)
where 𝐴𝐴 𝐴 𝐴𝐴𝐹𝐹 𝑏𝑏 is the standard deviation computed on a validation time series where F = F b .

To predict the estimated metric, 𝐴𝐴 𝐴 𝐴𝐴𝐹𝐹 𝑏𝑏 (𝑐𝑐) , again, a kriging metamodel is trained over a data set , for each value of F b . Whereas the loss values remain high for c = c 0 when F is strongly biased, the use of c*(F b ) significantly reduces the online loss value (Figure 3). As soon as F b ≠ F 0 , c*(F b ) ≠ c 0 . This means that the choice of c = c 0 is not optimal. It also means that at least part of the bias induced by the wrong parameter F b can be actually compensated by tuning another parameter of the model, c. This result suggest that, at least to some extent, even an imperfect model can be tuned to adjust long-term statistics.

Figure 3 clearly shows that, as soon as F b ≠ F 0 , the proposed method is a statistically significant improvement to the baseline (linear regression) physical parameterizations and the NN parameterization with c = c 0 . Confidence intervals associated with our estimates of the metric have been constructed by bootstrap through the sampling of two main sources of uncertainty. First, we build 9 different NN parameterizations, each of which has been fitted to a learning sample generated with different initial conditions. Second, each one of the 9 parameterizations are used to generate 3 validation orbits of length 10, 15 and 25 MTUs. Thus, a sample of n = 27 validation orbits is obtained and used to compute 95% confidence intervals for both the "standard" (c = c 0 ) and the "optimal" (c = c*) cases, and for each value of F b .

Conclusion -Discussion

To improve current ML based physics in climate models, the effort is usually only focused on improving both the learning sample and the offline fit of the ML model. But these ML models can and need to be further improved by using an online metric to their calibration [START_REF] Schneider | Earth system modeling 2.0: A blueprint for models that learn from observations and targeted high-resolution simulations: Earth system modeling 2.0[END_REF]. Thus, we propose a method to apply online calibration to ML parameterizations.

The key novelty of our approach is to include some of the physical parameters θ among the input variables of the NN. The NN model is fitted offline to a learning sample of outputs from an high-resolution climate model. In this way, the NN parameterization is able to emulate the physics not only for one single θ 0 , but for a range of values of θ. When the fitted NN model is plugged-in to replace the physical parameterization, the value of parameters θ is calibrated online, as to reduce errors on long-term statistics of the climate model. As a proof-of-concept experiment, our methodology is demonstrated using L63 and L96 models. We show that our method can be used to optimize the value of some parameters to compensate long-term errors due to biases carried by another parameter which cannot be calibrated.

In addition to the reduction of long-term model errors, including some physical parameters among the ML model inputs can also increase the confidence we have in ML parameterizations, the interpretability of which is often questioned. Supplementary parameters can also be used to estimate uncertainties related to some processes. However, although satisfying results were obtained using toy models, further research is needed to test our methodology into real climate models. The generation of a learning sample for subgrid scale parameterizations is already a challenging issue; this task is even more difficult in our new method as a sampling of θ is needed in addition to the sampling of x. The numerical cost of the generation of the learning sample is thus increased. If our approach was to be applied to a climate model, the methodology described for the L96 model could be applied. This would require a set of short high-resolution integrations for several values of θ, for example, taken from a LHS on θ only. In this way the number of model integrations would be kept small, while preserving a large learning sample on (x, θ) (as each integration provides a large sample of x values). In the L96 example, we did not discuss which minimum value of the number of integrations N i is sufficient to efficiently emulate the system. Finding such a minimum value will require careful examination in the case of a climate model. Once a satisfying learning sample is built, the next hurdle to overcome is the fit of the NN model. For more complex models, a simple feed-forward NN may be insufficient and the choice of MSE as the offline metric may not be relevant. Finally, finding the optimal online metric can also be a strenuous issue. (𝑐𝑐0) , are plotted in the top panel (dotted solid black line). The 95% confidence intervals are represented in gray shading. This method can be compared with output using a linear regression approach, fitted to approximate B (red, top panel). (blue) For each F b , we also compute the optimal value c* of c to minimize the online metric

Figure 1 .

 1 Figure 1. (left) The online metric m is computed over orbits of length 1000 Model Time Units (MTUs) generated with real L63 equations (Equation 4) for σ = σ 0 and for different regularly-spaced values of θ = (ρ, β). (middle) Values of 𝐴𝐴 𝐴 𝐴𝐴 in the kriging learning sample, that is, 𝐴𝐴 { 𝜽𝜽𝑖𝑖, m𝑚 (𝜽𝜽𝑖𝑖) } 1≤𝑖𝑖≤𝑁𝑁

Figure 2 .

 2 Figure 2. A Latin Hypercube Sampling (LHS) sample of (x, b) values, 𝐴𝐴 {𝐱𝐱𝑖𝑖, 𝑐𝑐𝑖𝑖} 1≤𝑖𝑖≤𝑁𝑁 𝑖𝑖 , is generated using LHS with c sampled in 𝐴𝐴 [7, 13] . The corresponding values of the long-term metrics, 𝐴𝐴 𝐴 𝐴𝐴 (𝑐𝑐𝑖𝑖) , 1 ≤ 𝑐𝑐𝑖𝑖 ≤ 𝑁𝑁𝑖𝑖 , are computed on 𝐴𝐴 [𝐱𝐱] B𝐵 (length: 3,5 Model Time Units (MTU)) from the LHS initial conditions. The resulting sample 𝐴𝐴 { 𝑐𝑐𝑖𝑖, m𝑚 (𝑐𝑐𝑖𝑖) } 1≤𝑖𝑖≤𝑁𝑁 𝑘𝑘 (light blue scatter points) is used to train the kriging model (solid blue line). For comparison, red dots represent the metric m computed over time series (length: 15 MTU) obtained by applying L96 equations (Equation 5), for discrete values of c. The fitted kriging metamodel, 𝐴𝐴 𝐴 𝐴𝐴 (solid blue line), is a smoothed version of the metric computed on time series generated with the neural network parameterization. The optimal value of c, c*, is computed by minimizing 𝐴𝐴 𝐴 𝐴𝐴 . c* = 9.922 (dashed blue line) can be compared with c 0 = 10 (dashed black line), which was used to generate the reference orbit.

Figure 3 .

 3 Figure 3. (black) Orbits computed with neural network parameterization 𝐴𝐴 B𝐵(𝑥𝑥; 𝑐𝑐) with c = c 0 for each F b in the interval 𝐴𝐴[9, 11] (black bars, bottom panel). The logarithm of kriged values of the metric corresponding to these time series, 𝐴𝐴 𝐴 𝐴𝐴𝐹𝐹 𝑏𝑏(𝑐𝑐0) , are plotted in the top panel (dotted solid black line). The 95% confidence intervals are represented in gray shading. This method can be compared with output using a linear regression approach, fitted to approximate B (red, top panel). (blue) For each F b , we also compute the optimal value c* of c to minimize the online metric 𝐴𝐴 𝐴 𝐴𝐴𝐹𝐹 𝑏𝑏 (blue bars, bottom panel). The 95% confidence intervals are represented in blue shading. The minimal value of each F b , 𝐴𝐴 𝐴 𝐴𝐴𝐹𝐹 𝑏𝑏 (𝑐𝑐 * ) , associated with orbits obtained with machine learning parameterization 𝐴𝐴 B𝐵 (𝑥𝑥; 𝑐𝑐 * ) remains close to zero (dotted solid blue line,

  Figure 3. (black) Orbits computed with neural network parameterization 𝐴𝐴 B𝐵(𝑥𝑥; 𝑐𝑐) with c = c 0 for each F b in the interval 𝐴𝐴[9, 11] (black bars, bottom panel). The logarithm of kriged values of the metric corresponding to these time series, 𝐴𝐴 𝐴 𝐴𝐴𝐹𝐹 𝑏𝑏(𝑐𝑐0) , are plotted in the top panel (dotted solid black line). The 95% confidence intervals are represented in gray shading. This method can be compared with output using a linear regression approach, fitted to approximate B (red, top panel). (blue) For each F b , we also compute the optimal value c* of c to minimize the online metric 𝐴𝐴 𝐴 𝐴𝐴𝐹𝐹 𝑏𝑏 (blue bars, bottom panel). The 95% confidence intervals are represented in blue shading. The minimal value of each F b , 𝐴𝐴 𝐴 𝐴𝐴𝐹𝐹 𝑏𝑏 (𝑐𝑐 * ) , associated with orbits obtained with machine learning parameterization 𝐴𝐴 B𝐵 (𝑥𝑥; 𝑐𝑐 * ) remains close to zero (dotted solid blue line, top panel).

, noted θ*.

  in the computation of m are typically the average and standard deviation values estimated over time series. The minimal value of m is reached at the optimal value of θ

	where 𝐴𝐴	 are a set of M long-term statistics computed over 𝐴𝐴	[𝐱𝐱] f𝑓 , and 𝐴𝐴	 ref a set of reference statistics. The long-
	term statistics involved

  We now optimize the value of c with respect to the online metric 𝐴𝐴 𝐴 𝐴𝐴𝐹𝐹 𝑏𝑏 . For each F b , we find the optimal value of c, noted c*(F b ), by minimizing the corresponding kriging metamodel, 𝐴𝐴 𝐴 𝐴𝐴𝐹𝐹 𝑏𝑏 (Figure3, bottom panel). Bias compensation can be evaluated by comparing the value of

			𝐴𝐴	{ 𝑐𝑐𝑖𝑖, m𝑚𝐹𝐹 𝑏𝑏 (𝑐𝑐𝑖𝑖)	} 1≤𝑖𝑖≤𝑁𝑁 𝑘𝑘
		where c i are N k = 200 values sampled in the interval 𝐴𝐴	[7, 13] . Hence, a kriging metamodel is obtained for each F b
	𝐴𝐴	and the corresponding functions are noted 𝐴𝐴 𝐴𝐴 𝐴 𝐴𝐴𝐹𝐹 𝑏𝑏 (𝑐𝑐) . 𝐴 𝐴𝐴𝐹𝐹 𝑏𝑏 computed on validation time series generated with 𝐴𝐴 B𝐵 (𝑥𝑥; 𝑐𝑐	B𝐵 (𝑥𝑥; 𝑐𝑐0)

* )

Acknowledgments

We are grateful for the insightful comments and suggestions made by Olivier Geoffroy.

Data Availability Statement

Code is made available at: https://zenodo.org/record/6141165.


Jean-Marie Lalande
email: jean-marie.lalande@meteo.fr

Guillaume Bourmaud

Pierre Minvielle

Jean-François Giovannelli

A kriging-based analysis of cloud Liquid Water Content using CloudSat data

Spatiotemporal statistical learning has received increased attention in the past decade, due to spatially and temporally indexed data proliferation, especially collected from satellite remote sensing. In the mean time, observational studies of clouds are recognized as an important step to improve cloud representation in weather and climate models. Since 2006, the satellite CloudSat of NASA carries a 94 GHz cloud profiling radar and is able to retrieve, from radar reflectivity, microphysical parameter distribution such as water or ice content. The collected data is piled up with the successive satellite orbits of nearly two hours, leading to a large compressed database of 2 Tb (http://cloudsat.atmos.colostate.edu/). These observations give the opportunity to extend the cloud microphysical properties beyond the actual measurement locations using an interpolation and prediction algorithm. In order to do so, we introduce a statistical estimator based on the spatiotemporal covariance and mean of the observations known as kriging. An adequate parametric model for the covariance and the mean is chosen from an exploratory data analysis. Beforehand, it is necessary to estimate the parameters of this spatiotemporal model; This is performed in a Bayesian setting. The approach is then applied to a subset of the CloudSat dataset.

Introduction

Clouds have a strong influence on weather and climate. They are a key element of Earth's hydrological cycle, bringing water from the air to the ground and from one region of the globe to another. They also dominate the energy budget of the Earth through their action on the exchange of solar and thermal radiation within the atmosphere and between the atmosphere, the hydrosphere, the land surface, the biosphere, and space. However, they still remain a major source of uncertainties in predicting the weather and climate change.

While the measurement of cloud occurrences and properties at useful spatial and temporal scales is notoriously difficult [START_REF] Marshak | 3D Radiative Transfer in Cloudy Atmospheres, Physics of Earth and Space Environments[END_REF][START_REF] Stephens | The Remote Sensing of Clouds and Precipitation from Space: A Review[END_REF], the proliferation of satellite platforms in the last decades (Stephens et al., 2002;[START_REF] Eriksson | Comparison between early Odin-SMR, Aura MLS and CloudSat retrievals of cloud ice mass in the upper tropical troposphere[END_REF][START_REF] Wu | Comparisons of global cloud ice from MLS, CloudSat, and correlative data sets[END_REF] is fostering a number of new approaches. One such satellite, CloudSat, whose payload is a Cloud Profiling Radar (CPR), has been dedicated to measure the cloud vertical structure and microphysical properties. It is part of the A-train constellation, who was originally set to comprise 7 satellites specifically designed to measure cloud and precipitation properties using different instruments. Since 2006 CloudSat has collected a large database of cloud properties, globally and over an extended period of time, despite some malfunctions.

1

In this study, we propose to analyze statistically a part of this database in order to perform interpolation and prediction of cloud properties. This kind of spatio-temporal estimation can be of major importance for the assessment of satellite cloud attenuation [START_REF] Lyras | Cloud attenuation statistics prediction from Ka-band to optical frequencies: Integrated liquid water content field synthesizer[END_REF]r, including Global Positioning System (GPS) radio occultation [START_REF] Yang | Assessments of cloud liquid water contributions to GPS radio occultation refractivity using measurements from COSMIC and CloudSat[END_REF]. It is required for improving the representation of cloud systems in numerical weather prediction [START_REF] Bodas-Salcedo | Evaluating cloud systems in the Met Office global forecast model using simulated CloudSat radar reflectivities[END_REF][START_REF] Chen | Partitioning CloudSat ice water content for comparison with upper tropospheric ice in global atmospheric models[END_REF], e.g. for data assimilation [START_REF] Qu | Evaluation of a high-resolution numerical weather prediction model's simulated clouds using observations from CloudSat, GOES-13 and in situ aircraft[END_REF]. It can be involved in the assessment of an aircraft icing detection system (Vivekanandan et al., 2001), the design of a satellite communication system [START_REF] Khan | Empirical Relations for Optical Attenuation Prediction from Liquid Water Content of Fog[END_REF], the systematic comparison with other cloud products from different instruments/satellites, etc.

As observations from different instruments are likely not collected at the same spatiotemporal positions, the procedure of interpolation-prediction usually involves regridding selected slots of data. A number of statistical methods dedicated to the analysis of spatial and spatio-temporal data have been developed over the years taking into account the spatial and/or temporal correlation of the observations [START_REF] Ripley | Spatial statistics[END_REF][START_REF] Cressie | Statistics for spatial data[END_REF]. In this study, we propose to use an approach based on a kriging estimator for the interpolation/prediction problem. The kriging estimator was initially introduced by [START_REF] Krige | A statistical approach to some basic mine valuation problems on the Witwatersrand[END_REF], from which it takes its name, to estimate the gold distribution at the Witwatersrand reef complex in South Africa based on samples from boreholes. It was then formalized mathematically by [START_REF] Matheron | Principles of geostatistics[END_REF] in the context of mining geology. Afterwards, the kriging estimator spread to many other areas of sciences (Wackernagel, 2013) (hydrogeology, geotechnics, agronomy, air quality, fishery, epidemiology, water and soil pollution, noise, etc.). Best-known kriging techniques are the Simple Kriging, which assumes stationarity of the 1 st -order with a known mean, the Ordinary Kriging, where the mean is unknown and to be estimated, and the Universal Kriging, which assumes non-stationarity of the 1 st -order. Since its first developmnent, kriging techniques have largely evolved, and a number of new kriging techniques have been developed [START_REF] Chiles | Geostatistics: modeling spatial uncertainty[END_REF][START_REF] Cressie | Statistics for spatial data[END_REF][START_REF] Cressie | Statistics for Spatio-Temporal Data[END_REF]. In the field of meteorology, the kriging estimator have primarily been used to estimate precipitation accumulation from rain gauges [START_REF] Nour | Geostatistical mapping of precipitation: implications for rain gauge network design[END_REF][START_REF] Belo-Pereira | Evaluation of global precipitation data sets over the Iberian Peninsula[END_REF] and in combination with satellite-derived precipitation (Jewell and Gaussiat, 2015;[START_REF] Verdin | Kriging and Local Polynomial Methods for Blending Satellite-Derived and Gauge Precipitation Estimates to Support Hydrologic Early Warning Systems[END_REF][START_REF] Varouchakis | Combining Geostatistics and Remote Sensing Data to Improve Spatiotemporal Analysis of Precipitation[END_REF] as well as aerosols concentration in the air from in-situ observations [START_REF] Park | Time-Series Mapping of PM10 Concentration Using Multi-Gaussian Space-Time Kriging: A Case Study in the Seoul Metropolitan Area, Korea[END_REF]. It has also been used for the estimation of temperature from in-situ measurements [START_REF] Heuvelink | Sampling Design Optimization for Space-Time Kriging, Spatio-Temporal Design: Advances in Efficient Data Acquisition[END_REF][START_REF] Didari | Enhancing estimation accuracy of daily maximum, minimum, and mean air temperature using spatio-temporal ground-based and remote-sensing data in southern Iran[END_REF], from satellite observations (Florio et al., 2004) or a combination of them [START_REF] Didari | Enhancing estimation accuracy of daily maximum, minimum, and mean air temperature using spatio-temporal ground-based and remote-sensing data in southern Iran[END_REF] as well as for the estimation of surface properties from remote measurements [START_REF] Der Meer | Remote-sensing image analysis and geostatistics[END_REF][START_REF] Zakeri | A review of geostatistical simulation models applied to satellite remote sensing: Methods and applications[END_REF].

Our approach is based on a second-order analysis of cloud Liquid Water Content (LWC) obtained from Level 2B product of the CloudSat ground segment data. We perform a detailed analysis of those properties and propose an adequate parametric model for their mean and covariance. The model parameters are obtained by maximizing the a posteriori probability density which comprises a likelihood term expressing how the model fits the observation conditionnally to the presence of clouds.

Finally we apply this model in the context of interpolation and prediction of CloudSat observations. To our knowledge, it is the first time the kriging estimator have been employed to infer cloud microphysical properties using CloudSat observations :

-We provide a comprehensive mathematical description of the kriging estimator. It allows us to mention that the parameters of the mean and covariance functions are commonly treated differently in the kriging literature. In our approach, parameters can enter both linearly and non-linearly in the mean, which makes it more flexible in accommodating the trend of the observations. Thus, we propose a global treatment of all the parameters of the model as the standard Universal Kriging is not compatible with non-linearity.

-We propose a model for the mean and covariance functions and thoroughly justify each choice we make (e.g., stationarity assumption, homogeneity assumption, etc.) -We perform a detailed analysis of both the estimated model parameters and the resulting interpolated/predicted cloud LWC. In this study, we use the level 2B Cloud Water Content product (2B-CWC-RO) that contains retrieved estimates of cloud Liquid and Ice Water Content, effective radius, and related quantities for each radar backscattered reflectivity profile [START_REF] Austin | Retrieval of ice cloud microphysical parameters using the CloudSat millimeter-wave radar and temperature[END_REF]. We focus on a region centered over Europe and under-sampled at a rate of 1/50 (cf. Fig. 1). In this region, we select data from June 16-th 2006 to June 14-th 2015 constituting a set of 239087 profiles. The CPR observations are classified according to their cloud types in the level 2B-CLDCLASS product (cf. Fig. 2). After analysis of the cloud types distribution inside the considered region, we decided to focus on the LWC labeled as Altostratus as they constitute a good balance between the amount of available data and the computational feasibility. Moreover, altostratus are less fractionated than other cloud types, thereby their physical properties are more continuous. 2.2 Mathematical description of the data 85

Background

We are considering the scalar function Φ of three spatial variables (x, y, z) and a temporal variable t mapping elements from D into R:

Φ : D = [0, 2π] × [-π/2, π/2] × R + × R → R s = (x, y, z, t) → Φ(x, y, z, t), (1)
where s ∈ D is spatiotemporal localization, x, y represent respectively the longitude and latitude at the Earth's surface and z 90 is the altitude. The function Φ can be modeled as a deterministic or stochastic quantity. The stochastic nature of the function Φ can be introduced to model some inner variability of the physical phenomenon under study or an incomplete knowledge of the phenomenon itself. The complete description of a stochastic process requires the construction of the joint probability density function for the continuous variables in space and time [START_REF] Chonavel | Statistical Signal Processing: Modelling and Estimation, Advanced Textbooks in Control and Signal Processing[END_REF][START_REF] Gaetan | Modélisation et statistique spatiales, Mathématiques et Applications[END_REF][START_REF] Brockwell | Time Series: Theory and Methods[END_REF], which can be cumbersome in a high-dimensional space. To alleviate this burden, it is possible to resort to the description of the mean and covariance of the random function. This proves very efficient when the distribution associated to the random function is not too complex (for instance when the random function is Gaussian, the mean and covariance are sufficient to completely describe the distribution).

In this study, we assume that the mean and covariance exist and represent sufficiently well the distribution of the random function Φ. The mean is a function of the spatiotemporal localization s:

m Φ (s) = E Φ(s) ,
where E • denotes the expectation of the random variable. The covariance is a function of 4 spatiotemporal variables and 4 spatiotemporal shift variables, noted δ = (δ x , δ y , δ z , δ t ). Thus we have,

R Φ (s, δ) = E (Φ(s) -m Φ (s))(Φ(s + δ) -m Φ (s + δ)) .
The covariance function describes the spatiotemporal dependency of the stochastic function Φ. Based on these definitions, it is useful to introduce certain stationary model to characterize the degree of homogeneity of the random function Φ.

Stationarity

In order to infer the distribution of the random function Φ, the structure of the spatiotemporal process is of great importance.

In this sense, stationary model enable to structure the spatiotemporal variability of the random function Φ. The process under study is said to be strictly stationary if the distributions [Φ(s 1 ),

• • • , Φ(s k )] and [Φ(s 1 + δ), • • • , Φ(s k + δ)]
are identical for all δ, s 1 , . . . , s k ∈ D, and for all k ∈ N. This is a very restrictive condition which is not usually satisfied in real-life applications.

The weak stationarity stipulates that:

the random function Φ is 1-st order stationary when its expected value does not depend on the localization s:

m Φ (s) = µ
the random function Φ is 2-nd order stationary when its covariance function depends only on the lag vector δ between 2 localizations s and s = s + δ: R Φ (s, δ) = R Φ (δ), for all s and δ .

Observations and uncertainties

We denote the successive spatiotemporal positions s n = (x n , y n , z n , t n ), n = 1, • • • , N and the targeted quantity Φ at these positions:

Φ n = Φ(s n ), n = 1, • • • , N .
We model the actual observations of the targeted quantity Φ as corrupted by an additive noise related to the measurement uncertainty. The observation at position s n is then denoted Ψ n and is written:

Ψ n = Φ n + B n ,
where Ψ n , Φ n and B n are random variables. We denote the values associated with a realization of these random variables by ψ n , ϕ n and b n . Thus, a realization of the random variable Ψ at the position s n is written:

ψ n = ϕ n + b n , (2)
and additionally, we write:

ψ = ϕ + b, where ψ = [ψ 1 , • • • , ψ N ], ϕ = [ϕ 1 , • • • , ϕ N ] and b = [b 1 , • • • , b N ]
are the collections of the related N quantities.

Objectives

The objective of this work is to determine the value ϕ 0 of the quantity of interest ϕ at a given location s 0 = (x 0 , y 0 , z 0 , t 0 ). This is an estimation problem that is tackled by designing an estimate φ0 of the quantity of interest ϕ 0 from the data ψ described in the previous section. On the one hand, spatially, this is an interpolation problem since the point x 0 , y 0 , z 0 is usually not on the grid of observations. On the other hand, temporally, this is a prediction problem since t 0 is naturally positioned in the future. The considered approach is to statistically learn from the large amount of available data both to interpolate/predict the targeted quantity and to quantize the uncertainty associated with the estimated value. In the next section, we introduce the kriging estimator, that is specifically designed to perform such a task.

The kriging estimator

The goal is to construct a statistical estimator Φ0 Φ0 :

R N → R ψ → φ0 = Φ0 (ψ)
in order to estimate ϕ 0 from the observations ψ. This section describes the strategy for the construction of the kriging estimator which has been widely used in geostatistics. This framework has been developed in a handful of handbooks [START_REF] Cressie | Statistics for spatial data[END_REF][START_REF] Chiles | Geostatistics: modeling spatial uncertainty[END_REF][START_REF] Diggle | An Introduction to Model-Based Geostatistics[END_REF] and later on by [START_REF] Montero | Spatial and Spatio-Temporal Geostatistical Modeling and Kriging[END_REF], in the case of a known mean and covariance. Our development differs slightly from the standards of geostatistical literature in two ways.

1. The parameters of the mean function m Φ (s) usually act linearly and consequently their estimation can be performed inside the kriging estimator while the estimation of the covariance parameters is classically performed in a previous offline procedure. Instead, in order to design a model with a higher capacity, we consider both a non-linear mean and a non-linear covariance functions and estimate their parameters in an offline procedure.

2. Our approach is fully parametric and consequently relies on a parametric covariance function instead of a variogram [START_REF] Banerjee | Hierarchical Modeling and Analysis for Spatial Data, Second Edition[END_REF].

Thus, in the following development, we consider that the mean and covariance functions are known and defer the (offline) estimation of their parameters to the next section.

We consider a linear estimator of the form:

Φ0 (Ψ) = a t Ψ + a 0 , (3)
where

Ψ = [Ψ 1 , • • • , Ψ N ] are the observations, a = [a 1 , • • • , a N ]
are scalar coefficients and N is the number of available observations. This can be alternatively expressed in terms of realizations of the random variables as:

φ0 = Φ0 (ψ) = a t ψ + a 0 . ( 4
)
Our goal is to determine the coefficients a and a 0 . The strategy is to minimize an error between the estimated quantity Φ0 and the true value Φ 0 . We choose the Mean Squared Error (MSE),

E(a, a 0 ) = E ( Φ0 -Φ 0 ) 2 . ( 5
)
The estimator Φ, defined by ( 3)-( 4) with:

(a, a 0 ) opt = arg min

a opt 0 = m Φ0 -a t m Φ • (6)
Then, plugin (6) into ( 5) and expanding leads to:

E(a) = E(a, a opt 0 ) = a t R Ψ a -2a t r ΨΦ0 + var Φ 0 (7)
where we introduced the following notations:

R Ψ = E (Ψ -m Φ )(Ψ -m Φ ) t r ΨΦ0 = E (Φ 0 -m Φ0 )(Ψ -m Φ ) var Φ 0 = E (Φ 0 -m Φ0 ) 2 ,
R Ψ is the covariance matrix of the observations, r ΨΦ0 the covariance between the observations and the quantity of interest, and var Φ 0 is the variance of the quantity of interest. The vector a opt that minimizes (7) satisfies

∂E ∂a a opt = 0,
and that leads to:

a opt = R -1 Ψ r ΨΦ0 .
Moreover, modeling noise B independent of the quantity of interest Φ, we have

R Ψ = R Φ + R B and r ΨΦ0 = r ΦΦ0 , (8)
so that,

a opt = (R Φ + R B ) -1 r ΦΦ0 . (9)
Finally, we deduce the estimated value φ0 of the quantity of interest ϕ at position s 0 by inserting equations ( 9) and ( 6) in ( 4):

φ0 = r t ΦΦ0 (R Φ + R B ) -1 (ψ -m Φ ) + m Φ (s 0 ) . (10)
Additionally, the MSE of the estimated value is given by:

E(a opt , a opt 0 ) = var Φ 0 -r t ΦΦ0 (R Φ + R B ) -1 r ΦΦ0 , (11)
which, in this context, equals to the variance of the estimation error var Φ0 -Φ 0 because the bias is equal to zero. This estimator is known as the kriging estimator in the geostatistical literature [START_REF] Montero | Spatial and Spatio-Temporal Geostatistical Modeling and Kriging[END_REF]. In a slightly different form, this estimator is also the one of Kalman and Wiener in the field of filtering.

4 Construction of a model for the mean and the covariance 4.1 Time-altitude and spectral analysis

In order to construct an adequate model for the mean and the covariance, we start with an exploratory analysis to extract some raw characteristics from the time-altitude observations represented in Fig. 3. We compute the biased empirical covariance for the 2-dimensional case (δt, δz) and the 1dimensional case at a fixed altitude of z = 6.711 km (cf. Fig. 4). This shows 1 year periodicity of the LWC observations in the time component, while the altitude component decreases to 0 for δz = 1 km. The blue spots located at δ z = +/ -6 km are due to zero-padding above and below the cloudy profile that introduces an artificial correlation. Another way of looking at the data ψ is to compute the empirical power spectral density Ŝψ (f ), i.e., the squared modulus of the Fourier transform of the observations ψ (a.k.a. periodogram) (cf. Fig. 5). There is clearly a peak at 1 year -1 .

For some altitudes, there are lower peaks at 2 year -1 or even at 3 year -1 , indicating that the signal cannot be modeled by a single sinusoidal component. At lower altitudes, the fundamental component becomes negligible as well as the harmonics.

Nevertheless, this suggests that CloudSat observations can be interpreted as the superposition of multiple periodic components.

Proposition of an adequate model

The observations Ψ are modeled as the addition of a random process Φ and a random noise process B, such as:

200 Ψ(t, z) = Φ(t, z) + B(t, z).
We dodge the complexity of handling 4 dimensions and choose to do an in-depth analysis of a time-altitude dependent model while omitting the latitude-longitude dimensions. This simplification can be justified when the spatial extension of the dataset is not too important. This is equivalent to assimilate the dataset to a single spatial location by averaging over the geographic components. All random processes are modeled as Gaussian [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF]: In order to take into account the periodic temporal trend observed in the observations, we consider a non 1-st order stationary model. Moreover, the mean is chosen not to be a function of the altitude in order to limit the number of parameters. Considering the time variations, we chose a periodic mean with 1 fundamental and 2 harmonics such that:

Φ ∼ N (m Φ , R Φ ), B ∼ N (0, R B ), Ψ ∼ N (m Φ , R Ψ ), (12)
m Φ (t) = β 0 + β 1 cos(α 1 + 2πt T ) + β 2 cos(α 2 + 2πt T /2 ) + β 3 cos(α 3 + 2πt T /3 )
where β 0 , β 1 , β 2 and β 3 are respectively the amplitude of the continuous component, the fundamental and the 2 harmonics, α i (with i = 1, 2, 3) are the phase parameters of the fundamental and the harmonics and T is the period.

In order to further simplify we resort to the widely used separable model for the covariance [START_REF] Genton | Separable approximations of space-time covariance matrices[END_REF]) so that the spatiotemporal covariance factors into a purely spatial and a purely temporal component. This reduces the number of parameters allowing for computationally efficient estimation and inference, it is then more tractable to run different test on the dataset.

Thus, the separable covariance model writes

R Φ (δ t , δ z ) = r Φ R Φt (δt) R Φz (δz),
where r Φ is the variance of the quantity of interest. We chose an exponential covariance model for the time covariance

R Φt (δt) = exp [-|δt|/l t ],
where l t is the correlation time. The choice of this covariance model is motivated by the shape of the empirical covariance described in the previous section (it decreases more rapidly than a Gaussian). The altitude covariance is set to be:

R Φz (δz) = exp -(δz/l z ) 2 ,
where l z is the correlation length in altitude. Finally, we use a white noise:

R B (δt, δz) = r B Dirac(δt, δz),
where r B is the variance of the noise. Finally, we have a model consisting of 12 parameters, including 4 parameters for covariance, and 8 parameters for the mean

θ = [r Φ , l t , l z , β 0 , β 1 , β 2 , β 3 , α 1 , α 2 , α 3 , T, r B ] .
Using this notation we can denote the parametric model for the mean m Ψ (δt; θ) and for the covariance R Ψ (δt, δz; θ). In the following section, we describe the strategy to estimate this set of parameters.

Model parameter estimation, optimization and final kriging equations

In this section, we present the strategy for the estimation of the model parameters θ. We choose the Maximum A Posteriori (hereafter abbreviated MAP) estimate for θ which is the mode of the posterior probability density function (pdf) in Bayesian statistics. The analysis of the estimated parameters and corresponding model is deferred to Sect. 6 for the 1-dimensional case and Sect. 7 for the 2-dimensional case.

The MAP estimator

In the previous section, we proposed a parametric model for the mean m Ψ (δt; θ) and for the covariance R Ψ (δt, δz; θ) of the random variables Ψ, which decomposes into a model for the quantity of interest Φ and the noise B. The problem is to find an estimate θ of the true parameters θ from the observations ψ. We introduce the posterior pdf:

π(θ|ψ) = f (ψ|θ)ρ(θ) f (ψ) ,
where f (ψ|θ) is the so-called likelihood function, which expresses the probability of the observations given the parameter θ, ρ(θ) is the prior distribution for the parameters θ and f (ψ) is the marginal pdf for the observations. Since we model the observations by a Gaussian pdf conditionally to the parameters θ (see ( 12)), the likelihood function f (ψ|θ) writes:

f (ψ|θ) = (2π) -N/2 (det R Ψ ) -1/2 exp - 1 2 (ψ -m Ψ ) t R -1 Ψ (ψ -m Ψ ) , (13)
where we have omitted the dependence in θ, δt and δz for R Ψ and m Ψ to weigh down the notations. Since the posterior pdf π(θ|ψ) is by definition positive, we consider the co-log posterior CLP(θ) = -log π(θ|ψ), which expands as follows:

CLP(θ) # -log f (ψ|θ) -log ρ(θ) # 1 2 ψ -m Ψ t R -1 Ψ ψ -m Ψ + 1 2 log det R Ψ -log ρ(θ) (14)
where # stands for "equality up to an additive constant" (not function of θ). The minimum of ( 14) is the maximum of the posterior pdf. Thus, minimizing the CLP gives the MAP estimate. To complete, we consider a uniform prior ρ(θ) on a given domain Θ:

ρ(θ) = U Θ (θ)
where Θ is a hyper-rectangle in R P describing a range of possible values for each component of θ. For our model P = 12.

Finally, we obtain:

CLP(θ) # (ψ -m Ψ ) t R -1 Ψ (ψ -m Ψ ) + log det R Ψ -log U Θ (θ) (15)
The function CLP is infinite when θ / ∈ Θ because of the term -log U Θ (θ). The minimizer θ, given by:

θ = arg min θ CLP(θ), (16)
is the MAP estimate of the true parameter θ .

Optimization procedure

An optimization procedure is required in order to determine the minimizer θ of the CLP. In this study, we use an Iterative

Conditional Mode strategy with a Golden-Section search (ICM-GS) for each component [START_REF] Kiefer | Sequential minimax search for a maximum[END_REF][START_REF] Press | Numerical Recipes in C (2nd Ed.): The Art of Scientific Computing[END_REF]. We then have an algorithm composed of an outer and an inner loop.

The inner loop optimizes w.r.t. one component at a time, say θ p , based on a golden-section search. It stops when the variation of θ p is smaller than a given threshold, say ε p . The outer loop repeats the scan of the P components and stops when the variation of θ becomes smaller than a second given threshold denotes by η. The values ε p and η have been fixed so as to reach stable results. The specific threshold values used in optimization procedure are reported in Table 1.

Table 1. Threshold values p used for each parameters and global threshold value η.

r Φ l t lz r B β 0 β 1 β 2 β 3 α 1 α 2 α 3 η 1e -5 1e -7 1e -1 1e -6 1e -3 1e -2 1e -2 1e -4 1e -2 1e -5 1e -2 1e -3
The interpolation/prediction solution given the observations ψ is performed using the kriging estimator developed in Sect.

3. The kriging equations were initially developed assuming a known mean m Φ and covariance R Ψ . In reality, we only have at hand the parametric mean mΦ and covariance RΨ that depend on the estimated parameter θ. Thus, the estimation of the quantity of interest Φ in t 0 is written:

φ0 = rt ΨΦ0 R-1 Ψ (ψ -mΦ ) + mΦ (t 0 ), ( 17
)
where mΦ is a vector of estimated mean at the positions of the observations. Finally, we can compute the minimum mean square error:

E = rΦ -rt ΨΦ0 R-1 Ψ rΨΦ0 . ( 18
)
This last expression is equal to the variance of the estimation error, which we note σ 2 = var Φ0 -Φ 0 . It gives an estimation of the forecast error associated with the kriging technique. However, it should be kept in mind that this variance does not account for the uncertainty in the model parameter estimate θ. Therefore, this expression tends to underestimate the true estimation variance. This feature has been well documented in [START_REF] Cressie | Statistics for spatial data[END_REF][START_REF] Montero | Spatial and Spatio-Temporal Geostatistical Modeling and Kriging[END_REF]. We leave the assessment of this impact for further development of this model. Note that the kriging stage is essentially a computational step that does not present any major difficulty except for the computational burden when dealing with large dimensions.

The 1-dimensional time case

We start our analysis with the 1-dimensional case at a single altitude. We give a detailed analysis of the estimated parameters according to the MAP criterion and the ICM-GS algorithm, both described in the previous section. We concentrate on the estimation over the west of Europe and analyse the impact of reducing the considered geographic area. We then proceed to the kriging of observations in the interpolation and prediction case.

Estimation over a European area at a single altitude

In order to analyze the results of the previously described parameter estimation procedure, we first consider observations associated to altostratus clouds over an European area extending from 10 o W to 10 o E in longitude and 30 o N to 60 o N in latitude (see Fig. 1) and use a subsampled database (at 1/50 rate) at a single altitude of z 6 km. This corresponds to a set of N obs = 3653 observations. The period parameter T has been readily set to 1 year so as to represent the seasonality observed in the data.

Fig. 6 shows the evolution of CLP during the optimization process which, as expected, decreases with the iterations and stabilizes after ∼ 5 iterations. The current value of the parameters changes at each iteration and stabilizes as shown in Fig. 7.

Even though the decrease of CLP along the iterative optimization process indicates convergence towards a minimum, there is no guarantee that it is a global minimum. This is highlighted in Fig. 8(a) which represents the isocontour plot of CLP as a ). The estimated value lt = 3.75×10 -6 corresponds to ∼ 2 minutes. This correlation time is very small compared to the 9-year period used to train the parameters. It indicates that there is essentially no correlation for observations later than 6 minutes apart. Indeed, the exponential covariance model reaches a 5% correlation at about 3l t [START_REF] Banerjee | Hierarchical Modeling and Analysis for Spatial Data, Second Edition[END_REF]. This feature can be either explained by the model, the observations or the phenomenon itself. However, it is then evident that such a model will not be very efficient in the perspective of long-term forecast.

As stated in [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF], each local minimum corresponds to a particular interpretation of the data, and a careful analysis may be required to choose one model instead of the other. In the next section, we will further discuss the presence of two minima of the CLP in relation with two competing models. An exhaustive search for all local minima requires to compute the criterion CLP(θ) for all possible combinations of parameters in the domain Θ, which is impossible. Instead, we performed a set of 10 optimizations with random initialization to track down additional local minima and converged toward the same minimizer. Hence we conclude there is no other local minima than the one previously pointed out (cf. Fig. 8(a)).

Discussion on the reduction of the geographic area

In order to assess the influence of the considered geographic area, the optimization procedure is run by reducing the geographic extent of the dataset. We consider 9 areas that are represented in Fig. 9, and the estimated parameters θ are summarized in Fig. 10.

The estimated values of the covariance parameters (i.e., rΦ , lt and rB ) vary in a complementary fashion, with approximately similar values for zones 1 to 3, a leap of several orders of magnitude for zones 4 to 8, followed by a drop for the 9-th zone.

These variations can be traced back to the presence of a local minimum in the domain Θ as it is shown in Fig. 11 where the variation of CLP -µ CLP with respect to l t has been represented for each geographic area. All CLPs show the presence of a local and global minima on the considered interval (except for the 9-th area). The variations of CLP w.r.t. the correlation time l t have similar behavior for areas 1 to 3 with locations of the local and global minimum roughly identical. The same is valid for areas 4 and 5. Concerning areas 6 to 8 (cf. Fig. 11(c)), they have a more complex shape with a well defined global minimum but less pronounced local minima at lesser correlation time.

It is also worth highlighting the sliding of the global minimum to the local minimum when switching from zones 1, 2 and 3 to zones 4 and 5 (cf. Fig. 11 (a) and (b)). The global minimum of zones 1, 2 and 3 is reached for a smaller correlation time l t ∼ 10 -6 year which means that the aggregation of spatially distant observations is interpreted as a less correlated process than for smaller geographic area (i.e., zones 4 and 5). In this case, the retrieved model tends towards a white Gaussian noise as the covariance function is close to a Dirac. However, the presence of a local minimum of CLP with a higher correlation time indicates that the observations could also be explained by a second model, provided some additional constraint on the correlation time (e.g., a different prior distribution for l t ). In particular, we note that the order of magnitude of the l t values associated with this local minimum is similar to the l t values obtained for zones 4 to 7. When we reduce the geographic area, the set of observations tends to be more homogeneous, which plays in favor of models with higher correlation times.

The result obtained for zone 9 has a different interpretation. The estimated correlation time is lt 7.13 × 10 -7 , which is several orders of magnitude smaller than every other geographic areas. This is in conjunction with a particularly low estimated noise variance rB rΦ (cf. Fig. 10 (c)). For this geographic area, the estimated rB sticks to the minimum bound of the prior density so the optimization procedure is stuck for this parameter. Thus, the MAP estimator compensates for the remaining signal variance by increasing rΦ and decreasing the associated lt , so that R Φ tends towards a Dirac function, i.e., a covariance

(1) (2) (3) (4) (5) (6) (7) (8) (9)
Figure 9. The 9 geographic areas used for the parameter estimation analysis.

associated to white noise. In other words, the estimator is no longer able to separate the quantity of interest Φ from the noise with this dataset. To account for computational burden and stability of the estimated model we will pursue our analysis with those obtained for the zones 1 to 3 and a greater noise variance r B . However, the correlation time corresponding to zones 1 to 3 is very small that the corresponding covariance tends toward a diagonal matrix (i.e., a Dirac function).

Results in interpolation

We start with the interpolation of the LWC over the 9-year period. Note that, despite the regular sampling rate of the CPR, actual observations of altostratus in the European area are sparsely distributed. The temporal positions t 0 at which kriging is applied are defined by t 0 = [0, 1/N 0 , • • • , 9-1/N 0 , 9], with N 0 = 4001. The result of this kriging interpolation is represented in Fig. 12. The variability of the estimated φ0 are smaller than the actual variability of the observations, which is consistent with the estimated model parameters. Indeed, the estimated noise variance rB = 2.2 × 10 -3 g 2 .cm -6 is greater than the estimated variance of the quantity of interest rΦ = 9.8 × 10 -4 g 2 .cm -6 . The periodic nature of the estimate φ0 is strong, especially when 345 there are no observations (for example during the 5-th year of operation of the satellite). In this case, the result of the kriging interpolation consists in the estimated mean mΦ . In the vicinity of observations, the kriging estimates take a more complex form because of the additional information brought by the neighbouring data. The behavior observed over this period of 9 years is not surprising in itself because of the relative weakness of the observed correlation time, which is only lt = 2.7 days. We perform another kriging interpolation over a period of ∼ 36 days centered on the beginning of the satellite's 3-rd year of operation (cf. Fig. 13(a)). In this figure, we observe that some observations are very closed temporally and have significantly different values of LWC which can explain the high estimated variance of the noise with respect to the variance of the quantity of interest. In this case, the kriging gives an average value of the observations. These observations generally correspond to successive profiles in the 1/50 database, which are separated by 8 seconds. It should be noted that between two successive observations, the satellite moves ∼ 54 km which can certainly explain a loss of spatial correlation not taken into account in our 355 model. It is therefore likely that the model will explain this loss of spatial correlation by noise. On the other hand, we notice the presence of peaks in the structure of φ that can be explained by the structure of the chosen temporal covariance, which is in fact characterized by an exponential decay corresponding to the decay observed on the curve of φ when moving away from an observation (cf. Fig. 13 (a)). On all kriging results, we represent the 4σ area, where σ 2 = var Φ0 -Φ 0 is the variance of the estimation error of (18).

360

This variance is decomposed into two terms, the variance var Φ0 and a quadratic term in rΨΦ0 (= rΦΦ0 ), which represents the covariance between the object of interest and the observations. At observation location (i.e., Φ 0 is chosen collocated with an observation), rΨΦ0 reaches its maximum value and so does the quadratic term rΨΦ0 R-1 Ψ rΨΦ0 . Consequently, at observation location the minimum mean square error given in (18) is minimal.

The 2-dimensional time-altitude case

In a similar fashion as Sect. 6, we examine the 2-dimensional time-altitude case starting from the parameter estimation followed by the kriging results.

Parameter estimation

We are considering here the estimation of model parameters including the altitude dimension. Following the conclusions of the previous section, we use the observations of the 6-th geographic area. However, for computational reasons, the training dataset is restricted to the first 3 years of observations, which consists of N obs = 4167 observations (cf. Fig. 17).

Fig. 16 represents the evolution of each parameter during the optimization process. We observe that all parameters have stabilized towards an estimated value. Moreover, none of the parameters have converged towards the bound defined by the prior.

405

The estimated values in the 1-dimensional (temporal) and 2-dimensional (time-altitude) cases are consistent from one model to the other (cf. Table 2). Overall, the order of magnitude for most quantities are close. However, we point out some noticeable differences:

Table 2. Comparison between estimated parameters for the 1d and the 2d model. Note that we have fixed T = 1 in the optimization procedure.

θ 1-D Model 2-D Model rΦ ∼ 9.84 × 10 -4 ∼ 5.3 × 10 -3 lt ∼ 7.5 × 10 -3 ∼ 5.6 × 10 -3 lz / ∼ 865 rB ∼ 2.2 × 10 -3 ∼ 1.7 × 10 -3 β0 ∼ 6.7 × 10 -2 ∼ 6.2 × 10 -2 β1 ∼ 4.2 × 10 -2 ∼ 7.9 × 10 -3 β2 ∼ 2.9 × 10 -3 ∼ 4.4 × 10 -3 β3 ∼ 9.7 × 10 -3 ∼ 3.4 × 10 -3 α1 ∼ 5.4 ∼ 5.5 α2 ∼ 6.1 ∼ 1.5 α3 ∼ 6.1 ∼ 3.6
there is a factor 10 between the estimated variance rΦ of the time model and the time-altitude model;

for the 1-dimensional case, rΦ < rB whereas in the 2-dimensional case we have rΦ > rB , which means adding a dimension to the model helps to find some structure in the data that were missing otherwise;

the amplitudes and phases of the first and second harmonics ( β2 , β3 , α2 and α3 ) differ substantially.

This last point explains the distinct mean behavior for the two models. The mean of altitude-time model has a lower amplitude which can be explained by the fact that the model must accommodate for a greater variability of observations due to an underlying altitude dependence of the dataset. It is an indication that future development should include this altitude dependence.

Kriging results

In this section, we present the kriging results obtained for the time-altitude model. The parameters corresponding to this model have been estimated in Sect. 7.1. The kriging estimator is applied on a regular grid of 20 × 500 positions in time and altitude, this corresponds to a total of 10000 positions. We take z 0 = [1, • • • , 8] km and t 0 = [2, • • • , 5.5] years. In addition, we excluded observations made after the 5-th year. Therefore, we used a set of N obs = 6112 observations. Fig. 17 represents the 2-dimensional map of the variance of the estimation error. We note that the variance of the estimation error is minimal when there are observations in the neighborhood, whereas it increases as we move away from the observations. This is an expected result that is consistent with what has been observed in the case of time kriging. Since the interpretation of kriging surfaces is complex, we represent these results at constant altitude (cf. Fig. 18). On constant height sections, the estimation structure seems more complex than in the 1-dimensional case, it can be explained by the interaction with observations at altitudes above and below the considered altitude. Moreover, we note that the model tends to the mean of the model when we move away from observations. This is especially true in the case of long-term prediction around the 5-th year.

Conclusion and perspectives

The interpolation, in space, and prediction, in time, of the cloud microphysics in medium and long term are of major importance in weather and climate analysis. Since a perfect estimation is obviously unattainable, it is an issue of uncertainty quantification.

In this original work, we develop a statistical spatio-temporal kriging-based approach that is able to interpolate/predict from the dataset and provide uncertainties. Beforehand, it requires in particular estimating the covariance model parameters; it is performed in a Bayesian setting, which allows for estimation and uncertainty quantification. The approach is then applied to a subset of the CloudSat dataset which shows promising results, especially in the 2-dimensional case where detailed structure appears in the quantity of interest. A natural extension to this approach would be to consider the latitude and longitude variables in order to interpolate horizontally the quantity of interest. Other extensions should consider the joint estimation of Ice Water Content and Liquid Water Content, the estimation of cloud types (cumulus, stratus, etc.) or some nonlinear functions of our actual quantity of interest (i.e., some parametric model that depends directly on Liquid Water Content or Ice Water Content like optical properties). The treatment of these problems would require to resort to suitable kriging estimators (co-kriging, disjunctive kriging, etc.), or, more likely, some adapted versions of them. Finally, the impact of the parameter uncertainty in the kriging results could be rigorously handle by developing a complete Bayesian hierarchical model.

2

  Data and Mathematical Description 2.1 The CloudSat data and the CPR instrument CloudSat has been flying in formation in the A-train with other satellites including Aqua, Aura and CALIPSO. CloudSat payload consists of a 94 GHz cloud-profiling radar (CPR) that was specifically designed to sense cloud-sized particles (i.e., cloud ice, snow, cloud droplets and light rain). It was declared operational on June 2-nd 2006 and has been flying in the A-train until February 22-nd 2018. It follows a Sunsynchronous orbit with an approximately 1330-LT equatorial crossing time. Since 2011 and a battery malfunction, it provides observations only during daytime. The satellite visits the same position of the globe after a period of 16 days corresponding to 233 orbits. Each orbit is achieved in about 1 hour and 58 minutes. The CloudSat radar samples profiles at 625 kHz and has an along-track velocity of approximately 7 km/s, which corresponds to a profile measured every 0.16 second with an along-track displacement of approximately 1.1 km. Each profile has 125 vertical bins of 240 m thickness.

Figure 1 .

 1 Figure 1. CloudSat ground track representing 14 orbits (a), zoom in on the European zone under study (b).

Figure 2 .

 2 Figure 2. Sample of CloudSat cloud classification over a period of 15 days. Clouds are classified according to Clear (no clouds), Cs (Cumulus), As (Altostratus), Ac (Altocumulus), St (Stratus), Sc (Stratocumulus), Cu (Cirrus), Ns (Nimbostratus), DC (Deep Convective).

  a 0 ), and the corresponding estimator is the so-called Minimum Mean Square Error estimator (MMSE). The value a

Figure 3 .

 3 Figure 3. Scatter plot of the altostratus LWC observations.

Figure 4 .

 4 Figure 4. Empirical covariance: (a) 2-dimensional covariance (δt, δz), (b) slice at δz = 0 km and (c) slice at δt = 0 year.

Figure 5 .

 5 Figure 5. Empirical periodogram at (a) z = 6.71 km and (b) z = 6.00 km.

Figure 6 .

 6 Figure 6. Evolution of CLP during the optimization process.

Figure 7 .Figure 8 .

 78 Figure 7. Evolution of each parameters along the optimization process. The x-axes represents the iteration index while the y-axes represent the current value for each parameter. Parameters associated to the covariance of the quantity of interest are represented in red (•); parameters associated to the mean of the quantity of interest are represented in blue (•) and the noise variance is represented in green (•).

Figure 10 .

 10 Figure 10. Comparison of the estimated parameters for 9 geographic area. Parameters associated to the covariance of the quantity of interest are represented in red (•); parameters associated to the mean of the quantity of interest are represented in blue (•) and the noise variance is represented in green (•).

Figure 11 .

 11 Figure 11. (CLP -µCLP) as a function of the correlation time lt around θ for the geographic zones 1 to 3 (a), 3 to 4 (b), 6 to 8 (c) and the 9-th (d). The results for the 9 geographic areas have been grouped in 4 sketches according to their similarities.

Figure 12 .

 12 Figure 12. Kriging based interpolation over a 9 year period. The training set observations are represented as (•), the kriging estimate as (-), and variance of the estimation error (shaded gray).

Figure 13 .

 13 Figure 13. (a) Kriging based interpolation over a period of 36 days centered at the beginning of 3-rd in the dataset. The mean is represented in as a dashed black line (--). (b) Kriging based prediction over a period of ∼ 36 days centered at the beginning of 3-rd in the dataset. The vertical black line (|) represents the correlation time lt after the last observation. In (a) and (b), the training set observations are represented as (•), the kriging estimate as black line (-), and the variance of the estimation error in shaded gray.

Figure 14 .Figure 15 .

 1415 Figure 14. Kriging results collocated to the observations. (a) Observations ψ (•), estimate φ0 (•), estimated mean mΦ(-) and variance of the estimation error in shaded gray. (b) Histogramms of ( φ0 -mΦ) and (ψ -mΦ). (c) Shaded gray representation of the variance of the estimation error and the variance of the quantity of interest Φ.

Figure 16 .

 16 Figure 16. Evolution of each parameters during the optimization process in 2-dimensional case (time-altitude). Parameters associated to the covariance of the quantity of interest are represented in red (•); parameters associated to the mean of the quantity of interest are represented in blue (•) and the noise variance is represented in green (•).

  (a) outlines the situation we have just defined. The obtained kriging surface represented in Fig.17(b) visually fits the observations from Fig.17.Fig. 17(c)

Figure 17 .

 17 Figure 17. CloudSat observations on geographic area 6. The red rectangle ( ) indicates the training dataset, the green rectangle ( ) indicates observations not used in the kriging estimate and the plain black lines represent the spatialtemporal positions where the kriging estimates are computed (a). Result of the 2-dimensional kriging on zone 6 (b) and its associated 2-dimensional variance of the estimation error (c).

Figure 18 .

 18 Figure18. Kriging results at constant altitudes, z 2.17 km (a), z 3.62 km (b), z 5.08 km (c), z 5.96 km (d). We represent the observations ψ (•), the estimations φ0 (•), the mean (-) as well as 2σ around the estimation φ0.

Acknowledgements. The authors would like to thank Audrey Giremus for her contribution to the preliminary step of the work and the two anonymous reviewers for their suggestions and comments.

Data availability. The level 2B-CWC-RO CloudSat data used in this study have been downloaded from the CloudSat Data Processing Center (http://www.cloudsat.cira.colostate.edu/).

the dataset corresponding to the 6-th geographic area.

Discussion on the 1-dimensional model

To conclude, it seems more interesting to carry on with geographic areas associated with higher correlation time, especially those obtained for zones 4 to 7 which corresponds to l t ∼ 1 -6 days, compared to l t ∼ 2 minutes for zones 1 to 3. Considering the exponential covariance model reaches 5% of its total variance at δ t ∼ 3l t [START_REF] Banerjee | Hierarchical Modeling and Analysis for Spatial Data, Second Edition[END_REF], the latter tends towards its 335 mean m Φ after ∼ 6 minutes. Note that the models corresponding to zones 4 to 7 are associated with variances r Φ lower than In order to study the result of kriging prediction, we selected a period that includes observations until a certain date and analyzes the behavior of the prediction beyond the last available observation (cf. Fig. 13(b)). In order to facilitate the interpretation, we have subtracted the estimated mean mΦ from the observations ψ and the estimates φ0 and we have centered the x-axis on the last observation available. This clearly shows a decrease in the kriging result after the last observation towards 0, which indicates that the result tends towards the estimated mean. We plot the value of the correlation time l t as an indication. In the case of the exponential covariance, we get 5% of the variance r Φ for a time ∼ 3l t which corresponds to the time from which the kriging estimation tends towards the estimated mean. It is therefore a result consistent with the chosen covariance model and the estimated correlation time. This means that the prediction will be different from the average when an observation is available in a period of approximately 8.2 days (because l t 2.7 days) before or after the position where the estimation is performed.

Error analysis

In this section, we interpret in more details the kriging results through the analysis of the variance of the estimation error var Φ0 -Φ 0 and the estimated noise. This is accomplished with a kriging estimation at temporal positions t 0 collocated to the positions of the observed data (cf. Fig. 14 (a)). Fig. 14 (b) represents the corresponding histograms for quantities ( φ0 -mΦ ) and (ψ -mΦ ). Although the model used to describe the observed data is Gaussian, the histogram of the zero-mean quantity, ( φ0 -mΦ ), is not exactly Gaussian. It is actually skewed towards zero. However, its range of variation is smaller than the corresponding histogram for the observations, (ψ -mΦ ), which is an expected result as some of the observed variability is attributed to the presence of noise.

The variance of the estimation error var Φ0 -Φ 0 gives an approximation of the variability around the estimated value Φ0 .

Note that it is impossible to compute exactly the difference ( Φ0 -Φ 0 ), but we have estimated the variance of the noise knowing the observations, it is then possible to compute the variance var Φ0 -Φ 0 (cf. Sect. 3). This last variance term can be used to compute posterior realizations of the quantity of interest Φ. Fig. 14 (c) represents the variance var Φ 0 and the variance of the estimation error var Φ0 -Φ 0 . The latter has a lower magnitude than var Φ 0 which is consistent with (11) as the term r t ΨΦ0 R -1 r t ΨΦ0 is positive. Finally, we computed the difference (ψ -φ0 ) (cf. Fig. 15(a)). According to (2), this distribution must correspond to the noise distribution. We represented the normalized histogram of (ψ -φ0 ) in Fig. 15 (b) as well as the pdf of the noise B ∼ N (0, rB ).

We observe that the two distributions have slightly different shapes but similar ranges of variation. The differences observed are mainly due to the fact that the parameters of a Gaussian model are estimated from observations that do not have a strictly Gaussian distribution. Indeed, we notice that the observations have an asymmetric distribution which extends towards higher values of LWC than the Gaussian distribution.

Author contributions. PM and JFG designed the preliminary study, with continuous discussions from all co-authors along the project. JML 445 developed the model code and performed the simulations. JML prepared the manuscript with contributions from all co-authors.

Competing interests. The authors declare that they have no conflict of interest.