=== Summary and general comments ==

This is a well-organized and well written paper describing the application of physics-informed neural networks (PINNS) for further improving the generation of sub-glacial bed topography datasets, building on previous “BedMachine” efforts that have been ongoing for the past decade or so. Two conservation equations – the continuity equation and a momentum balance equation for ice flow – are considered and introduced as additional constraints on the loss function, akin to their introduction as constraints in the “cost function” of PDE-constrained optimization (a good analogy to consider including for readers more familiar with the language of glaciology modeling / optimization?). The two constraints are considered on their own and in combination and the resulting bed topography datasets are compared and contrasted with previous BedMachine results for three glaciologically distinct regions. Overall, the authors argue convincingly that the new approach has merit and demonstrates potential for improving inferred bed topography in regions where the traditional BedMachine approach begins to break down.

Overall, the work seems very worthy of publication and readers of The Cryosphere will find it a worthy contribution. My suggestion would be to accept for publication with minor revisions, noting that these are those suggested revisions (detailed below and identified by their line number in the submitted version) are largely editorial in nature.

My one more substantial suggestion – not necessarily for this publication but possibly for a future effort – is that I think it would be very useful to redo this exercise for a single region (I realize the computational cost could be a challenge, so pick a single region, like the most challenging one discussed here) but using L1L2 (Blatter/Pattyn) for the momentum balance model as opposed to SSA (which you’ve already done and could reuse as the baseline). Because the former allows for internal deformation, a possibly very different vertical velocity profile (and hence depth-averaged velocity and flux divergence) might be implied in regions of slower moving ice (noting that the modeled 2d surface velocity field could still be used as the velocity constraint, so that there should not necessarily be any significant reformulation of the loss functions discussed herein). It would be interesting to see if the more accurate stress balance constraint helped to alleviate any of the remaining problems discussed below.  

=== Detailed comments ==

21: Are these refs now considered the definitive source for defining the amount of potential SLR locked up in the ice sheets? If not, maybe consider adding one or two more from other authors for the sake of diversity?

23: “numerical ice sheet modeling”

25: It seems like a summary-level reference might also be appropriate here (?), e.g. something from one of the recent IPCC reports (that integrates results from a large number of individual publications).

28-29: It should probably be noted here that these experiments were assuming a marine ice sheet with a significant over-deepening inland (since this configuration would necessarily be more sensitive than say, an ice sheet grounded above sea level).

37-38: Is the ~2km limit proposed here coming from the Durand et al. (2011) paper? While I more-or-less agree with this idea, I don’t know that this single reference is adequate to support the precision implied by this statement. Maybe consider softening it a little bit to something less precise, e.g. “order km-scale spatial resolution”?

63: “three regions in Greenland”; maybe add a few words of clarification here that they are glaciologically distinct / different? E.g., presumably you mean regions where velocity occurs primarily via fast sliding, a region where it occurs via a mix of sliding and deformation, etc.

Figure 1 caption: “The loss function is comprised of …” or “The loss function includes data loss …”

89: “fully connected layers”, maybe use “fully connected (‘dense’) layers …” ?

101: Should it be “the apparent mass balance residual” ?

103-105: I am guessing that maybe this is discussed further below (?), but it seems like you are already potentially limiting the usefulness of this approach by restricting the momentum balance to SSA. I.e., if one of the main interests here is in improving the inference in regions of slower moving ice flow, which is presumably due to less sliding and more internal deformation, then SSA doesn’t seem like the right assumption to make for the model dynamics. I know that ISSM has higher-order approximations available (e.g., L1L2 or “Blatter-Pattyn”). Has that also been explored (acknowledging the obvious additional computational burden) and compared against the approach using SSA?

116: “… from THE regional climate model RACMO…”

Section 2.3: It sounds like the basal mass balance term in equation 2 is assumed to be ~0? If so, it would be good to note that explicitly here in the discussion of the apparent mass balance term.

152: “to prevent from taking” (omit “from”?); “or diving by zero” (“dividing by zero”)

176-177: Would it be worth commenting on the choice of median vs. mean? Is the median chosen because of the small number (5) of samples, such that the mean could be easily biased?

180: By “challenging to implement”, do you mean where the traditional / previous mass conservation approach does not perform well? Implementation sounds more like the approach is challenging, but I imagine the approach is just as easy to implement in these regions, it’s more the prior / baseline result that you are not happy with.

242: It’s not clear here exactly what “Fig.2(1)” is referring to.

Figure 5: In the caption for this and figure 4 it would be helpful to remind the reader which dataset is subtracted from which and shown in panels d-f (e.g., PINN minus original BedMachine product or vice versa?).

Table 4: It might also be useful to provide some percent / fractional metrics here? E.g., for the apparent mass balance RMSE, how does that number compare to the average apparent mass balance over the same area? Such a table could be added to the SI if it’s not deemed important enough for the main text.

3.2.2. – It’s left hanging a bit as to the significance of the differences in u, apparent mass balance, and sfc. elevation when using the different approaches. For example, how do these differences compare to those that arise when using the original BedMachine approach? Would it make sense to include those metrics (differences in u, apparent mass balance, and sfc elevation) somewhere here for comparison? It’s a bit unclear to me what the broader implications are of these secondary metrics w.r.t. using the derived datasets for modeling. If the authors have additional thoughts on this they would be welcome in the supplementary information.

326: If the discussion starting in 4.1 is intended to be specific to Deception, then maybe that should be noted earlier in this paragraph? Alternatively, if the discussion in 4.1 up to line 326 where Deception is mentioned is supposed to be generic, then perhaps line 326 should be something more like, “… a far more realistic bed topography map, particularly for Deception.”

329-330: Would “…slightly higher RMSE SUGGESTS …” be more appropriate here than “indicates”? I think the speculation in this sentence makes sense, but it seems like it is perhaps speculation as opposed to a concrete fact.

336-348: W.r.t. the prediction of thinner ice – is it also possible that this could be the result of the chosen stress balance model? E.g., in order for the SSA model to match surface velocities, it would need to assume a depth-averaged velocity profile that is larger than would be assumed in a model that allowed for internal deformation (E.g., L1L2), because SSA can only accommodate velocity via a change in the sliding component (unless I’m misunderstanding the model used here). If that is indeed the case, then it seems like the optimization process might necessarily bias the ice thickness on the thin side; if the depth averaged velocity is too large, the same flux (constrained by continuity equation and the apparent mass balance terms) can only be accommodated by reducing the ice thickness.

345: “These reasons imply that …” is a bit awkward. “This implies that …”? “These arguments imply that …” ?

353: “… and exceeds their INDIVIDUAL limitations …” ?

398: “We observe that the PINN better captures…” --> “We observe that the PINN captures observable features better with …”

403: “… state variable predictions”. (remove plural on “variable”)

415: “… mass-conserving approach, AS CONFIRMED BY THE DISCOVERY OF new bed features beneath Narssap and Deception.”

417: “… we recommend USING this approach …”

A last thought / general comment: The implied “geomorphology” of the three focus areas studied here look very different from one another. E.g., The Upernavik and Narsaap beds look very smooth when compared to Deception. In the areas where there are no troughs, they almost look like high-resolution DEMs from past, heavily glaciated regions of Canada. Is there any published work on previous Greenland glaciations that might provide some more insight into this? I’m not suggesting it should be part of this paper, but it could be interesting to look into whether or not the “smoothness” that your methods are implying about the bed in different regions is in line with current glacial geological / geomorphological understanding. It would seemingly be a further testament to the power of the methods used here if you were resolving that level of information about the bed through hundreds / thousands of meters of ice.

This paper by Krishna et al. introduces a neural-network-based approach to infer bed topography constrained by both conservation of mass and momentum. The authors find that this approach can infer more physically realistic bed topography in slower-moving regions than BedMachine, which is only mass-conserving. The paper is well written and presents a valuable comparison between the two approaches. As PINN is a new method, such methodological studies are important for building an understanding of their benefits and current limitations. I have a few comments, summarized below:



First, it’s a good idea to hide part of the thickness data for independent validation in Appendix C. As b = s - H, the error reporting in the paper for b is a direct result of the errors in s and H, for which you have training data. My understanding is that the error reported in the main text between the prediction and the “ground truth” topography, along locations where H and s are provided to the NN, simply reflects how well the NN fits the H and s data. It is not relevant to the quality of the physics-based interpolation. A small error along the flight tracks does not imply good b prediction between the tracks, so I’d be cautious about mixing the interpretation of data error in b with the success of the interpolation. That said, in Appendix C the prediction of b at locations without H data is the true demonstration of how effective the physics-informed interpolation is in regions where data are not directly available.

Second, the prediction of velocity from SSA is the depth-averaged velocity, but the observed velocity in the loss function is the surface velocity. Can the authors comment on when this distinction is important and why it is justifiable in the paper? Additionally, is SSA a good approximation of Stokes in all regions studied by the authors?

Third, why does the mass-conserving PINN (MS) produce results so different from BedMachine (Fig. 4), if they both conserve mass? It is mentioned that PINN (MS) predicts isolated, unrealistic crater-like features along radar flight tracks due to “overfitting” (line 323), but why does the mass-conserving BedMachine not suffer from the same issue in the same region? It appears that PINN is using less data than BedMachine, but are they using the same thickness data? Can this difference in input data between PINN (MS) and BedMachine be made more explicit?

Fourth, As PINN solves the equations weakly, the only measure of success, apart from data misfit, is the equation residual. In addition to the various data errors, I believe the paper needs to show map views of the equation residual to demonstrate convincing training success, and the PDE residual should be evaluated on a higher density of collocation points than the training collocation points to check whether the PDE is satisfied between the training collocation points.

Fifth, in the comparison of errors between different inversion results, it is only meaningful to say A is lower than B if the difference between them is larger than the uncertainties in A’s errors. There are many discussions of error comparisons between PINNs and BedMachine, but the PINN errors are averaged over an ensemble of PINN predictions. It is important to consider the spread of errors among PINN predictions. In error reporting such as Table 2, I highly recommend including not only the error of the mean PINN prediction, but also error bars representing the range of error for each PINN prediction. Comparisons between errors are only meaningful after including the uncertainties in PINN errors due to the ensemble.

Finally, would it be possible for the authors to comment on the feasibility of enforcing both momentum and mass balance in the classical adjoint method? If this would be difficult, that would also strengthen the paper’s narrative. This is a natural question that readers are likely to have and will be eager to hear the authors address.

Minor comments:



Line 57-59: Literature review: Bolibar et al., 2023 did not use PINN. It is correct that Riel et al. (2021) is the first PINN study in glaciology. But Riel and Minchew (2023) appeared after some of the other cited studies, and thus I recommend moving it into the sentence along with other papers.

Line 85: “satisfy the PDE residuals” -> “satisfy the PDEs”



Figure 1: Nice figure. Subscript “data” is missing in “L = L_{data} + L_{\phi}”

Line 86-87: Are your collocation points fixed in location or changing throughout iterations? Changing throughout iteration is highly recommended; if not doing so you’ll likely overfit the physics on the fixed discrete collocation points.



Line 115-117: As the apparent mass balance contains both ice thinning rate and the surface/basal mass balance, can you explicitly say how RACMO and ICESat-2 data are combined to give \dot{a}? Does ICESat-2 give a thinning rate dH/dt without the effect of \dot{M}_s,b?



Table 1: Can you use the same time units (either year or second) between the weight values and the variable values for comparisons?

Line 167: Regarding “the PINN output variables will have different values for different regions of the GrIS”, are you using the same weights across the three different regions where the velocities can be very different?



Line 172: I like the fact that different random seeds allow you to sample different solutions that can solve the ill-posed inverse problem. Given the ill-poseness, and the importance of training an ensemble of PINNs, could you elaborate on why 5 time is sufficient, and if you expect different medians if you can train more PINNs?

Eqn 12: Why do you not need data loss for surface elevation s in the loss?



Line 205: Regarding “PINN (SB) is exposed to the ice velocity data along the boundaries of the region of interest, thereby satisfying the stress balance boundary conditions” How many velocity data points do you have along the boundary? I thought the velocity data is 400 points sampled within the domain, meaning that along the boundary of the ROI the velocity data points would be sparse, not truly “satisfying the stress balance boundary conditions”.

Additionally, as you’re also solving for H and s, can you comment on the boundary conditions involving s and H? As they are within spatial derivatives in the momentum equation they also theoretically require BCs just like velocities.



Somewhere in the paper write down \sigma_{SSA} in stress balance in terms of velocities; this will make the discussions of the velocity boundary conditions more clear.

Figure 4 d-f: Does positive denote higher Bedmachine or PINN topography?



Line 259: Regarding “the PINN (MB+SB) RMSE of 137 m is 260 lower than that of BedMachine (see Table 2).” I think you really need error bars for the PINN RMSE to compare which one is lower.

Line 278: Regarding “Lastly, the PINN predicts a ‘disconnected’ trough beneath the southern fork of Upernavik Isstrøm North, while BedMachine suggests that this trough is continuous.” Do you also see disconnected troughs at each PINN result, prior to the averaging across PINN results? When taking the mean between different PINN results do you remove some features that were apparent in each individual PINN results?

Table 4: Why is thickness RMSE not included?