| 16 Jun 2026
Efficient ensemble estimation using an MCMC sampler: a reconstruction of the Mediterranean low-frequency variability combining observed and simulated sea level
Abstract. The skill of climate projections depends on the ability of models to reproduce the long-term and low-frequency variability of the system. It is thus important that low-frequency model statistics can be checked against observations. In this paper, a method is proposed to estimate directly the low-frequency component of the ocean variability from native observations using statistics from a prior long-term simulation. It is designed to account for possible model biases and to provide an estimate of the correction required to fit observations. The result is obtained by an MCMC sampler (modified to include localization of the model covariance), which provides an ensemble description of the solution, so that uncertainties can be properly assessed using independent data. This algorithm is shown well suited to work with MPI and GPUs, and efficient enough to solve large-size problems (about 108 variables and 107 observations). The approach is illustrated by the reconstruction of the low-frequency variability of the Mediterranean sea level, using statistics from a 1/12° resolution ensemble model simulation. The resulting ensemble is assessed against independent observations (by cross-validation), showing good reliability (flat rank histogram). The method also produces a consistent estimate of the model bias and of the observation error variance (mainly representativity error), while the missing prior ensemble variance is shown to be less controlable by the observations, and thus rather computed as a diagnostic. Overall, this application shows the importance of reliable model statistics, and thus the importance of enhancing model simulations to represent all main sources of uncertainty.
The authors describe a novel method to deduce sea level over the Mediterranean Sea from altimetry observations and from the model results, accounting for model deficiencies in terms of bias and underestimated ensemble spread.
I think that the authors provide important results and the application here will in future lead to future improvement in the context of data assimilation. The problem that the authors address belongs to the wide class of inverse problems (including for instance also data assimilation).
In this context, virtually all approaches assume that the error covariance is perfectly known, bias free and Gaussian distributed. Which leads to a quadratic cost function which can be minimized in closed-form by solving a linear system. The paper illustrates the issue very well that when the error covariances and biases are themselves considered unknowns, closed-form solutions (e.g. Kalman Filter) are no longer available, the problem must be solved using more general techniques like the MCMC solver.
I think this paper is very well suited for Ocean Science because it includes significant theoretical advances but at the same time provides a concrete application for SSH estimation from observation and model results. They also include a description of the MCMC algorithm in the appendix to make their approach easily understandable for a wide audience.
The authors also show limitations of the initial procedure (difficulty for the convergence of the adjustment factors) and the use of additional constraints is discussed and tested.
One might question if the application of this method has not been done on a (borderline unrealistic?) hard test case where the model simulation did not include any perturbation from the forcing fields and boundary conditions and thus significantly underestimate the error variance. The method had thus to rely significantly on the adjustment of the priori error covariance. Apparently the model ensemble was not created with the purpose of providing a realistic error variance. However, if the method works with a significantly underestimated error variance it will also work for a moderate underestimation when the underlying ensemble is more realistic.
This manuscript is very well written and any question that I had while reading it was answered in the subsequent paragraph. Therefore I recommend publication with minor clarification listed here below.
Minor comment:
Line 107: ...subsequently driven by the same surface and lateral fluxes throughout the 39-year integration period (from January 1979 to December 2020). Ensemble dispersion is generated by tiny perturbations of the initial condition (as explained by Héron et al., 2026), so that the spread of the ensemble results only from the intrinsic variability of the nonlinear system
Can you comment on the lack of variability due to uncertainty in the surface fluxes?
Line 196: For future studies, could the parameter l_3, l_loc , tau_loc also be a random variable to reflect their uncertainty?
Equation 11: it seems that log and ln are both the natural logarithm. I would suggest making the notation uniform.
Line 272: For instance, with m = 30 and 6 Schur products, the number of different directions of perturbation is: N = 30×29×28×27×26×25×24 ≃ 1.026×10^10 , which is sufficient to span all relevant degrees of freedom of our system.
The fact that here m=30 is it related to choice of the number of ensemble members which is also 30r? It seems that these two are unrelated, or not? The number of directions is substantial, but their usefulness also depends on how the initial perturbations are chosen? Can you give more details on this? Would a rough lower bound be the area (or volume) of the domain divided by the effective support of the localization function, adjusted by the sqrt(6) parameter? (by effective support I mean, loosely speaking, the area where localization function is significantly different from nonzero)
Line 305: From Fig. fig:map:posterior…: The figure number is missing