This paper describes an interesting and potentially very useful methodology for realistic interpolation of sub-glacial topography (one of many potential applications).  Overall I think that this is solid and important contribution to the field.  The methodology is quite complex, however, with many steps involved.  The text is a bit dense with undefined jargon, and I feel the authors could do a much better job at explaining these steps, and particularly in explaining basic concepts.  For example, I was never sure what the authors meant by “distance” between two training images, and that set me at a big disadvantage in trying to comprehend the rest of the methodology.  Another example: the authors never define how the values of MDS1 or MDS2, key parameters in the methodology, are determined.  There are many more such examples noted in my marked-up pdf file.

I hate sounding like the aggrieved reviewer, but really, the authors scant mention of my own paper on the conditional simulation of nearly the same data set had me at a loss.  The two methods are extremely different, but the ultimate product and goals are identical in trying to produce a realistically rough surface conditioned on existing radar soundings and accounting for a high degree of spatial heterogeneity.  Of particular note, my method spent a considerable effort on ensuring the continuity of fjord-like channels beneath the glacier, which are obviously very important factors in flow simulation and likewise are poorly reproduced by standard interpolation schemes like kriging.  How does this method perform in that measure?  I suspect it actually does quite well – that the highest probability deglaciated terrain training images do a good job in conditioning the data interpolation to that geometry.  But the authors do not explore that property.  The authors also did not do a good enough job distinguishing the superiority of their method over SGSIM.  The latter images actually looked quite good.

As noted in my returned pdf, the figures and captions could use some work.  A few of the issues: A lot of the training images were just reduced from larger versions, meaning that the annotation was too small to read.  On several images the color white is used both to indicate areas of no data and Z values >500 m.  This ambiguity needs to be resolved.  Many of the captions were far too brief and failed to explain what is going on in the figure.  

# General comments

The MPS-BedMappingV1 paper focuses on the subglacial topography of the Amundsen Sea Embayment (ASE), which is one of the most important parts of Antarctica to resolve so that ice sheet models can better estimate ice loss and hence sea level rise projections. Using a spatial statistical method, the authors have matched suitable digital elevation models (DEMs) of bed topography from deglaciated areas (in the Arctic and Antarctic off-shore regions) with sparse groundtruth ice-penetrating radar bed elevation data over Pine Island Glacier and Thwaites Glacier. These matched DEMs are training images (TIs) which are used as the basis for a gap-filling algorithm which can be used to fill in the bed topography of areas without data coverage in between radar flightline. The product of this exercise is multiple realizations of subglacial topography over the ASE region, which the authors then use to quantify subglacial hydrological flow uncertainty.

The paper itself was technically challenging, but the authors have done a great job at walking through the procedure step by step, with generous use of figures to illustrate their methodology. One of the core contributions and fascinating pieces of this work is the probabilistic framework of finding training images (from a large pool of deglaciated bed DEMs) that match best with sparse radar data observations over a glaciated area. The amount of work behind sourcing the relevant datasets and building the model framework is no small feat, and I applaud the authors for tacking the time to implement this very important piece of work.

In general, the paper is very readable, though I found that some of the terminology could be made more precise (e.g. high-resolution could be quantified as 500 m spatial resolution). Using subglacial hydrology as the motivation for this subglacial topography modelling work seemed a little weak in my opinion, but I think the authors can do better job to persuade us on what other benefits a high-resolution bed topography model entail. The authors could definitely be a lot more generous with their citations to mention other related interpolation or gap-filling techniques (e.g. Graham et al 2017), and list some more newer and key publications in the cryospheric science field. Following this will be more specific comments and technical suggestions on ways to improve the paper.

# Specific comments and technical corrections

## General point on figure presentation

- In general, the quality of the figures in this paper is great and gets the point across, but being a perfectionist, I have some nitpicks on missing units and coordinates labels. I encourage you to check out Rougier et al., 2014 for some tips on making better figures, both for this paper, and in your future publications.

- Note that the current colour map used in this study may distort the data and is not accessible to people with colour-vision deficiencies. I encourage you to use one of the freely available Scientific colour maps provided at https://www.fabiocrameri.ch/colourmaps or see Crameri et al., 2020 for more details. You can use the package at https://github.com/callumrollo/cmcrameri which integrates with your matplotlib plots.

## Section 0: Abstract

- pg1, L1: Suggest changing "West Antarctica" to "Amundsen Sea Embayment" in the title, since the study doesn't cover some parts of West Antarctica like the Siple Coast or Weddell Sea.

- pg1, L13-14: "However, mapping of subglacial topography is subject to high uncertainty". Please quantify the order of the uncertainty, is it 10 metres? 100 metres?

- pg1, L14-15: "This [high uncertainty] is mainly because the bed topography is measured by airborne ice-penetrating radar along flight lines with large gaps up to tens of kilometers". Note that ice penetrating radar measures bed elevation accurately, I think what you mean to say is that the gaps (without radar coverage) is where the uncertainty lies. Please reword sentence.

- pg1, L19: "We collect 166 high-resolution topographic training images". Please quantify high-resolution, i.e. write high-resolution (500 m).

- pg1, L23: "then provide candidate" -> "then provides candidate".

- pg1, L24: "traditional MPS" -> "traditional multi-point geostatistics (MPS)".

- pg1, L25: "demonstrates significant improvement". Please quantify improvement over previous baseline, e.g. "significant improvement of 10 metres"

- pg1, L27: "We use the multiple realizations to investigate the impact of basal topography uncertainty on subglacial hydrological flow patterns". This sentence should be expanded to mention what the results of the investigation are exactly. From reading Section 4.2, besides seeing that multiple subglacial hydrological realizations were made, I could not quantify how the method here improves upon an existing water flow model (e.g. from the subglacial water flux map of Le Brocq et al., 2013).

## Section 1: Introduction

- pg1, L29: "The topography beneath the Greenland and Antarctic ice sheets is essential for nearly every ice sheet investigation". Need to cite more examples, the following sentences (L30-L32) only include citations for Antarctic papers, and are not representative of the cryospheric literature. I have listed a few suggested citations in the next bullet points, but you are welcome to add others too.

- pg1, L30: "modeling subglacial hydrology (MacKie et al., 2021)". Please cite more, e.g. Siegert et al., 2016; De Fleurian et al., 2018; Sommers et al., 2018; etc

- pg1, L30: "interpreting geologic conditions (Holschuh et al., 2020)". This sentence seems incorrect, how can bed topography be used to interpret geologic conditions? The cited article by Holschuh et al., 2020 mentions "a likely role for preexisting geology in glacial bedform shape", i.e. that geology affects topography, not the other way round. Please reword the sentence, and add more Greenland/Antarctic subglacial geology papers e.g. Anandakrishnan et al., 1998; Aitken et al., 2014; Bell et al., 1998; Lowry et al., 2020; etc

- pg1, L31: "estimating ice volume and sea level rise contributions (Fretwell et al., 2013)". Please also cite BedMachine papers by Morlighem et al., 2017; Morlighem et al., 2019.

- pg1, L31: "ice sheet modeling for sea level rise projections (Le clec’het al., 2019; Seroussi et al., 2017)". Please cite more, e.g. Schlegel et al., 2018.

- pg2, L35: "radar along flight lines separated by up to tens of km (Fretwell et al., 2013; Herzfeld et al., 1993)". Please cite more of the underlying radar data sources, e.g. Robin et al., 1970; Shi et al., 2010; etc. Also, you should mention ice penetrating radar surveys with closer sub-kilometre spacing e.g. over iSTAR data over Pine Island Glacier by Bingham et al., 2017 and Rutford Ice Stream data by King et al., 2016.

- pg2, L36: "using methods such as kriging (Fretwell et al., 2013)". BEDMAP2 was not created using kriging, they used the ArcGIS Topogrid algorithm (see pg 381 of Fretwell et al., 2013). Please find an alternate paper citation that uses kriging, or change kriging to Topogrid.

- pg2, L37: "model inversions (Huss and Farinotti, 2012; Morlighem et al., 2017, 2020)". Recommend citing ITMIX paper by Farinotti et al., 2017 too. Also note that ice sheet model inversion techniques are not deterministic, I suggest you separate this sentence from the previous L36 "generally interpolated deterministically using methods such as kriging" part.

- pg2, L37: You may also want to mention other spatial statistical papers by Graham et al., 2017 and Goff et al., 2014 somewhere here for completeness.

- pg2, L49-50: "contain a globally wide range" -> "contain a range". Make this sentence more concise.

- pg2, L59: "non-stationary bed topography is directly measured using high-resolution remote sensing data such as satellite images". Bed topography is not measured directly by optical satellite images (e.g. Landsat/Sentinel-2), they can be measure indirectly using photogrammetry. Direct measurements would require using satellite altimeters such as Cryosat-2 or ICESat-2. The ArcticDEM dataset by Porter et al., 2018 you cited uses stereo-photogrammetry on WorldView imagery. Please reword this sentence, and also mention which version of ArcticDEM was used in this study (is it v3.0?).

- pg2, L60-61: "reveal glaciated morphologies resembling the topography beneath the ice sheets" -> "reveal glaciated morphologies resembling the topography beneath the contemporary ice sheets". Need to be clear that you are using deglactiated paleo-ice-sheets that resemble the current present-day ice sheet.

- pg2, L61-62: "They therefore bear significant information on the subglacial topography". Clarify exactly what is meant by 'information', something to do with the texture or structure of the topography? Or something else?

- pg2, L63: "satellite imagery of deglaciated topography has not been explored to stochastically simulate subglacial topography". Again, optical satellite imagery does not contain direct topography information without extra processing. Be clear that you mean satellite imagery derived DEMs.

- pg3, L66: "high-resolution training images" -> "high-resolution (500 m) training images".

- pg3, L76-77: "Then MPS then uses" -> "MPS then uses".

- pg3, L83: "More importantly," -> "Moreover,". Be concise.

- pg3, L88-89: "there are no available training images because" -> "there is limited training images because".

- pg3, L89-90: "Satellite altimetry observations from deglaciated areas in the Arctic offer a potential source of training imagery". Note that satellite altimetry products (e.g. from Cryosat-2/ICESat-2) are typically point-based measurements, not raster grid images. Suggest rewording this sentence to "Satellite-based observations of deglaciated areas...".

- pg3, L90: "training imagery. However, training images retrieved from the Arctic would be" -> "training imagery, but these training images would be". Be concise.

- pg3, L92: "significant computational burden". Could you quantify what significant means?

- pg3, L96: "We will address" -> "We address". Be concise.

- pg3, L98-99: "We first collect a large amount of topographic images to serve as the training images. These images are taken from the deglaciated areas in the Arctic and Antarctica" -> "We first collect topographic images from deglaciated areas in the Arctic and Antarctic regions to serve as training images (Fig 2)". Be concise and link to Fig 2.

- pg4, L103: "We will demonstrate our method using the entire" -> "We demonstrate our method over the entire".

- pg4, L106-107: "We use the topographic simulations to" -> "We then use these simulated topographic images to"

- pg4, L107: "in order to investigate the impact of topographic uncertainty on hydrologic" -> "in order to investigate how topographic uncertainty impacts hydrologic"

## Section 2: Radar data set & training images

- pg4, L110-118: Recommend to provide a table of all the data sources used since the list is long. At a minimum, the table columns should include the placename of the data (e.g. Thwaites Glacier, Laurentide Ice Sheet, etc), the data citation, and the paper citation.

- pg4, L110: "seafloor bathymetry measurements". Fig 1 map doesn't show bathymetry, i.e. areas outside the grounding zone.

- pg4, L111: "subaerial topography" -> "ice surface topography". Also, could you show where the REMA DEM is being used over the ASE region? There are not many deglaciated places in the ASE on land.

- pg4, L113-114: "gridded at a 500-meter resolution" -> "gridded at a 500-meter spatial resolution". Specify the gridding algorithm used, and how data gaps (i.e. NaN values) are treated.

- pg4, L114-115: "(Holschuh et al., 2020) (provide" -> "(Holschuh et al., 2020) provide"

- pg4, L115-116: "we augment the available training data with deglaciated subaerial topography from the ArcticDEM (Porter et al., 2018)" -> "we increase the available training data with deglaciated subaerial topography from ArcticDEM data (Porter et al., 2018)". Augmentation can have a slightly different connotation in machine learning, such as synthetic data generation, suggest to just use 'increase'.

- pg4, L119: "may have experienced additional depositional processes" -> "may have experienced additional erosional or depositional processes".

- pg4, L119-120: "any topographic alterations are likely minimal at a 500 m resolution". How is it likely minimal? Need to provide a citation or have some numbers to back up this claim.

- pg4, L120: "geologic settings" -> "geological settings".

- pg4, L121: "100x100 km^2" -> "100 km x 100 km".

- pg4, L121-122: "training image data repository is publicly accessible from Zenodo repository (https://zenodo.org/record/5083715#.YQT2JI5Kiiw, DOI 10.5281/zenodo.5083715)" -> "training image data is publicly accessible on the Zenodo repository (https://doi.org/10.5281/zenodo.5083715)".

- pg5, L123: Figure 1 colorbar. Suggest using Scientific Colour Map instead of "terrain" because values at 1500 m are white, which is the same as the NaN values. You may want to select a different colour for NaN values and/or indicate with a small box that NaN values are of a certain colour.

- pg5, L123: Figure 1 axes. Please show x/y (or longitude/latitude) coordinates for this map with units.

- pg5, L127: Need to indicate that the inset map's red box shows the study region.

- pg5, L128: "The topography patches" -> "The two topography patches".

- pg6, L131: Figure 2a axes. Please show x/y (or longitude/latitude) coordinates for this map with units.

- pg6, L131: Figure 2a map. What do the red boxes with the checkerboard patterns mean? Mention in the caption.

- pg6, L131: Figure 2b axes. Please show x/y (or longitude/latitude) coordinates for this map with units.

- pg6, L131: Figure 2b map. Why are there some black lines on the main map? Black does not seem to be in the colorbar.

- pg6, L132: Figure 2c axes. Please show x/y (or longitude/latitude) coordinates for this map with units.

- pg6, L132: Figure 2c map. T1-165 and T1-166 look like DEMs from Holschuh et al., 2020, are they rotated 90 degrees? Where is North? Ideally there would be x/y or lon/lat coordinates to show the location and orientation of these training images.

- pg6, L132: Figure 2 caption. "(b) Antarctica". Are training images gathered from the whole of Antarctica, or just the Amundsen Sea Embayment region as shown in Fig 2b?

## Section 3: methodology

### Section 3.1: Multiple-point geostatistics

#### Section 3.1.1: Overview

- pg7, L142-143: "Both MPS and covariance-based methods have the ability to interpolate data exactly". Please clarify what "exactly" means, e.g. is the vertical precision 1 metre? 5 metres?

- pg7, L149: "section 0" -> "section 3.1.2".

#### Section 3.1.2: Direct Sampling (DS)

- pg7, L164: "a unknow" -> "an unknown".

- pg7, L165: "to find the similar structure in TI. The similarity within" -> "to find a similar structure in the TI pool, Similarity within".

- pg8, L166: "defined by a certain distance metric". What are examples of distance metrics that can be used? Suggest change sentence to "defined by a distance metric (e.g. Euclidean distance, Hausdorff distance, etc)".

- pg8, L170: Figure 3b. Put a dashed box (width: 4 pixels, height: 3 pixels) around the [37], [?][86][80] cells, otherwise they seem disconnected. It took me a while to understand what was happening.

- pg8, L170: Figure 3 colorbar. Add units for the colorbar.

- pg8, L174: "Based on the explanation above, there are mainly three important parameters in DS" -> "Based on the explanation above, there are three main parameters in DS to tweak".

- pg8, L178: "The TI pattern with mismatch distance below t will be accepted and pasted to the simulation grid". What happens when 2 or more TI patterns are found below the threshold? Is the one with the lower distance metric value (i.e. more similar) picked? Or is it done on a first comes first serve basis? Please add a sentence to clarify in text.

- pg8, L182: "avoid verbatim copy, an recommended value of" -> "avoiding verbatim copies, a recommended value of"

#### Section 3.1.3: A metric space for training images

- pg9, L192-193: "we rescale each TI to range between 0 and 1 by min-max normalization (Han et al., 2012)". Are the min and max values used for normalization obtained per image (i.e. local) or derived from the entire training dataset (i.e. global). Depending on the implementation, I have some concerns on how suitable this normalization is. Consider a training image with high mountains and low valleys, and a training image which is relatively flat and lacks topographic features. A local min-max normalization would stretch both images into the same 0-1 range, i.e. little bumps on the flat terrain training image would become mountains. In effect, this could lead the subsequent model to use training images from flat areas to inform the gap filling of mountainous areas.

- pg9, L194: "from each TI with a" -> "from each TI using a".

- pg9, L195-196: "into a finite number of groups". How many groups? Please state the number here.

- pg9, L197-198: "We referred to the commonly used distance threshold in DS approach to set it as 0.1 (Meerschman et al., 2013) of the maximum pattern distances of the TI". Sentence not clear, maybe reword to "We set the distance threshold for DS to 0.1, as per Meerschman et al., 2013"? Also, the Meerschman et al., 2013 paper uses 0.05 as a default (though they have a range of recommended values), could you elaborate why 0.1 used here instead, and not another value like 0.2?

- pg9, L200: "The distance used in the clustering is the normalized Euclidean distance". The agglomerative hierarchical clustering scheme seems to typically works on points, but your training images are raster grids (with a z value), so how is the Euclidean distance computed? Are you treating each pixel as a point with x,y,z values and running the clustering on those, or are you computing the pixel-wise distance (i.e. elevation difference) between the grids and taking the mean? Please clarify.

- pg9, L202: Figure 4 colorbar. Add units for colorbar (e.g. Normalized Elevation). Also, does the min-max normalization turn NaN values to 0 (because I what looks like data gaps as dark blue, which is 0 on the colorbar), or are NaN values still treated as NaN?

- pg9, L204: "training images are now represented expressed" -> "training images can now be expressed".

- pg9, L207: "used to define the" -> "used to quantify the".

- pg9: L211: "x_A is any representative pattern in A". The term "any representative pattern" can sound a bit ambiguous, suggest rewording to clarify that "x_A is any one of the possible representative patterns of A".

- pg9, L211-212: "d(. ) is the Euclidean distance between any two representative patterns". Again, see my point above on pg9, L200. Is the Euclidean distance computed by vectorizing the pixels to points first, or pixel-wise on the images?

- pg9, L212-213: "the modified Hausdorff distance". My understanding of Hausdorff distance is that it is a method originally developed for x, y point cloud measurements, but your training images are raster grids. Depending on whether you are treating each pixel in your training images as x, y, z points, or computing the Euclidean distance pixel-wise (i.e. distance is calculated for pixel-values in the same image coordinates), this can result in very different results, so you will need to clarify and possibly justify your approach.

- pg10, L219: Figure 5 colorbar. Add units for the colorbar.

- pg10, L220: Figure 5 caption. "TI" -> "training image (TI)".

### Section 3.3: Formulation of the problem through probability aggregation

- pg11, L245. Equation 3. Need to state what is i and j.

- pg12, L258: "additional weight term" -> "additional weight term (w_{ij}) as follows"

### Section 3.4: Probability of training images given radar line data.

#### Section 3.4.1 Most probable set of training images

- pg12, L268: "challenging because the d_{Ai} are very high-dimensional" -> "challenging because the flight radar line data d_{Ai} is high-dimensional".

- pg12, L269: "low-dimensional" -> "low-dimensional problem"

- pg12, L277: Equation 8. Does this function yield n? Maybe do "n = argmin..."

- pg12, L279: "size n of TI in the Appendix" -> "size n of TI via Particle Swarm Optimization (PSO) in the Appendix".

- pg13, L283: Equation 9. This does not look like a Hausdorff distance formula like in Equation 1, more like a directed Hausdorff with no max argument. Please check that this formula is correct.

- pg13, L285: "We use flexible sized templates". Please state the range (min and max) of template sizes used, e.g. 10x10, 20x20, etc.

- pg13, L286: "randomly chosen between the maximum radius up to 15 pixels to include 40 measurement points". This sentence carries ambiguity. Firstly, radius implies a circle, is that what was used, or was a rectangle used as in Figure 3? Secondly, does 'include 40 measurement points' mean that only 40 measurements points are needed and no more (i.e. =40), or that it is a minimum (i.e. >=40).

- pg13, L288: Figure 7 colorbar. Add units for the colorbar.

- pg13, L289: Figure 7 caption. "training image and radar line data" -> "training image (TI) and radar line data (d)".

- pg13, L291: "We use a particle swarm optimization (PSO) to minimize the distance function dis(.)" -> "We use particle swarm optimization (PSO) to minimize the distance dis(.) (see Appendix)".

- pg13, L292: "PSO has its specific advantages in requiring less parameterizations". Are the parameters you used in the particle swarm optimizer for this Amundsen Sea Embayment study applicable to other parts of Antarctica, or would people need to re-run things to find new optimal parameters? I.e. how stable are the parameters across different study areas?

- pg14, L298: Figure 8b colobar. Add units for the colorbar.

- pg14, L298: Figure 8b axes. Please show x/y (or longitude/latitude) coordinates for this map with units.

- pg14, L298: Figure 8b map. Is A1=T1-99, A2=TI-88 and so on? Not sure why the TI numbers are on the bottom of the map.

- pg14, L299. Figure 8 caption. "for each area displayed on MDS plots" -> "for each subarea (A1, A2, A3, A4) displayed on multi-dimensional scaling (MDS) plots".

- pg14, L299-300. Figure 8 caption. "The red dots highlight the estimated TI" -> "The blue dots represent a TI, and the red dots highlight the estimated TI". Also, why do some subareas have two red dots while some has three?

#### Section 3.4.2: Kernel density estimation of P(TI(A_i )|d_{A_i})

- pg15, L303: "We assume that the TIs" -> "We assume that TIs".

- pg15, L304: "(Figure 4)" -> "(Figure 5)".

- pg15, L305: "we can observe the spatial patterns of nearby TIs look similar" -> "we can observe that spatial patterns of nearby TIs look similar". Please define how it looks 'similar'. Is it a qualitative similarity, i.e. the topographic structure looks similar, or something else?

- pg15, L306: "predict the probability to each TI". Probability of what to each TI? Please state.

- pg15, L311: "dis(TI, TI_k) is distance between a TI and" -> "dis(TI, TI_k) is the modified Hausdorff distance between a TI and".

- pg15, L314: Figure 9 colorbar. Add units for the colorbar.

- pg15, L316: "(a) Estimated TI for A1" -> "(a) Estimated TI for subarea A1".

- pg15, L316: "plotted in MDS space" -> "plotted in multi-dimensional scaling (MDS) space".

- pg15, L316: "The red dots are" -> "The blue dots are training images (TI) and red dots are".

### Section 3.5 Aggregation by weighting log-ratios

- pg16, L323: "two areas" -> "two subareas".

- pg16, L327: "between area" -> "between subarea"

- pg16, L328: "x_l and x_l' are the radar data locations". Clarify what this means. Are they point (x,y) or pixel (image) coordinates?

- pg16, L328: "Epanechnikov kernel". Need to provide a citation, or better, state the formula.

- pg16, L329-330: "the weights w_{ij} are simply" -> "the weights w_{ij} are"

- pg16, L332: Equation 13. "where i, j = [1,2,3,4]". Does this mean i or j can be 1,2,3 or 4? I.e. 4x4 = 16 possibilities? I feel like there is a better way to state this mathematically. Alternatively, you could take this out of Equation 13 and describe i and j in-text at L334.

- pg17, L338: Figure 10. (Optional) Would be nice to plot the A1-A4 images on the top right corner of each panel, or at least refer to the subarea images in the caption, i.e. link to Fig 8b.

- pg17, L338: Figure 10 axes. Please spell out PDF in full.

- pg17, L339: Figure 10 caption. "aggregated TI" -> "aggregated training image (TI)"

### Section 3.6 Direct sampling with TI sampling

- pg17, L344: "tend to have higher elevations". How many metres higher? 2000m? 3000m?

- pg17, L344: "larger scale low-elevation valleys". How many metres deep, and how many metres wide?

- pg17, L345: "that can be related to warm water routing". This is contentious, the valleys could just be due to ice stream action.

- pg17, L346: "At the end, multiple realizations of topographical models are generated with" -> "At the end of the simulation, multiple realizations of topography are generated, each with"

- pg18, L351: Figure 11 caption. "Examples of sampled TIs from the posterior distributions" -> "Examples of sampled training images (TIs) from the posterior distribution corresponding to each subarea (A1, A2, A3, A4) for each realization". You might need to rephrase this sentence a bit to explain the figure better.

- pg19, L353. Figure 12 axes. Please show x/y (or longitude/latitude) coordinates for this map with units.

- pg19, L353: Figure 12 colorbar. Add units for the colorbar.

- pg19, L355: Figure 12 caption. "DS" -> "direct sampling (DS)"

- pg19, L355-356: Figure 12 caption. "Model realization number corresponds to the TI realization number in Figure 11". I think you mean "Model realization number (#n) corresponds to the model realization number in Figure 11"?

### Section 3.7 Comparison with traditional MPS modeling and two-point geostatistical modeling

- pg19, L360-372: This subsection seems to be interpreting some results, and you may want to consider moving it into a 'Discussion' section for the final publication.

- pg19, L361: "scanning all the 166 TIs" -> "scanning all 166 TIs".

- pg19, L361-362: "Figure 13a shows one realization of the simulated result. It is obvious that the conventional approach results in a much noisier topographical model. There are significant" -> "Figure 13a shows one realization of the conventional MPS simulated result. The conventional approach appears to show a much noisier topographical model, and there are significant"

- pg19, L362-363: "we take a cross-section A-A’ on the Pine Island glacier and plot the comparison" -> "we take a cross-section A-A' across the trunk of Pine Island Glacier and plot the elveation comparison".

- pg19, L364-365: "We can observe that the DS without TI sampling creates a significant amount unrealistic" -> "We observe that the DS without TI sampling methods generated more unrealistic"

- pg19, L365: "Pine Island glacial" -> "Pine Island Glacier".

- pg19, L366: "where mostly radar data" -> "where dense radar data".

- pg19, L367-368: "when using all the 166 TIs without proper sampling, the DS finds too large a set of patterns likely many incompatible with the sparse data" -> "when using all 166 TIs without proper sampling, DS without TI sampling finds too large a set of patterns that is likely incompatible with sparse data".

- pg19, L369: "thereby improving the result" -> "thereby improving the result (see Fig 12)".

- pg19, L369-370: "avoiding the channel artifacts is critically important for modeling subglacial hydrological flow (see section 0)" -> "avoiding topographic artifacts in subglacial channels is important for modeling subglacial hydrological flow (see section 4.2)".

- pg19, L371: "to simulate one realization. When using our TI sampling approach, it took less than 1 hour" -> "to simulate one realization, compared to 1 hour for our TI sampling approach".

- pg19, L372: "are run on a PC with an Intel i9-11900 of 2.5GHz processor and 32GB of RAM" -> "were ran using an Intel i9-11900 (2.5GHz) processor with 32GB of RAM".

- pg20, L374: "compare to the two-point geostatistical modeling with kriging and Sequential Gaussian Simulation (SGSIM)" -> "compare the two-point geostatistical modeling results (Fig 12) with kriging (Fig 13b) and Sequential Gaussian Simulation (SGSIM; Fig 13c)".

- pg20, L375: "Figure 13b and Figure 13c plot topographic modeling results from kriging and SGSIM" -> "".

- pg20, L375-376: "We can observe that kriging produces the most smoothed topographical model" -> "We observe that kriging produces the smoothest topography".

- pg20, L376: "are very clear" -> "are clear".

- pg20, L377-378: "After all, kriging is a deterministic modeling approach. Thus, it cannot capture location scale elevation variations and quantify the spatial uncertainty" -> "After all, kriging is a deterministic modeling approach which cannot capture variations in local scale elevation, nor can it quantify the spatial uncertainty".

- pg20, L378: "quantify the spatial uncertainty". Note that kriging can output variance, which could be a measure of uncertainty.

- pg20, L378: "Our SGSIM" -> "The SGSIM".

- pg20, L379: "limited the neighborhood" -> "limiting the model to the neighborhood".

- pg20, L379-380: "The limitation of SGSIM, an approach based on spatial covariances, lies on its limitations in capturing" -> "However, SGSIM, an approach based on spatial covariances, is limited in its ability to capture".

- pg20, L382: "using the four different approaches". Technically the radar is the groundtruth, and you are comparing three different approaches. Please reword the sentence.

- pg20, L382: "the DS using sampled TIs has" -> "the DS using sampled TIs method has"

- pg20, L384: "the DS without TI sampling has" -> "the DS without TI sampling method has"

- pg21, L386: Figure 13 map. It would be nice if you can combine Fig 12 and Fig 13 so that the different methods (MPS, Kriging, SGSIM) could be directly compared in one plot. Even better if you can plot the groundtruth data lines (e.g. from Fig 15a) too.

- pg21, L386. Figure 13 axes. Please show x/y (or longitude/latitude) coordinates for this map with units.

- pg21, L386: Figure 13 colorbar. Add units for the colorbar.

- pg21, L387: Figure 13 caption. Please mention that the black ovals on Fig 13c refers to areas of sparse data.

- pg22, L388: Figure 14 top map. Please include a colorbar, and refrain from using jet as a colormap, choose a Scientific Colour Map instead.

- pg22, L388: Figure 14 top map. Can you put the dashed box from Fig 14 bottom (the transect plot) in the top map too?

- pg21, L388. Figure 14 top axes. Please show x/y (or longitude/latitude) coordinates for this map with units.

- pg22, L388: Figure 14 bottom plot. Don't use red and green for the lines as they are not good for people with Deuteranopia. Also, suggest changing the black colour for the kriging line to a different colour to avoid confusion with the radar data points. Try https://colorbrewer2.org for suggestions on what to use.

- pg22, L389: Figure 14 caption. Mention that the dashed box shows main channel of Pine Island Glacier.

- pg22, L390: Figure 15. Don't use red and green for the lines.

- pg22, L390: Figure 15. For the SGSIM line (blue), which realization are you using? 1 or 2?

- pg22, L390: Figure 15. Could you describe in text, why the lines deviate away from the radar data more at 70-80km? Is it because of sparser radar data coverage?

- pg22, L390: Figure 15 axes. Why does the x-axis stop at 80km, and not continue to 100km like in Fig 14?

## Section 4 Application to the entire Amundsen Sea Embayment (ASE)

### Section 4.1 Training image sampling and DS simulation

- pg23, L395: "We apply the methodology to" -> "We apply the methodology in Section 3 to".

- pg23, L395: "high-resolution topography" -> "high-resolution (500 m) topography".

- pg23, L397: "diving" -> "dividing".

- pg23, L398: "Equally divide the ASE area into four subareas". By four subareas, do you mean quadrants? Or 4 long strips? Are the subareas always square, or can the be rectangular?

- pg23, L399-401: The loop in Step 2 and 3 says "continue to divide it into four equally sized areas... until the threshold", does this mean that subarea sizes are not fixed at 100km x 100km? Does the model work with different sized subareas?

- pg23, L400: "Repeat step 2 until amount of data" -> "Repeat step 2 until the amount of data".

- pg23, L401: "into totally of 56 subareas" -> "into a total of L=56 subareas".

- pg23, L403: "Figure 16 shows the final ASE subareas with the corresponding radar data density". Please specify how radar data density is calculated. Is it the number of non-NaN pixels over an area?

- pg23, L409: Figure 16 map. Need a colorbar for the radar elevation data, and I highly encourage using a Scientific Colour Map instead of jet.

- pg23, L409. Figure 16 axes. Please show x/y (or longitude/latitude) coordinates for this map with units, instead of using image coordinates.

- pg24, L411. Figure 17 axes. Please show x/y (or longitude/latitude) coordinates for this map with units, instead of using image coordinates.

- pg24, L411: Figure 17 colorbar. Add units for the colorbar.

- pg24, L411: Figure 17 caption. "TIs assigned to the entire ASE area" -> "training images (TIs) assigned to the entire Amundsen Sea Embayment (ASE) area"

- pg25, L413. Figure 17 axes. Please show x/y (or longitude/latitude) coordinates for this map with units.

- pg25, L413. Figure 17 colorbar. Add units for the colorbar. Also, find a way to represent NaN values properly instead of as white colour.

- pg25, L411: Figure 17 caption. "ASE" -> "Amundsen Sea Embayment (ASE)".

- pg25, L413: Figure 17 caption. "circles" -> "ovals".

## Section 4.2 Uncertainty in subglacial hydrological flow

- pg26, L419: "model was applied to 20 realizations to model" -> "model was applied to the 20 topography realizations in Section 4.1 to model".

- pg26, L419:420: "The direction of water flow is determined by calculating hydraulic potential, phi, using the Shreve equation (Shreve, 1972)" -> "The direction of subglacial water flow can be determined by the gradient of hydraulic potential, phi, which is calculated using the following equation (Shreve, 1972)".

- pg26, L424: "100 kg m-3" -> "1000 kg m-3"!!

- pg26, L424: "gravitational acceleration" -> "gravitational acceleration (9.8 m s-1)". Or whatever value you used.

- pg26, L425: "The hydrological model is implemented" -> "The subglacial hydrological model was implemented".

- pg26, L425-426: "These functions use the hydraulic potential gradient to compute flow accumulation, or the number of pixels that flow into another pixel". Describe the algorithm briefly for someone unfamiliar to the tools, e.g. is it using a D8 routing model?

- pg26, L428: "We assume spatially uniform basal melt rates". Specify the basal melt rate quantity used.

- pg26, L428: "water pressure" -> "subglacial water pressure".

- pg26, L432: "grounding line" -> "grounding zone".

- pg26, L433: "These tributaries are located near a system" -> "Some of these tributaries are located over a system".

- pg26, L434: "(Smith et al., 2017)". Please also cite more recent Thwaites Glacier active subglacial lake papers by Hoffman et al., 2020 and Malczyk et al., 2020.

- pg26, L434-435: "with a drainage and refill cycle" -> "with a drain and refill cycle".

- pg26, L435-436: "Lake drainage events are associated with increases in ice velocity (Stearns et al., 2008), making it important to characterize the connectivity of active lake systems". This is stated too simplistically. Lake drainage doesn't always cause ice speedup, see Section 4.2 of Smith et al., 2017. You can however, say that "Lake drainage events are sometimes associated with increases in ice velocity", though I recommend mentioning the caveat that lake drainage may not necessarily lead to ice speedup.

- pg26, L436: "(Stearns et al., 2008)". Please consider reading and citing other relevant papers such as Wright et al., 2014 and Fricker et al., 2016.

- pg26, L436: "Our results" -> "Our topography modelling results in Section 4.1".

- pg26, L436-437: "additional observational constraints" -> "additional radar data observational constraints".

- pg26, L436-438: You will need to better quantify and justify why a better bed topography is needed, since ice surface elevation is 11x more important than bed elevation for hydropotential calculations (which can be derived from Equation 14). Otherwise it will be hard to argue your case that additional bed elevation observations are needed for deriving better hydropotential maps when ice surface elevation is the more important driving force.

- pg27, L439: Figure 19 axes. Please show x/y (or longitude/latitude) coordinates for these maps with units.

- pg27, L439: Figure 19 maps. Mention in the caption where you obtained the glacier boundaries, provide a data citation if possible.

- pg27, L439: Figure 19a map. Is the Mass conservation map derived using BedMachine Antarctica data? If so, please mention in the caption and cite accordingly.

- pg27, L439: Figure 19e map. You have shown the mean flow accumulation which is good, but maybe show a standard deviation plot as well to see how the different realizations differ, and where they differ the most. Also, it would be nice to show a difference map between Fig19d and Fig19e to highlight the uncertainty of the subglacial water flow model between the mass conservation map and your different realizations.

## Section 5 Conclusions

- pg27, L446-447: "fill large-scale geophysical data gaps and applied it to map high-resolution subglacial topography in West Antarctica" -> "fill geophysical data gaps and applied it to produce a high-resolution (500 m) subglacial topography map of the Amundsen Sea Embayment".

- pg28, L448: "high-resolution topographic" -> "high-resolution (500 m) topographic".

- pg28, L449-450: "from the Arctic and Antarctica. These training images represent" -> "from deglaciated regions in the Arctic and Antarctica, which represent".

- pg28, L450: "We have placed them in" -> "We have placed the training images (TIs) in".

- pg28, L450: "publicly available repository". Link to the repository here, or point to the Data and code availability section.

- pg28, L457: "based on the distance" -> "based on the modified Hausdorff distance".

- pg28, L459: "posterior TI distribution allowed us" -> "posterior TI distribution then enabled us".

- pg28, L459: "direct sampling" -> "direct sampling (DS)".

- pg28, L461: "Such non-stationary TI sampling framework avoided the use" -> "Such a non-stationary TI sampling framework avoids the use".

- pg28, L461-462: "It significantly improved" -> "It has significantly improved".

- pg28, L462: "reduced the DS running time" -> "reduced the DS running time from 21 hours to 1 hour".

- pg28, L470: "across realizations" -> "across different realizations".

- pg28, L470: "tributaries are near a system" -> "tributaries cross a system".

- pg28, L471-472: "and could have the potential to influence ice sheet velocity. The high hydrological uncertainty in this area highlights the need for additional measurement constraints". Similar to my point made above for pg26, L436-438, you will need to argue your case better to arrive at this conclusion, because ice surface elevation is 11x more important than bed elevation for calculating hydropotential, so additional measurement constraints (assuming you mean radar bed data) will not affect hydropotential as much as better ice surface elevation data from satellite altimeters. My suggestion is that you can also mention why high-resolution (<500 m) fine-scale topography matters, see e.g. Kyrke-Smith et al., 2018 and related papers.

## Section Appendix

- pg30, L520: Figure 20a axes. The y-axis needs units for distance, is it metres? Also, the x-axis says "size n of TI", but I think you mean "swarm size n" according to the caption.

## Section Code and Data availability

- pg30, L524: "training image database". The training image (TI) database of 166 images was stored as a single text file. In order to be more user friendly (see FAIR guidelines by Wilkinson et al., 2016), I strongly advice that you 1) separate the TIs to be individual files, one for each image, so that people don't need to figure out how to reshape the arrays; 2) use an alternative file format like NetCDF or GeoTIFF, so that you can have geographical coordinates and other metadata associated with each TI, especially important for the sake of data provenance (where was the TI sourced from? The Arctic or Antarctic?). Also, you have not included any of your output topography or subglacial flow model realizations on Zenodo, but I expect those to be released in time for the final publication.

- pg30, L525: "modelling codes" -> "modelling code".

- pg30, L526: "https://github.com/sdyinzhen/MPS-BedMappingV1". I highly appreciate the source code being made available. For the sake of reproducibility, could you provide a dependency specification (e.g. a pip requirements.txt or conda environment.yml) file so that others know what libraries you have used?

## References

Aitken, A. R. A., Young, D. A., Ferraccioli, F., Betts, P. G., Greenbaum, J. S., Richter, T. G., Roberts, J. L., Blankenship, D. D., & Siegert, M. J. (2014). The subglacial geology of Wilkes Land, East Antarctica. Geophysical Research Letters, 41(7), 2390–2400. https://doi.org/10.1002/2014GL059405

Anandakrishnan, S., Blankenship, D. D., Alley, R. B., & Stoffa, P. L. (1998). Influence of subglacial geology on the position of a West Antarctic ice stream from seismic observations. Nature, 394(6688), 62–65. https://doi.org/10.1038/27889

Bell, R. E., Blankenship, D. D., Finn, C. A., Morse, D. L., Scambos, T. A., Brozena, J. M., & Hodge, S. M. (1998). Influence of subglacial geology on the onset of a West Antarctic ice stream from aerogeophysical observations. Nature, 394(6688), 58–62. https://doi.org/10.1038/27883

Bingham, R. G., Vaughan, D. G., King, E. C., Davies, D., Cornford, S. L., Smith, A. M., Arthern, R. J., Brisbourne, A. M., De Rydt, J., Graham, A. G. C., Spagnolo, M., Marsh, O. J., & Shean, D. E. (2017). Diverse landscapes beneath Pine Island Glacier influence ice flow. Nature Communications, 8(1), 1618. https://doi.org/10.1038/s41467-017-01597-y

Crameri, F., Shephard, G. E., & Heron, P. J. (2020). The misuse of colour in science communication. Nature Communications, 11(1), 5444. https://doi.org/10.1038/s41467-020-19160-7

De Fleurian, B., Werder, M. A., Beyer, S., Brinkerhoff, D. J., Delaney, I., Dow, C. F., Downs, J., Gagliardini, O., Hoffman, M. J., Hooke, R. L., Seguinot, J., & Sommers, A. N. (2018). SHMIP The subglacial hydrology model intercomparison Project. Journal of Glaciology, 64(248), 897–916. https://doi.org/10.1017/jog.2018.78

Farinotti, D., Brinkerhoff, D. J., Clarke, G. K. C., Fürst, J. J., Frey, H., Gantayat, P., Gillet-Chaulet, F., Girard, C., Huss, M., Leclercq, P. W., Linsbauer, A., Machguth, H., Martin, C., Maussion, F., Morlighem, M., Mosbeux, C., Pandit, A., Portmann, A., Rabatel, A., … Andreassen, L. M. (2017). How accurate are estimates of glacier ice thickness? Results from ITMIX, the Ice Thickness Models Intercomparison eXperiment. The Cryosphere, 11(2), 949–970. https://doi.org/10.5194/tc-11-949-2017

Fricker, H. A., Siegfried, M. R., Carter, S. P., & Scambos, T. A. (2016). A decade of progress in observing and modelling Antarctic subglacial water systems. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 374(2059), 20140294. https://doi.org/10.1098/rsta.2014.0294

Goff, J. A., Powell, E. M., Young, D. A., & Blankenship, D. D. (2014). Conditional simulation of Thwaites Glacier (Antarctica) bed topography for flow models: Incorporating inhomogeneous statistics and channelized morphology. Journal of Glaciology, 60(222), 635–646. https://doi.org/10.3189/2014JoG13J200

Graham, F. S., Roberts, J. L., Galton-Fenzi, B. K., Young, D., Blankenship, D., & Siegert, M. J. (2017). A high-resolution synthetic bed elevation grid of the Antarctic continent. Earth System Science Data, 9(1), 267–279. https://doi.org/10.5194/essd-9-267-2017

Hoffman, A. O., Christianson, K., Shapero, D., Smith, B. E., & Joughin, I. (2020). Brief Communication: Heterogenous thinning and subglacial lake activity on Thwaites Glacier, West Antarctica [Preprint]. Glaciers/Subglacial Processes. https://doi.org/10.5194/tc-2020-80

King, E. C., Pritchard, H. D., & Smith, A. M. (2016). Subglacial landforms beneath Rutford Ice Stream, Antarctica: Detailed bed topography from ice-penetrating radar. Earth System Science Data, 8(1), 151–158. https://doi.org/10.5194/essd-8-151-2016

Kyrke-Smith, T. M., Gudmundsson, G. H., & Farrell, P. E. (2018). Relevance of Detail in Basal Topography for Basal Slipperiness Inversions: A Case Study on Pine Island Glacier, Antarctica. Frontiers in Earth Science, 6, 33. https://doi.org/10.3389/feart.2018.00033

Le Brocq, A. M., Ross, N., Griggs, J. A., Bingham, R. G., Corr, H. F. J., Ferraccioli, F., Jenkins, A., Jordan, T. A., Payne, A. J., Rippin, D. M., & Siegert, M. J. (2013). Evidence from ice shelves for channelized meltwater flow beneath the Antarctic Ice Sheet. Nature Geoscience, 6(11), 945–948. https://doi.org/10.1038/ngeo1977

Lowry, D. P., Golledge, N. R., Bertler, N. A. N., Jones, R. S., McKay, R., & Stutz, J. (2020). Geologic controls on ice sheet sensitivity to deglacial climate forcing in the Ross Embayment, Antarctica. Quaternary Science Advances, 1, 100002. https://doi.org/10.1016/j.qsa.2020.100002

Malczyk, G., Gourmelen, N., Goldberg, D., Wuite, J., & Nagler, T. (2020). Repeat Subglacial Lake Drainage and Filling beneath Thwaites Glacier. Geophysical Research Letters. https://doi.org/10.1029/2020GL089658

Morlighem, M., Williams, C. N., Rignot, E., An, L., Arndt, J. E., Bamber, J. L., Catania, G., Chauché, N., Dowdeswell, J. A., Dorschel, B., Fenty, I., Hogan, K., Howat, I., Hubbard, A., Jakobsson, M., Jordan, T. M., Kjeldsen, K. K., Millan, R., Mayer, L., … Zinglersen, K. B. (2017). BedMachine v3: Complete Bed Topography and Ocean Bathymetry Mapping of Greenland From Multibeam Echo Sounding Combined With Mass Conservation. Geophysical Research Letters, 44(21), 11,051-11,061. https://doi.org/10.1002/2017GL074954

Morlighem, M., Rignot, E., Binder, T., Blankenship, D., Drews, R., Eagles, G., Eisen, O., Ferraccioli, F., Forsberg, R., Fretwell, P., Goel, V., Greenbaum, J. S., Gudmundsson, G. H., Guo, J., Helm, V., Hofstede, C., Howat, I., Humbert, A., Jokat, W., … Young, D. A. (2019). Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the Antarctic ice sheet. Nature Geoscience, 13(2), 132–137. https://doi.org/10.1038/s41561-019-0510-8

Robin, G. D. Q., Swithinbank, C. W. M., & Smith, B. M. E. (1970). Radio echo exploration of the Antarctic ice sheet. In A. J. Gow, C. Keeler, C. C. Langway, & W. F. Weeks (Eds.), International Symposium on Antarctic Glaciological Exploration (ISAGE) (Vol. 86, pp. 97–115). International Association of Scientific Hydrology. https://iahs.info/uploads/dms/086017.pdf

Rougier, N. P., Droettboom, M., & Bourne, P. E. (2014). Ten Simple Rules for Better Figures. PLoS Computational Biology, 10(9), e1003833. https://doi.org/10.1371/journal.pcbi.1003833

Schlegel, N.-J., Seroussi, H., Schodlok, M. P., Larour, E. Y., Boening, C., Limonadi, D., Watkins, M. M., Morlighem, M., & van den Broeke, M. R. (2018). Exploration of Antarctic Ice Sheet 100-year contribution to sea level rise and associated model uncertainties using the ISSM framework. The Cryosphere, 12(11), 3511–3534. https://doi.org/10.5194/tc-12-3511-2018

Shi, L., Allen, C. T., Ledford, J. R., Rodriguez-Morales, F., Blake, W. A., Panzer, B. G., Prokopiack, S. C., Leuschen, C. J., & Gogineni, S. (2010). Multichannel Coherent Radar Depth Sounder for NASA Operation Ice Bridge. 2010 IEEE International Geoscience and Remote Sensing Symposium, 1729–1732. https://doi.org/10.1109/IGARSS.2010.5649518

Siegert, M. J., Ross, N., & Le Brocq, A. M. (2016). Recent advances in understanding Antarctic subglacial lakes and hydrology. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 374(2059), 20140306. https://doi.org/10.1098/rsta.2014.0306

Sommers, A., Rajaram, H., & Morlighem, M. (2018). SHAKTI: Subglacial Hydrology and Kinetic, Transient Interactions v1.0. Geoscientific Model Development, 11(7), 2955–2974. https://doi.org/10.5194/gmd-11-2955-2018

Wright, A. P., Young, D. A., Bamber, J. L., Dowdeswell, J. A., Payne, A. J., Blankenship, D. D., & Siegert, M. J. (2014). Subglacial hydrological connectivity within the Byrd Glacier catchment, East Antarctica. Journal of Glaciology, 60(220), 345–352. https://doi.org/10.3189/2014JoG13J014

Review of « Mapping high-resolution basal topography of West Antarctica from radar data using non-stationary multiple-point geostatistics (MPSBedMappingV1) » by Z. Yin et al.

This manuscript addresses an important question: the interpolation of sparse data in a non-stationary context. It develops a framework that is seen as generic and demonstrates it in the specific context of interpolating subglacial bedrock topography. The novelty is the use of a large number of training images coming from a variety of deglaciated areas, and to devise a scheme to select a subset of them that is specific to each simulated region.

My overall assessment is that the manuscript is scientifically sound and will have an impact in the glaciology community, and also possibly more generally in other applications of MPS to large-scale problems, including space-time simulation. The tests in section 4 clearly show the value of the approach. Importantly, this study improves our understanding of the usage of MPS in the case of very large models that are fed by lots of data and training images. It demonstrates the need to select a specific subset of training images for each region rather than using the entire training image database for all areas. This strategy provides gains in terms of quality and computation, which is very valuable. Also, it could fuel the discussion on using machine learning approaches to address such problems (e.g. GANs) which would require a new training for each sub-area, making them inefficient.

This said, I do have remarks on some technical aspects of the proposed approach that I outline below. Despite these comments, I find the method convincing and I believe the manuscript deserves publication in GMD:

• While the approach is novel, some of the existing literature was missed. For example, on l.63-54, it is mentioned that this is the first application of MPS to subglacial topography. A recent paper doing just this appeared in The Cryosphere: https://tc.copernicus.org/preprints/tc-2021-161/
• The proposed method involves a number of modeling steps and choices. This is fine when justified, but here I did not see a clear justification for many of the choices made. I detail some instances of this below:
  • I do not see a clear motivation for using a modified Hausfdorff distance in eq. 1. Why using a distance between patterns rather than some statistical similarity between patterns, e.g. in terms of integral scales?
  • There exists a literature on methods to select one or more TIs based on conditioning data (e.g. Pérez, C., G. Mariethoz, and J. M. Ortiz (2014), Verifying the high-order consistency of training images with data for multiple-point geostatistics, Computers and Geosciences, 70, 190-205, 10.1016/j.cageo.2014.06.001, or Abdollahifard, M. J., M. Baharvand, and G. Mariéthoz (2019), Efficient training image selection for multiple-point geostatistics via analysis of contours, Comp. & Geosci., 128, 41-50, https://doi.org/10.1016/j.cageo.2019.04.004). This is the same problem that is addressed here with the probability aggregation approach. Such methods are not considered or mentioned. It is fine that the authors propose a new methodology, but the reason for not using existing approaches should be given.
  • The general approach proposed is rather sophisticated (probability aggregation - PSO optimization – kernel density estimation), which complexifies the implementation. From a user point of view, these steps should be justified. Basically, it is pretty clear how things are done but the why is not always explicit.
  • 303-304: this is a strong assumption because these distances are considered in a high-dimensional space. Can it be justified?
• The procedure for the choice of the TIs is rather complex. This is justified by the dependence between neighboring TIs which is modeled by probability aggregation. However, since there are local data everywhere in the domain, simply selecting a set of TIs statistically similar to the data might also perform well. Has this been tested? For instance, looking at figure 10, I am wondering whether a similar ranking could be obtained with a simpler approach, for instance considering similarity in terms of histogram and variogram similarity.
• It is likely that the training images need to be optimally oriented prior to simulation. Or possibly, one could use a rotation-invariant distance. Has this been tested?
• Section 3.1.3 and the legend of figure 4 are quite unclear to the uninitiated reader and should be improved.
• 227: This statement about boundaries seems incorrect. DS is affected by the boundaries of the TI, as the edges of the TI cannot be sampled, especially for large data events.

Minor comments/typos:

l.22: sentence grammatically incorrect

l.23: provides

l.24: flight lines data

l.86-87: here the work by G. Pirot could be mentioned as it does exactly what is mentioned in this sentence (https://doi.org/10.1016/j.geomorph.2014.01.022, and also doi:10.1002/2015WR017078).

l.96-97: Awkward sentence, rephrase

l.115: remove extra parenthesis

l.149: section 0 is a mistake

l.182: verbatim copy should be explained

l.204: remove word expressed

l.286.287: this reference seems out of place. This is the TI selection, not DS simulation

l.287: erroneous reference: there is no section 3.13

l.292: parameterization … convergence

l.293: make PSO

l.320: “is to reflect” -> should be rephrased

l.370: section 0 doesn’t exist

l.397: rewrite sentence (diving ->dividing)

Figure 18: problem of color scale where white areas appear in the simulation domain

l.424: density of water is 1000 kg/m3

Check the citations: several are duplicated.