MCMC sampler fuses satellite altimetry and ensemble simulation to reconstruct Mediterranean sea level variability

Sea level is one of the most consequential variables in climate science, and its slow, interannual swings in a semi-enclosed basin like the Mediterranean Sea are notoriously hard to pin down. Satellites measure it constantly, but models simulate it too, and the two rarely agree perfectly. A team of French researchers at Univ. Grenoble Alpes, CNRS, IRD, Grenoble INP and IGE has now unveiled a computational method that fuses these two worlds with unprecedented efficiency, reconstructing the low-frequency variability of Mediterranean sea level over 26 years while honestly quantifying its own uncertainty. The work, published in Ocean Science, is a technical tour de force that could reshape how oceanographers and climate scientists blend observations with simulations.

The core challenge the researchers tackled is one that has long frustrated the community: climate projections are only as trustworthy as the models behind them, and the skill of those projections depends on whether models can reproduce the long-term, low-frequency variability of the system. Checking that ability against real observations is essential, but doing so in a statistically rigorous way is fiendishly difficult when the model simulation carries systematic errors and the observations are sparse along satellite ground tracks. The new method, developed by Jean-Michel Brankart, Damien Héron, Lisa Weiss, Thierry Penduff and Pierre Brasseur, estimates the low-frequency component of ocean variability directly from native observations, using statistics from a prior long-term ensemble simulation as the starting point.

At the heart of the approach lies a Monte Carlo Markov Chain (MCMC) sampler, a statistical engine that explores a probability distribution by constructing a stochastic process whose attractor is the distribution one wants to sample. The team used a modified variant of the classic Metropolis/Hastings algorithm, one that includes localization of the prior covariance. Localization is a trick borrowed from ensemble data assimilation: because the prior ensemble contains only 30 members, its estimated correlations over long distances are unreliable, so the method drops those small long-range correlations using a Gaussian function of distance and time difference, with a correlation length scale of roughly 100 kilometres and a time scale of about six months. What remains is a local, trustworthy correlation structure that the sampler can exploit.

The estimation problem itself is formulated in a Bayesian framework, where the target is the posterior probability distribution obtained by conditioning the model-based prior distribution on the observations. Crucially, the researchers did not assume the model was perfect. They augmented the state vector with additional free parameters: a time-constant bias field for sea surface height, temperature and salinity; a multiplicative factor to adjust the prior ensemble variance; and another factor to adjust the observation error variance, which is dominated by representativity error from unresolved time scales. Including these auxiliary parameters makes the prior distribution non-Gaussian, which is precisely why linear estimation methods such as standard Kalman filters cannot be applied directly and why the MCMC sampler is needed.

The computational innovation that makes the whole thing feasible at scale is a clever way of generating random perturbations. In a naive implementation, sampling a large-dimension proposal distribution would require computing a linear combination of basis vectors, a cost that grows prohibitively with the size of the state vector. Instead, the researchers exploit the structure of the localized covariance: they generate a moderate-size sample of vectors whose covariance, when raised to the sixth Schur-power, matches the localized correlation matrix. By taking six different vectors from that sample and computing their Schur product, together with one vector from the prior ensemble, they obtain a pseudo-random perturbation direction with exactly the right localized covariance, at the cost of just six element-wise products. With 30 vectors in the localizing sample, the combinatorial effect yields more than ten billion distinct perturbation directions, enough to span all relevant degrees of freedom of the system.

This algorithmic trick is not just elegant; it is fast. The researchers show that the method is well suited to parallel computing with MPI and GPUs, and efficient enough to solve problems with about 10^8 variables and 10^7 observations. In their application, they achieved a speedup of roughly a factor of 100 compared to running on CPUs alone. The overall cost is linear in the size of the state vector and the observation vector, and the authors note that the main practical limitation was the memory required to store the prior ensemble on the GPUs. A classic implementation of the Ensemble Kalman Filter on the same problem would have been substantially more expensive, given the density of observations and the cubic scaling of local updates with the number of observations used.

To demonstrate the method, the team reconstructed the low-frequency variability of Mediterranean sea surface height using along-track altimetric observations from 12 satellite missions between 1993 and 2018, obtained from the Copernicus Marine Environment Monitoring Service. The prior statistics came from a 30-member ensemble simulation of the Mediterranean Sea based on the NEMOMED12 configuration of the NEMO ocean model, run at 1/12-degree resolution over 75 vertical levels and driven by realistic atmospheric forcing between 1979 and 2020. The ensemble dispersion was generated by tiny perturbations of the initial conditions, so that the spread reflects only the intrinsic variability of the nonlinear system. A low-pass Lanczos filter extracted the interannual component, discarding periods shorter than about 2.5 years.

The resulting posterior ensemble estimate was assessed against independent observations using cross-validation, with the altimetric data split into two subsets so that observations not used in the estimation could serve as a check. The rank histogram of the independent observations in the updated ensemble came out essentially flat, indicating that the ensemble solution is globally unbiased and only slightly overdispersive. Self-consistency checks comparing the full-observation experiment with the two partial-observation experiments suggested that the posterior uncertainty might be slightly underestimated, possibly because the statistical assumptions, such as uncorrelated observation errors, give too much confidence in the observations. The posterior ensemble mean also became remarkably similar to the Aviso mapped product, despite the two being obtained through very different processing routes, which the authors interpret as evidence that both approaches made consistent use of the same along-track data.

Beyond the reconstruction itself, the method produced consistent ensemble estimates of the model bias and of the observation error variance. The bias estimate revealed that the model simulation has a too weak cyclonic circulation in the Balearic Sea and the South-Western Algerian basin, and too weak Iera-Petra and Pelops gyres in the Eastern basin, features already noted in earlier work. The observation error standard deviation, dominated by representativity error from unresolved scales, was found to be larger in regions of intense mesoscale activity, such as where eddies detach from the Algerian current. However, the adaptive estimate of the missing prior ensemble variance proved difficult to control by the observations; the various Markov chains spread from their background value without converging towards a common pattern, so the researchers computed the missing variance as a diagnostic instead, revealing that the model largely underestimates its own uncertainty along coasts and in shallow basins.

The broader implications of this work reach well beyond the Mediterranean. The authors argue that producing reliable ensemble reanalyses, forecasts and projections requires that the main sources of modelling uncertainty be adequately represented, and their diagnostic framework provides a concrete way to identify where variance is missing in the model world. They suggest that moving towards a stochastic formulation of the ocean model, as has been proposed for NEMO and recently developed in the CROCO model, would not only improve the solution of modelling and estimation problems but also help understand which observations are most needed to fill the information gaps. With the code and data openly available, the method stands as a significant step forward in combining numerical ensembles with satellite data, and it is likely to influence future improvements in data assimilation for climate and ocean science.

Subject of Research: Ensemble estimation of Mediterranean low-frequency sea level variability using an MCMC sampler combining observed and simulated data

Article Title: Efficient ensemble estimation using an MCMC sampler: a reconstruction of the Mediterranean low-frequency variability combining observed and simulated sea level

Article References: Brankart, J.-M., Héron, D., Weiss, L., Penduff, T., & Brasseur, P. (2026). Efficient ensemble estimation using an MCMC sampler: a reconstruction of the Mediterranean low-frequency variability combining observed and simulated sea level. Ocean Science, 22(5), 3207-3224. https://doi.org/10.5194/os-22-3207-2026

Image Credits: AI Generated

DOI: 10.5194/os-22-3207-2026

Keywords: MCMC sampler, Mediterranean Sea, sea level, satellite altimetry, ensemble simulation, Bayesian inference, data assimilation, low-frequency variability, ocean modelling, uncertainty quantification, GPUs, NEMO ocean model