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

The authors propose an original ensemble-based approach designed to improve the representation of the long-term and low-frequency variability of the climate system (ocean compart)by providing an appropriate probabilistic initial condition for carrying out climate projections.

The authors formulate a Bayesian inverse problem in which the main unknown is the Sea Surface Height (SSH) in the Mediterranean Sea (MED) and where the objective is to reconstruct the interannual SSH variability between 1995 and 2019 by using along-track satellite altimetry observations to constrain the problem. For that purpose by using the NEMO hydrodynamic model, the authors construct a prior ensemble of simulations ( 30-members ) by perturbing the only initial condition, after a one-member spin-up period to asses the prior probability distribution (pdf) of the inverse problem. The prior pdf is mainly described by the ensemble mean and error covariance matrix. To reduce the sampling errors related to the limited size of the ensemble, the correlation of the prior ensemble members are localized in space and time using a Schur product involving Gaussian kernels. The authors introduce additional statistical parameters in the control vector. These parameters are also identified along with the low-frequency SSH time evolution. The existence of a prior ensemble bias motivates the authors to also control stationary in time biases in the Dynamic Height, Temperature and Salinity fields. Also, the spread of the prior ensemble being insufficient to explain the misfit with observations, the authors decide to also control the error variance both on the prior SSH estimate and the satellite altimetry observations. Due to these auxiliary statistical parameters introduced to control the bias and the spread, the overall prior probability distribution becomes non-gaussian. According to the Bayesian framework, the authors combine the prior probability distribution with the likelihood based on the altimetry observations and derive a cost function which is non-quadratic since the problem is non-linear. They propose to apply a Markov Chain Monte Carlo solver (MCMC) to sample the posterior probability distribution, rejecting or accepting the Markov Chain candidate iterate based on the increase or decrease of the cost function. This MCMC method applies a variant of the Metropolis/Hastings algorithm : the drawing of the random direction of perturbation to construct the Markov Chain candidate iterate involves the use of Schur products designed in such a way that the covariance of the pseudo-random direction of perturbation is the covariance of the localized prior ensemble. With such a variant, the MCMC rejection/acceptance of the candidate iterate is only driven by the likelihood function. 

At the end, the method enables to obtain an updated ensemble with members sampling the posterior probability function with a posterior estimate of the Dynamic Height evolution in the MED sea. The authors discuss their results and propose different approaches to validate the Dynamic Height reconstruction. 

First, they compare the interannual variability of the reconstructed Dynamic Height maps to the AVISO L4 product : compared to the prior ensemble, the spread of the posterior ensemble appears to be reduced at time period and location where dense along-track altimetry measurements are available. The authors also propose a cross-validation framework using part of the observations to constrain the likelihood and keeping the other part for validation. In this manner, they can assess the reliability of the updated ensemble by using rank-histogram with independent observations. They discuss the reconstructed Dynamic Height bias as well as the updated ensemble spread and tuned observation error variance.

The article is very well written. The authors provide many details so that the reader can understand the chosen methodology. The authors tackle a complicated estimation problem with non-gaussian distributions, biased datasets and ensemble spread issues. Outside the Gaussian and linear context, solving inverse problems is a challenging task. They apply an original variant of the MCMC Metropolis/Hastings algorithm which enables to optimize the MCMC solver, re-using the ensemble localization operator. They test this algorithm on a realistic application, which contributes to methodological progress.

I recommend the publication of this article. I have only few minor comments/questions that are listed below.

Minor comments/questions :

• L30 simulationsimulationss -> simulations
• Eqn 8 : gamma_loc depends on r and delta_t, for the sake of clarity I suggest you had the delta_t dependence on the LHS of the equation.
• L183 + Eqn 5, L215 + Eqn 10 : you decide to tune the prior ensemble error variance (Pr,1 matrix with the alpha diagonal matrix) because the spread of the ensemble is not sufficient to explain the misfit with altimetry observations. Why also tuning the observation error variance (Rr matrix with the beta diagonal matrix) ? You explain that it should account for the observation representativity error.  Could the tuning of both alpha and beta vectors lead to introduce unnecessary degrees of freedom in the inverse problem ? Can this affect the MCMC convergence ? Regarding the alpha parameter (spread) ?
• Eqn 11 : in the quest for accuracy, I would suggest adding a constant value to the LHS of the cost function observation term, so that equality between the LHS and RHS of equation 11 holds, even if at the "minimization" stage this extra constant term won't play any role.
• In the discussion, you did not comment on the Temperature or Salinity biases reconstruction. Are those two biases consistent with the Dynamic Height bias ?
• To produce the prior ensemble, can you explain further why you choose to only perturb the initial condition ?
• Figure 5 : for the reliability diagrams, I would suggest using the same scale for the Y-axis of the top pannel plots on one hand, and also for bottom pannel plots on the other hand. This way we would clearly see (top pannel) that SSH-B experiment assimilates more observations than the SSH-A experiment. For the bottom pannel we expect the same amount of verifying observations for the two rank histograms. Could this help to compare the distorsion of the reliability in a more quantitative manner  ?