Please find my review as an attached pdf file.

See attached file for review.

Review of Hill et al; "Computationally efficient subglacial drainage modelling using Gaussian Process emulators: GlaDS-GP v1.0"

This paper describes a Gaussian Process emulator of the GlaDS subglacial drainage model and its testing on a synthetic ice-sheet margin setup. Modelling subglacial drainage is starting to become an important aspect of ice dynamics simulations as that system impacts the basal boundary condition significantly. However, subglacial drainage models are relatively costly to evaluate and in particular operate on different, shorter time scales compared to ice flow. Thus running coupled ice-flow drainage simulations is typically difficult and costly at the moment. Emulating the subglacial drainage model using a statistical representation is likely an important step in making these types of coupled models readily applicable.

Whilst emulations of GlaDS with neural network based emulators have been achieved over the last few years, this is the first time a Gaussian Process based emulator has been put forward. The advantage of GP emulators is their greatly reduced number of parameters to fit compared to a neural network as well as built-in capability to quantify uncertainties of the emulation.

The manuscript lays out the procedure to construct the GP emulator; of note is that this construction is relatively involved as it also entails, for instance, decomposition of the GlaDS training data into principal components, fitting of hyperparameters using Bayesian schemes, etc. The emulator is then tested extensively on a synthetic setup and the authors discuss the pros and cons relative to neural network based emulators.

The study and manuscript are carefully constructed. As I am not an expert in statistical emulators, I cannot judge the appropriateness and correctness of the approach to implement the GP emulator. The testing and assessment of the emulator is certainly fine and the discussion is interesting and relevant. Thus, with above caveat, I recommend to publish this manuscript in GMD with the minor corrections outlined below.

Comments

I think it would be useful to discuss a bit more how this emulator could be used for inversions or for coupled ice-flow & drainage simulations as, in my opinion, this are the most sought after usages of such tools.  This can just be in the Discussion and/or Introduction, no need for more simulations or an implementation.

The construction of the emulator has many steps.  Looking through the manuscript, I can see:

- training data construction using parameter design matrix

- running the simulations with GlaDS

- principal component decomposition and component selection or (reduction of variables to scalars)

- fit the GP emulator to the data using an MCMC scheme

Then using the GP in different ways for predictions and analysis is then yet another step. Would it make sense to somehow graphically represent this, flow-chart or some such? Or maybe a numbered list?

Irrespective of the lack of such a graphical overview, I struggled to understand the GP emulator from the description. I am not sure whether I should expect to understand GP emulation from reading about it in such a publication or whether I should just need to go elsewhere to learn it. I see that the authors try to keep the reading smooth by moving quite a bit of the explanations to the appendix but I wonder whether that makes it even harder to understand as now the content is disjoint? Maybe if this layout is kept, then make it even more high level in the main text and have the full description in the appendix which then could be in one place; or, alternatively, move all into the main text? In fact, I think that would be my preferred option and, I think, would fit GMD well as this journal is mostly about methods and not science. As it is, I think it is a bit of a difficult split.

The authors state the principal component decomposition will make the representation necessarily smooth (line 170). Around the channels the hydraulic potential is often not smooth but has the channel as a kink, is that a problem (i.e. a spatial non-smoothness)? Also related to smoothness: in setups like the one presented, where there is no lateral variation in topography, channel position is not necessarily stable with parameter variation but they can jump around (and, for certain, channels move if the mesh is varied). Is that a problem for GP?

Line-by-line

L4: "the combination of the number" is not clearly formulated. Reword.

L8: "daily representation" is not clear to me. Maybe "diurnally averaged"?

L66: I would cite the ISSM GlaDS implementation here too, I think that is Ehrenfeucht&al 2023.

L83: "see B" -> "see Appendix B"

L84: "fast predictions" is a bit sloppy, they are fast to run but not fast themselves.

Table 2: r_b is not defined in the original GlaDS paper nor in this manuscript.  Needs to be defined, at least in Appendix A.

L108: state here that theta is what is fitted and maybe also state the (approximate) size of theta.

L110: "The second choice" really needs a clear statement above of what the first choice is (namely k), otherwise the reader will stumble over this.

L124: $x$ is not defined, or if its definition is "preditction input", then that is not clear enough.

L127: the "posterior distribution" comes out of the blue here

L130: are there d+1 hyperparameters for any k? Couldn't it be less as well? Or more?

L223: A negative floation fraction implies negative water pressure, right? But how can the water pressure go negative in the presented setting? I don't think it can drop below the value of the Dirichlet BC which corresponds to zero water pressure.

L274: RMSE is not defined yet. But then it gets defined in L296.

Fig 2: state to which fields the PCs are encoding

L283: It would be nice to have some snapshots of the PC fields and the GlaDS fields side by side (probably in the appendix). So similar to Fig 2 panel b, but not width averaged but instead just a few instances in time. This would allow to get a bit of a feel on how accurate the spatial fidelity of the PCs are.

Fig 4: Eyeballing the convergence of the two error metrics (panel a,b,d,e), it looks like that the errors do not go to zero but approach some non-zero asymptote. Is that expected? If so, why? Maybe this could be briefly mentioned in the text.

L394: "supporting the interpretation of PC1 as representing water pressure in the absence of surface melt inputs": to me Fig2b1 shows that PC1 has a clear seasonal signal which the basal melt does not. So, I'm not sure this statement is correct or at least needs some more information.

L444: Formulating more clearly what "in ice-flow modelling" means would be helpful

Tab5: here the typesetting seems a bit off: in the fields spanning multiple lines, the line spacing should be less than between different rows.

L552: ideally a DOI and stable archived version of SEPIA and ISSM should also be provided. At the very least the version of ISSM used needs to be stated.

L554: the air-temp dataset needs to be clearly specified. The provided link points to very many datasets.  Do note that this data-repository provides DOIs for each dataset.

Eq A2: I would expect r_b to feature here.

L613: "maximing" -> "maximising"