The authors present a novel snow microwave radiative transfer model for SAR altimetry signals including surface scattering and snow-volume scattering. Overall, the manuscript is well written but sometimes lacks some clarity. Especially, the volume and interface backscattering function implementation on DDMs was hard to understand from the equations and text presented. Therefore, I believe the manuscript requires major revision.

Author’s Response (AR): We would like to thank the reviewer for their comments and the chance to improve the accuracy of our text. 

We will acknowledge reviewer’s contribution in the manuscript.

**** RC: Page 2 line 31: I personally would write cm/yr instead of cmyr^{-1}

AR: We have changed this unit as proposed.

**** RC: Page 3 line 59: I am not sure about the formulation “an inverse method”. I guess you mean a nonlinear optimization approach fitting a measured waveform with a modelled waveform maximizing the likelihood?

AR: Inverse method is common but improper, we have replaced by “2) a nonlinear optimization approach fitting the measured waveform with the modelled one ”

**** RC: Table 1: I guess AN means analytical (okay it is probably analytical numerical based on table 2, but then I do not understand why Ray is AN and Dinardo is A)? The antenna pattern in Ray et al. (2015) was really a free function? I thought it was Ellip. Gaussian as well. Okay in the paper he did not seem to define the antenna pattern to be Gaussian, but in the SAMOSA based retrackers I am very sure that the Ellip. Gaussian approximation is used. The surface backscatter in SAMOSA should be Gaussian as well? Please check the Halimi retracker as well for the surface backscatter

AR: We have replaced AN by “Mix”.

We confirm that Ray et al. 2015 present their model with the f0 and f1 functions without providing the analytical form in case of a Gaussian antenna pattern. For this reason we classified it as AN (now Mix) and it is implemented with a numerical integration in SMRT. To our understanding Dinardo et al. 2018 is the first publication to provide the analytical form.

Despite the interest, the retrackers derived from the published models are not reviewed in our paper because it is more difficult to track technical notes than published papers. To our opinion, table 1-3 makes it sufficiently clear that the information relates to the paper, in the header, while the name of the model / retracker is only as the last row. We do not see how to improve this distinction without overloading the paper.

**** RC: Table 2: The caption should be changed to match with figure 1. Figure 1 is early models and Figure 2 is models. Maybe say something like newer or recent approaches? Or make one table for ocean and one for ice-sheets and one for sea-ice and soil (not sure to be honest)?

AR: We had to split the table mainly for formatting reasons, but we tried to find slightly different captions. We have added “more recent” in Table 2 caption as suggested.

**** RC: Buchhaupt (2018) is a bit outdated but fine to use. Maybe add in the text that an update for stack retracking (10.3390/rs15174206, 10.1016/j.asr.2022.12.034 ) and an antenna pattern update of that approach ( 10.1016/j.asr.2025.02.056 ).

AR: We have added some of these references in two places:

In the Speckle section, we have added:

“When averaging the delay-Doppler map to compute the waveform, this noise distribution approximates to a gamma distribution \citep{halimi_2014} under stationarity approximation, which, with a large number of Doppler beams, tends to resemble a Gaussian distribution \citep{wingham_2004}.

However this approximation is in general not valid for a moving target, as on the ocean \citep{buchhaupt_2023} or for a complex terrain, as on the ice-sheet \citep{aublanc_2025}. The averaging process (stacking) is more complex and advanced strategies can be deployed. This is not considered in the new development.”

and

“The Gaussian assumption remains convenient because it is well-suited for analytical calculations. For this reason, complex antenna patterns can be described as a sum of Gaussian assumption as in \citet{buchhaupt_2025} an advanced version of \citet{buchhaupt_2018}.”

**** RC: The missing yes or no in terrain slope means that e.g. Boy and Dinardo don’t consider the terrain slope?

AR: It is a subtle nuance between “no” = the slope was mentioned as neglected or the model makes assumptions that prevent to consider slope without modification, and the intermediate case where the the slope is not addressed in the publication but the model could take it into account easily (e.g. when a DEM is used as input). We have added “-” and explained this nuance in the caption of table 1.

**** RC: The PTRs in Buchhaupt 2018 is btw arbitrary the paper only uses the sinc^2 PTR as an example.

AR: We have updated the table accordingly

**** RC: Table 3: Jack Landy is a coauthor, so I guess he checked it, but wasn’t the facet based retracker allowing an arbitrary antenna pattern and PTR or do you refer to the final LARM implementation here? If yes, I think that pitch and roll only had a limited support which might be worth mentioning? Maybe check with Jack.

AR: We have checked and corrected the tables.

**** RC: Page 9 line 178 needs some clarification as I do not understand why they do not account for slope variations. For open ocean and sea-ice those are considered in the backscatter function as far as I know. Or do you mean non random variations?

AR: Yes, we meant “deterministic” slope variations, but it is better to remove this statement, it is too imprecise.

The new text is “Consequently, these models cannot simulate heterogeneous terrain with spatial variations in roughness within the footprint; the roughness is assumed to be uniform and Z is a true random variable rather than a random field Z(x, y). “

**** RC: Page 9 line 180: The ocean PDF for SAR altimetry is actually 2D (one dimension is the elevation and one dimension the vertical velocity) and the 1D PDF is an approximation yielding SWH errors for open ocean surfaces. However, for sea-ice it is indeed 1D. Not sure if you want to discuss this here of not since your focus is sea-ice.

AR: We appreciate that for the ocean point of view, these statements are not precise enough, but we prefer not to address this complexity given the cryosphere focus of the paper.

**** RC: Line 185: The height and slope are correlated for non-Gaussian surfaces? Please provide a reference for this claim.

AR: Our claim is not correct, thank you for spotting this logic problem. It is correct that the height and slope are independent for Gaussian surfaces (Eq 5.5.6 in Adler and Taylor, 2009) but we can not conclude for non-Gaussian surfaces.

“This is an important consideration for any non-Gaussian target, as it affects the expected surface scattering response. For Gaussian surfaces, the local height and slope are independent \citep[eq 5.5.6 in ][]{adler_2009}, so the surface topography effects can be included by convolving the scattering response with the elevation distribution of scatterers. For non-Gaussian stationary surfaces, the height and slope may not be independent,”

**** RC: Section 3.1: The DDM and \sigma^0 terms need a clearer definition, preferably with equations defining each. Otherwise, it is unfortunately very hard to follow the authors intention in this section

**** RC: Eq. 3: Is \sigma^0_{interface i} a function of the incident angle? In the text below q. 3 it is not. It would also be helpful to define \sigma^0_{volume} and DDM_{volume}. In  my (still unpublished) work I just assumed DDM_volume == DDM \conv \sigma^0_{volume}(\tau)

AR: For these two comments, we rewrote entirely Section 3.1 to introduce more formally the equation computed in step 3. For this we start by explicitly assuming first order scattering, write the DDM as a sum of the contributions from each depth z in the snowpack (as an integral) and from the interfaces (as a discrete sum). Then, assuming that volume backscatter is virtually independent of the zenith angle, it is possible to rewrite the volume term as a linear convolution, as in Wingham et al. 2004 (eq 17) and as suggested by the reviewer. For the interfaces, this latter approximation is not possible because the backscattering coefficient depends on the zenith angle strongly (and every interface is potentially different due to roughness and dielectric contrast), so that the DDM computation is more complex than a linear convolution.

The remainder of the section (description of the steps 1, 2, 3) is adapted to avoid duplication, and to use the symbols now introduced in the beginning. In the new version we avoid the terms volume DDM and interface DDM.

**** RC: Line 517: A horizontal correlation length of 10 cm sounds very short. Are you sure that this is not a typo?

AR: Tens of centimeters is a correct order of magnitude for the snow surface on the ice sheets.

We have added the citation: Lacroix P, Legrésy B, Langley K, et al. In situ measurements of snow surface roughness using a laser profiler. Journal of Glaciology. 2008;54(187):753-762. doi:10.3189/002214308786570863

**** RC: Line 557: A more detailed description of your roughness scales would be helpful if I did not overread it. Is the roughness w.r.t. the air-snow interface?

AR: L 478 gives the definition and role of the two roughness scales considered here.

To clarify this important distinction, we decided to use “radar-scale roughness” everywhere it relates to roughness that influences the radar backscatter intensity (formerly “small-scale roughness”) and use “large-scale roughness” for the height distribution (topography) that influences the spread of the range of the echo.

L478 we modified the text: “Key measurements for altimetry are the characteristics of the surface roughness at the radar wavelength scale and the large-scale height distribution within the footprint. Both are related to height variations of the surface, but differ in the spatial scales. The former (hereinafter called radar-scale roughness) controls the backscatter intensity of the surface (and the interlayer interfaces), while the latter (hereinafter large-scale roughness) controls the spread of the delay of the echo and influences the local incidence angle. “

We propagated and uniformized this terminology throughout the paper.

We also clarified that the surface roughness parameter (MSS) is used for the internal interface as well.

**** RC: Figure 6: What extinction coefficient would that correspond to?

AR: The extinction coefficient is 0.10 m^-1 corresponding to a significant penetration of several meters. We have added this information in the caption.

Dear Chief Editor,

Multiple iterations with Copernicus staff have been done to conform the Code and Data policy between the submission and the publication in preprint.

The “code and data availibility” section, page 35, indicates:

- The code published under the license LGPL-3.0-or-later and the documentation are available from the archive at https://doi.org/10.5281/ZENODO.17808241 . The link is working.

- The in-situ measurements are available at: https://doi.org/10.18709/PERSCIDO.2022.05.DS367. As of today, the link is temporarly not working, but the manuscript provides an alternative link for an easier access, and this link is working. In addition, I provide here the zenodo link to this repository is https://doi.org/10.5281/zenodo.6519037

Regarding the BC-005 ESA Land Ice Thematic Products data, they are publicly available from ESA. The extracted data are now copied to the archive: https://doi.org/10.5281/zenodo.19231174

**** Reviewer Comment (RC): This is a great paper that tries to provide a simulation tool for the altimetry DDMs. I have several comments in the attached pdf.

Author’s Response (AR): We acknowledge reviewer’s positive contribution to our work, and will add a statement in the manuscript accordingly.

**** RC: This paper aims to provide a Delay Doppler Map simulation tool chain for the cryosphere, especially for the radar altimeter observations over the ice sheets. The simulator not only wishes to include the effects of the layering ice sheet structures, but also the elevation changes in the resolution cells. Wave statistics are also considered such that the coherent and incoherent waves are included in the model. This is indeed an ambitious work, as the effects from those many factors are very hard to be evaluated correctly. The simulation majorly relies on the incoherent model simulations. Coherent effects are included but do not seem to be evaluated in a rigorous way. This is however acceptable given the limited ability to perform full wave simulations to Ku band density media. Overall, the paper provides a valuable tool for high frequency microwave altimeter DDM simulations. There are several details that need to be clarified.

**** RC:1. The authors reviewed several models that calculate the altimeter wave forms. It will be great if the authors can provide some math equations in the main body or in the appendix. This would help understanding the simulation algorithm.

AR: Regarding the model review, we have considered this question from the beginning and decided not to include equations for the following reason. The parts common to the eight delay-Doppler models are small, and despite they all start from the classical radar equation, their formulation diverges quickly depending on assumptions and choices. In these conditions it is not possible to present a few overarching equations that would help readers understand how they work and differ. Even the general equations 1 and 2 presented in our paper do not strictly apply to all the models.

On the other hand, providing all the equations as given in each paper would require introducing numerous symbols and miss the clarification objective. Instead recasting all the equations of each model in a common symbol system would be useful but this is a very significant task. This could be the topic of a dedicated review paper. Moreover, it is not how we implemented these models in SMRT, as each model is implemented as close as possible to the original paper, with references to the original equations at every line (L415-416). Our approach has been to describe qualitatively the differences in the Background section and let the readers and SMRT users refer to the original papers for the equations.

Related to this comment, the other reviewer asked more formalisms for a better description of the algorithm. For this, we rewrote section 3.1 with mathematical symbols for the scattering, extinction coefficients and effective permittivity as well as the general equations to combine the backscatter coefficients and the delay Doppler map. We have also added the effective permittivity symbol in the workflow diagram (Figure 2).

**** RC: 2. In the implementation of volume scattering in the DDM, how is the delay for the snow to be calculated? Or alternatively speaking, what is the speed of light in the ice sheet when calculating the delay?

AR: It was mentioned L 421 but we added in the new version of the Section 3 more explicit information: “In addition we compute the delay and attenuation at depth $z$ relative to the surface, represented by the generalized function $\Delta(z, \tau)$ using the extinction $K_e(z)$, transmittivity of the interfaces and the wave speed deduced from the refractive index $\sqrt{\epsilon_{eff}(z)}$”

**** RC: 3. Coherent effects between boundaries. No doubt there will be multiple scattering effects between the internal boundaries. Please clarify in the main body whether this effect is considered in the waveform and DDM simulations and also tell why.

AR: In the review L 148, it is written “They adopt a geometrical approach to energy transport (radiative transfer) and only account for first-order scattering (i.e. single reflection at the surface or in the volume), which is certainly sufficient when the surface scattering dominates”.

In L397, about our implementation, we mention that first order scattering is assumed: “This relative decoupling is made possible first because we account for only first-order scattering”.

We propose to further add this assumption around L434:

“The third step combines the backscatter and DDM as follows:

\begin{equation}...

\end{equation}

neglecting multiple scattering between the layers and multiple reflections between the interfaces.”

**** RC: 4. Rough surface effect. Maybe I missed it, but could the authors tell which kind of topography is used? The random modelled one or the deterministic one form DEM? For the DEM, what is the size of the DEM “patch”?

AR: It was not explicit. In the section “comparison of DDM models”, we propose to add the following text (in bold):

“Figure \ref{fig_waveform_comparison} shows the waveforms calculated with all the delay-Doppler models for the Sentinel 3 Ku-band altimeter parameters for a Gaussian surface with a RMS height of 40\,\unit{cm}, without any volume underneath, the satellite perfectly at nadir with a circular antenna, no apodisation and the delay window widening of 2 (i.e. 256 gates are used before slant range correction). These simplified conditions are chosen to use only the common capabilities across all DDM models. For the models working with a deterministic DEM as input (Boy17, Landy19), independent samples of the Gaussian distribution were drawn for every 10\,m$\times 10\,$ cell in the footprint, thus allowing comparable behavior to the other statistics-based models. “

**** RC: 5. This question is an open discussion; I just want to know your opinion. In the paper, the rough surface is decomposed into microwave roughness and DEM roughness. Based on what could this partition be made? Partitioning the scales is a rather difficult problem which now is determined empirically just like the two-scale ocean surface scattering models.

AR: We share the same opinion, it is a difficult problem, and unfortunately given the importance of the roughness in the results, this problem has some consequences.

The most problematic is to find the cut-off of the radar-scale roughness in the case roughness parameters are estimated from observations (from laserscan or photogrammetry) as in our results in Figure 9. In Larue at al. 2021, the original work, from which we got our parameters, the process of computation is described in details. The cut-off is simply the size of the original cloud obtained by photogrammetry in the field (~10×10 m2 for ASUMA sites and ~ 5×5 m2 for EAIIST). Only the general linear trend is removed from the cloud, the residuals are considered as the radar-scale roughness. Choosing smaller clouds would reduce the RMS height and affect the correlength length, with probable impact on the MSS. Despite this relatively arbitrary cut-off choice, the results presented in Figure 9 show that the order of magnitude is fair.

Better results are obtained after optimizing the radar-scale roughness parameters, as expected. In this case, the cut-off scale is implicit and there is virtually no consequence on the results, because the optimization hides the problem.

**** RC: 6. Figure 6 applied slant range corrections to the DDM. Could the authors elaborate more on how the correction is performed? Also please provide a color bar to the map, it helps with the evaluation of the relative strength of each DDM bin

AR: We propose to add a figure in Section 3.1 to illustrate the principle of the slant range correction. We also explain how

In this stage 2, the slant range correction (a.k.a delay migration), the critical processing step providing the advantage of the SAR over LRM mode, is part of the analytical formulation from some models, yielding directly the corrected DDM. However, in other cases (Halimi14, Wingham04, Boy17, Landy19), the model first computes the uncorrected DDM, and we then migrate each cell in delay, with a left translation as depicted Figure \ref{fig_slant_range_correction}. The delay depends on the square of the Doppler frequency $f$ as $\tau_0(f)=

f^2 \alpha * H * \lambda^2 / (8 * V^2)$ \citep[(eq 7 in ][]{raney_1998}) with $\alpha$ the Earth curvature correction, $H$ and $V$ the satellite altitude and velocity, and $\lambda$ the radar wavelength. This ensures that, in the slant range corrected DDM, the echoes for all the Doppler frequencies are aligned at the range of the zero-Doppler echo (dark green line) by compensating for the additional range sensed at different along-track distances from the target on ground (dark red parabola).

The color bar has been added as suggested in Figures 6 and 7 (now 7 and 8).

**** RC: 7. I would like to recommend the authors provide a DDM area map, maybe in the appendix that shows the area each DDM bin corresponding to. Since the simulation is at Ku band, I would expect the difference in intensity is majorly due to the difference in the area that each DDM bin covers. Also, providing the delay and doppler weighting map on the ground surface would also help interpret the results.

AC: To our understanding our Figure 1 serves this purpose. It shows how each location in the footprint contributes to the DDM map. The equi-delay is concentric and the equi-doppler is rectangles along the x axis.

**** RC:8. Also on figure 6, could the authors explain why the simulated wave form using SAR mode is sharper than the LRM simulation?

AR: As noted in the text, this sharper waveform is typical of SAR mode, and it is one of the main inherent motivations why Raney 1998 proposed this new technique. To avoid repetition, we now point the reader to this reference:

“The simulated waveforms have a typical shape of UF-SAR altimeter waveform on a flat surface \citep{raney_1998, aublanc_2018} with a sharp rise (leading edge) and quasi-exponential decrease (trailing edge). This decrease is quicker than in LRM waveforms \citep{larue_2021} because the surface area effectively contributing to each time gate (Fig. \ref{fig_impulse_function}) is decreasing with time instead of being constant for the LRM as explained in \citet{raney_1998}.”

This information is briefly addressed in the text “The simulated waveforms have a typical shape of UF-SAR altimeter waveform on a flat surface \citep{raney_1998, aublanc_2018} with a sharp rise (leading edge) and quasi-exponential decrease (trailing edge)”. The original SAR altimetry paper by Raney 1998 describes the difference between LRM and SAR altimetry.

**** RC:9. Section 5.3 delay window widening. Could you elaborate more on why extending the gates from 128 to 256 makes the DDM wider in the doppler domain? My understanding is that using longer gates only allows longer delay signals to enter the DDM but it should not provide energy to small delay bins

AR: Before correction, it is right that only longer gates are impacted by cutting the acquisition at 128 gates. However, the slant range correction consists to bring high doppler (positive or negative) with long delays to small delays. This is why some “long delays” are converted to “small delays”.