Please see the supplement, thanks.

The paper Bayesian inversion of satellite altimetry for Arctic sea ice and snow thickness from René-Bazin et al. describes the application of a novel inversion scheme to simultaneously derive sea ice thickness and snow depth from multi-sensor (Ku/Ka-Band radar altimeter and Ku-Band radar/Laser altimeter) data. One of the novel aspects of the method is that is uses a Bayesian inversion and a data sampling strategy based on random triangulation of surrounding data and thus can work with gappy input data.

A lot of work has undoubtly went into this study and the method has been to my knowledge never applied to snow depth and sea ice thickness estimation. Nevertheless, from seeing the results I am not convinced that this method improves upon the results from much simpler methods.

There are several points leading to my conclusion: Firstly, the chosen forward models include penetration factors. These factors are dependent on snow depth and have a limited capability to describe the freeboard errors of laser and radar altimeters (General Comment 1). Equally important, the penetration factors also seem to introduce multiple solutions, which cannot be resolved with 2-dimensional data vector (General comment 2). The validation is also very limited (General Comment 3) and may not be decisive enough to demonstrate that the CS-IS-2p SIT outperforms the AWI SIT. The CS-IS-3p and CS-IS-4p are definetely not outperforming the AWI SIT. On top of these concerns, I find the description of the methodology could be expanded substantially as key characteristics of the spatial sampling, such as the node and ensemble statistics (General Comment 4) are missing.

I have provided more information in the General Comments below and added specific questions and recommendations in the pdf file attached to this review.

My recommendation for going forward with this would be remove the penetration factors from the forward model (I don’t see how the ambiguity they introduce can be resolved) and to add another sensor sensititive to snow depth (e.g. a simple forward model used in passive microwave retrieval of snow on sea ice). Then the ensemble statistics and residual inversion errors could be used as indication for snow depth uncertainties. Potentially sea ice density could also be inverted for and not prescribed by an external sea ice type mask.

But as the manuscripts is now, I cannot recommend publication at this point.

Best Regards,

Stefan Hendricks

Penetration factors are a central piece of this work. This factor is used to relate a freeboard estimate from any altimeter type to the sea ice freeboard (alpha=1) of snow freeboard (alpha=0), with the difference being the snow depth. Sea ice freeboard is the target variable for Ku-Band radar altimetry, while snow freeboard is the target variable for Ka-Band radar and laser (infrared and green) altimetry. Deviations from the 0 and 1 values are used here and in other literature to describe the freeboard error with respect to the corresponding target variable.

The issue with describing errors of radar-derived freeboard as a direct function of snow depth is that there are many error sources for which snow depth is not a suitable proxy. The authors write in L44+ as their motivation for using the penetration factor that the known backscatter from the snow at Ku-Band raises the backscatter elevation distribution and thus bias the radar freeboard value high. While snow backscatter is undoubtly relevant at Ku-Band, there are issues with the assumption that the freeboard error must be a direct function of snow depth.

The main issue is that sea ice surfaces are very rarely 1D layers at SAR altimeter footprint sizes and different sea ice surface have shown strong differences in backscatter values. Not only elevation distribution of backscatter matters, but also its azimuth distribution. Or more plainly, a patch of thin ice or open water with substantially larger backscatter coefficient than the surrounding snow-covered sea ice a bit off-nadir can dominate the main peak of the waveform and the corresponding retracker range is invariant of snow backscatter or the ice surface elevation distribution.  The backscattered power arrives at later times and biases the retracker range high. And since the off-nadir angle is unknown (except for CryoSat-2 SARin with very limited coverage) the freeboard in these case willbe be biased low, often by factors several times the snow depth thus providing a strong opposing bias to snow backscatter.

Additional issues are sea level anomaly errors directly affect the freeboard errors as well and that also for ice surface backscatter may have a relationship with elevation, resulting in a bias in waveform models akin to the sea state base in ocean altimetry. For gridded freeboard values there might be additional impacts based on the surface type classification and possibly the rejection rates of waveforms.

My main take away message that freeboard errors cannot be realistically modelled only from snow depth. Penetration factors may be used as an empirical bias correction, but they absolutely should not be used to infer the actual backscatter mechanism at local scale as the errors of gridded freeboards.

Inverting for 3 or 4-dimensional model space with only a 2-dimensional data space is an ill-posed problem. But the additional challenge I see is the promotion of false solutions by including the penetration factors. E.g. the observed freeboard differences between the two altimeters can be easily explained by a range of suitable combinations of snow depth and penetrations factor(s). The ice thickness than just can be chosen within the substantial valid range to account for the varying snow mass and to place the freeboards at the correct magnitude. With just two freeboard measurements there is no possibility to resolve this ambiguity.

At least it is my guess that a wrong minimum is the reason why CS-IS-3p and CS-IS-4p do not generate competitive results compared to the other products and not against the CS-IS-2p setup without penetration factors.

One option to reduce the likelihood for false minima would be to define the prior distribution from the climatological values of snow depth and sea ice thickness and a realistic deviation from the mean value. The penetration factor range could also be narrowed to prevent the negative IS penetration factors seen in figures 12 & 14.

But the main take-away message from this analysis is, when taking the results from from CS-IS-2p/3p/4p as presented here is: The inclusion of penetration factors substantially degrades your retrieval. I am not sure that this is the message the authors intend to make.

Airborne data from NASA Operation IceBridge is a good validation data source for snow depth on sea ice but not so much for sea ice thickness. The parameter is derived from the laser scanner and snow radar data, and the estimation of sea ice thickness uses one of the forward model functions. There are more suitable sensors for sea ice thickness validation such as upward looking sonars or EM-induction sounding sensors.

It is also my understanding from the manuscript that the authors use the OIB quick look product as is, and do not apply the sea ice density values for FYI and MYI as in the forward model functions. In this case the OIB sea ice thickness becomes inconsistent, since the OIB sea ice thickness is based on a single value for sea ice density (915kg/m3, Kurtz et al. 2015). An exact match between the inversion result and the OIB thicknesses would than rather point to an error in the freeboards/snow depth input.

I recommend here to use ULS observation, for example BGEP, for sea ice thickness validation. The required data should be available within the author team. That would provide a more robust validation, also in other month than Arctic spring.

Kurtz, N., Studinger, M., Harbeck, J., Onana, V., & Yi, D. (2015). IceBridge L4 Sea Ice Freeboard,Snow Depth, and Thickness (IDCSI4, Version 1) [Data set]. Boulder, Colorado USA. NASA National Snow and Ice Data Center Distributed Active Archive Center.

I find it difficult to assess the maturity of sampling the observation space with the Delauny Triangulation. My understanding is that the number of triangles changes are variable variable, but what are the numbers of nodes/triangles in the ensemble? The only indication of triangle properties is in Figure 1, where the triangles look quite coarse and there seems to be extrapolation to the entire sea ice cover. Some triangles look substantially larger than correlation length scale of sea ice surface types and I am not sure that sea ice in triangles center can be accuratetely described by the node positions.

Another part of my confusion comes from the statement that “the dimension of the model space (the number of nodes n_Hi) will be treated as an unknown variable (L166)” and “At each iteration of the random walk, the algorithm makes a random choice between four perturbations of the model parameters (L196)”. If birth and death of nodes are chosen randomly than my assumption is that the number of nodes does not change substantially in the ensemble. Therefore, does the inversion for spatial resolution means to sacrifice spatial resolution in one area for another?  

My recommendation is to describe the method in more detail and show the results more prominently. How are the nodes distributed and what are interpolation distances in the aggregation to the regular grid? What are the statistics of the ensemble spread?