Reviewer 1

This manuscript presents an application of the Graphical Subglacial Sediment Transport Model (GraphSSeT) to the Finnish Lake District Ice Lobe (FLDIL), using GlaDS simulations from ISSM as hydrological forcing. The study advances beyond the synthetic test cases of Aitken et al. (2024) by applying the model to a realistic glacial system with well-preserved geomorphological constraints. The analysis of sediment conditions, flux rates, grain size dynamics, and provenance yields interesting results, particularly the strong seasonal cyclicity characterized by fluvial erosion in spring/summer and deposition in autumn, and the extensive sediment mobilization near the glacier terminus.

I have some small(ish) concerns regarding presentation, model coupling and its implications that could be a useful addition, if addressed:



We appreciate the reviewer’s comments with responses in Italics



One-way coupling limitations. The sediment model is forced by GlaDS but does not feed back to the hydrological system as channel and conduit geometry remain unaffected by sediment presence.

The physical formulation in GlaDS uses semi-cylindrical channels with a time-varying diameter. For reasons explained below, we do not believe, in this case, that it is an important factor for accuracy of the model outcomes. 

Given the presence of eskers in the FLDIL, which provide evidence of channel infill, this limitation warrants discussion. Sediment accumulation would reduce the hydraulic radius and effectively choke the drainage system (see Hewitt and Creyts, 2019; doi:10.1029/2019GL082304).

We agree that a sufficiently thick deposition occur within or at the outlet of a channel, both the area and aspect ratio of the channel may be changed significantly. As noted, this would affect hydrological conditions locally and possibly also at network scale. For the model results this is not a significant factor for the following reasons:

1. High enough rates of deposition do not occur in the conditions we model

The peak positive dH/dt of ~1e^-9 m/s yields, over a 13 week “autumn” deposition event, a total thickness deposited of 8 mm or so. Our channels have sectional areas of ~1 m² up to 25 m² (i.e. heights of 0.8 to 4 m).  The implied changes to channel area and shape from sediment deposition and remobilisation are small and not likely to have a significant effect on hydrology.

1. The deposited sediment, having been transported to site in a prior season, is easily remobilised.

With the  forcing applied, conditions are approximately the same each season so we do not get any net buildup over time. This could change significantly for scenarios with a retreating margin and/or waning hydrology forcing however this is not tested here.

1. We do not model the end-point deposition process, which lies outside the model domain.

The outlet edges have a free-flow boundary condition (at transport capacity or rate of supply). By implication, the depositional environment is assumed to accept all sediment discharged without limitation. A further limitation on discharge-rate could, in principle, be applied in GraphSSeT, but we do not do so here.

I recommend the authors expand the discussion to address how the hydraulic system is likely to respond to the modelled sediment dynamics, and how such feedbacks might influence the spatial zonation between erosion, mobilization, and transport. While the authors discuss time-transgressive margin behaviour, the likely implications for hydraulic system evolution would strengthen the manuscript.



With these comments in mind, we have revised the discussion to address the interplay between sediments and hydrology considering a) the validity and significance of the model results b) potentially important processes operating that are not represented in the model.



Temporal resolution. GraphSSeT receives input from GlaDS at four-day intervals, and ice geometry is held fixed. I recommend the authors consider in the discussion what higher-frequency diurnal dynamics in meltwater transport might imply for sediment transport patterns.

The input model (Hepburn et al., 2024) does not include diurnal or daily melt forcing. Melt is derived from an interpolated PDD model, and total monthly melt was routed to the bed via a series of moulins. Sampling the input GlaDS model at a finer temporal resolution would not meaningfully alter our results.

Diurnal or other short-term meltwater dynamics could significantly increase transport capacity with high velocity flow in unenlarged channels. Previous models indicate that diurnal variation can have a net positive effect on sediment discharge, broadly in line with daily peak-flow, see Aitken et al., (2024).

With diurnal cyclicity, the net capacity for sediment transport may be increased, and the effective area susceptible to transport increased. The seasonal system may be affected (e.g. earlier onset of transport) and the sediment transport system may reach further inland.

Despite this, all the model scenarios in this study rapidly become supply limited, and a major impact on net transport on multi-annual timescales is not expected. We include a brief comment on the expected effects of including diurnal variation.

Sensitivity to transport initiation (possibly out of scope for this paper). The inverse relationship between peak sediment flux and grain size is intriguing, explained as a consequence of the upper drainage system mobilizing high volumes of finer sediment. In this context, the results are likely sensitive to parameterization at the onset of sediment transport under relatively low flow velocities. Some discussion of this sensitivity would be valuable.

This is a valuable point and there is important work to be done here with more sophisticated sediment transport laws, incorporating shielding, improved physics for coarse and fine grain sizes and multimodal distributions. Such work is on the agenda for GraphSSeT but is not currently available.

The Engelund and Hansen (1967) formulation is empirically derived for sands, and the ‘transportability’ of sediment with grain size is represented only by d50. In GraphSSeT, many mixing operations are performed, with possibly diverse and/or multimodal distributions resulting. To avoid propagating degrees of complexity, from these mixtures the new d50 and a unimodal distribution are calculated.

Consequently, the results of our model do not contain enough information to comment with precision on mechanisms of initiation or the details of differing transport rates for different grain sizes which are lost in the mixing.

We have observed in earlier synthetic model testing, see Aitken et al. (2024) - that a batch of finer-than-normal sediment will cause an enhanced discharge to propagate downstream, while coarser-than-normal sediment batch will have the opposite effect.

We have improved the discussion here around several points that warrant attention

• In the context of increasing water flow, we infer, but cannot demonstrate, that the observed drop in mean grain size is due to the mobilisation of a (randomly occurring) finer-grained component, rather than the deposition of a courser-grained component.
• With limited sediment availability, sediment will survive longer in areas of the model where transport is only possible for finer grain sizes. In supply limited conditions this will exacerbate grain-size selectivity – which is seen in the relatively muted grain size variation in the Mixed Bed models.
• The d50 values in this work are dominantly focused in the fine sand range however in cases they do (locally) reach down to very fine sands, which is approaching the lower limit of applicability for Engelund and Hansen (1967).

 

Specific Comments



- L2: I suggest adding that the FLDIL is interpreted as a re-advancing lobe in the Fennoscandian Ice Sheet during the Younger Dryas.

Yes, we can add something like this: “...to model subglacial sediment processes. The Finnish Lake District Ice Lobe (FLDIL) is a major ~500 km long ice stream of the Fennoscandian Ice Sheet that has been interpreted as being oscillatory and re-advancing during the Younger Dryas. The FLDIL presents an opportunity to...



- L26-27: Integration of subglacial hydrology in ice-sheet models is not recent. However, recent developments include more sophisticated hydrology modeling, like GlaDS with a two-component drainage system and discrete channels.

Edited to reflect the significant recent progress in this building on the long history

- L160-169: From this description, I do not fully comprehend how the initialization run is set up. Obviously, the model needs sediment to be present or added somewhere in order to mobilise it in the drainage system during the initialization.  Hepburn et al. 2024 ran FLDIL-GlaDS for a hard-bed geometry, where topography was constructed by removing sediment (L134-135).  For the initialization of the three scenarios with GraphSSeT, did you apply a uniform initial sediment thickness (low and high-H) for the entire area, before letting GraphSSeT evolve it under a winter-state from GlaDS, or something else entirely? If so, how do the two scenarios compare to the volume removed in Hepburn et al., 2024 (i.e. the volume in the mixed-bed case)?

Revised along with the preceding text to make clear that for the low-H and high-H scenarios the sediment begins with a small initial thickness (5±2.5 cm), and is dominantly sourced from erosion. Erosion potential is constant through the model (see line 153) in line with the constant velocity used in Hepburn et al (2024).

For the mixed bed model, we effectively ‘add in’ the sediment that was removed in Hepburn et al., (2024). So as to maintain consistency of forcing across the scenarios we do not adjust the bed elevation. This leads to some discrepancy with the real situation, but the main point of this model is to show the impact of increased supply, not the small change in bed elevation.

Sensitivity testing done in Hepburn et al., (2024) showed that this change in elevation from removing the sediments does not drive major change in the hydrology system (e.g., see Figures A22—A23 in the Appendix of Hepburn et al., 2024).

- Fig. 3: I think it would be an aid to the reader if the labels of panel a are also included in panels b-d, i.e. default 1, default 2, detrital, high-sigma, as well as Mixed-bed, High-H and Low-H. Currently, the eyes have to find the corresponding label for a given run by running from each panel with the physical result and the ensemble tree structure, and it is easy to get lost during the travel.

We now include these labels in Fig 3



- Fig. 3: I'm a bit lost in the explanation between parameterization between fig. 3, L.174-181. I do not understand why there are two runs for each configuration, i.e. there are two "Low-H + detrital" nodes, etc. Also, what is the difference between "default 1" and "default 2"? L180-181 mention a cascade of runs for each leg, does this imply some kind of perturbation since the results are not equal?

GraphSSeT is stochastic with respect to grain size, and several other model behaviours,  and so multiple runs are made to indicate the intrinsic variation under identical conditions, before we vary parameters. We have tried to explain this logic more clearly



- Fig. 3d: Is the proportion of basal sediment referring to the proportion of sediment in active transport?

We have reversed this definition for clarity – bedrock-derived is sediment that has never been in the basal sediment layer, it comes straight from erosion of bedrock. We have explained this more clearly in the caption also we updated to bedrock instead of basement.



- L176: I suggest adding the explaining text (high-sigma) to mode b, so it is easier to understand the equivalent modes explained in text and table.

Yes, we do this



- Fig. 4: Consider making the "Grain size" label text on the secondary y axis red, as a visual aid. Would it be meaningful to include the channel water flux also, to show sediment starvation?

Grain size colouring is edited as indicated.

We refrain from directly including the water flux at this stage and adding a third axis to the plot, but we do change the x-axis to ‘absolute’ years to allow easier comparison with Fig 2a and we mark on the times of water flow maxima and minima.

- Fig. 4 (cont.): Is it correctly understood that the basal sediment flux is the locally eroded and mobilized sediment along the channel (edge), and the total sediment flux is the local contribution and upstream flux combined? I am a bit confused about the "basal" terminology, as I understand the Engelund and Hansen 1967 to account for both bed-load and suspended sediment transport.

The basal vs bedrock indicates whether the sediment came directly from erosion or from mobilisation of basal sediment. As in Fig 3d this point is made more clear in the caption to aid understanding  



- Fig. 4b and c: Could you please consider inverting one of the colorbars? In panel b, black is low and yellow is high, and the opposite is the case in panel c.

We have reversed the scale for the standard deviation in all plots, to be light colours up



- L190-205: These findings are very interesting, in particular the drop in grain size during peak discharge. Did you try and visualize the spatial patterns in grain size in these figures (4 and 5) or in the transient maps (fig. 7)?

The spatial patterns in grain size are recorded in the supplementary material for all models, showing both the active sediment and the basal layer. These do show significant and systemic variations. I recognise the link was not accessible for reviewers due to a Zenodo problem.

While the grain size is a little harder to interpret, this is an important model behaviour pointed out by both reviewers, and so we have added a figure showing d50 in the active layer to complement Figure 7 (the new Figure 9).

This manuscript applies a coupled ice-sheet, hydrology, and sediment transport model to reproduce the geomorphological features of the Finnish Lake District Ice Lobe (FLDIL) utilizing the Graphical Subglacial Sediment Transport model (GraphSSeT). The GraphSSET model evolves based on subglacial hydrology inputs from the hydrology model GLaDS coupled with ISSM. The manuscript highlights the strong seasonality of the sediment system across multiple model ensembles under initial sediment/bed conditions. In particular, coarse grains are transported through the system in the early spring and then fine sediment dominates the summer melt season. These results show that coupling sediment transport with subglacial hydrology produces sediment deposition patterns that resemble geomorphological features consistent with the ice-marginal deposits of Salpausselkä II. 

The authors address a central challenge in glaciology: reproducing observed glacial landforms within a process-based modeling framework. By linking sediment transport dynamics informed by subglacial hydrology with preserved geomorphological features, this study provides insight into how subglacial hydrology can influence the formation of ice-marginal features. This manuscript demonstrates how intertwined processes, such as seasonality of subglacial hydrology, ice dynamics, and topography can inform landscape evolution. The research presented here would be of interest to those who study sediment in subglacial systems, glacial geomorphology, and those seeking to use numerical models to interpret or reconstruct past ice sheet behavior. Given this broad interest, I believe the manuscript is well suited for publication in The Cryosphere, with the minor revisions suggested below. 

We thank the reviewer for their review and comments 

I have a few suggestions on how to improve this manuscript, in particular, where to expand on certain modeling choices.

1. As a reader, I would have appreciated a short summary or schematic of equations/variables that show the sediment and hydrology coupling. This manuscript follows up the idealized GraphSSeT model publication (Aitken et al. 2024) by applying the GraphSSeT model to the FLDIL. Readers are referred to the original manuscript that does include a schematic, but I suggest the coupling between GraphSSeT and GlaDS (named FLDIL-GLaDs simulation in line 133) can be briefly described with a new schematic that shows the application of these models to this system. In the appendices, important parameters discussed in the results (e.g., volume, mean flux, mean grain size, and proportion of basal sediment) vary within the model time. A schematic that shows how these variables are derived within the coupled sediment-hydrology model would be appreciated. 

To address a potential point of confusion, the ‘coupling’ is the same as it was in the earlier study – a one- way use of GlaDS outcomes to force sediment transport (see also comments above).  This application is of course much larger and more complex than the previous study, but the process is identical.  

We add a summary figure showing some key model behaviours, focusing on the network scale dynamics and seasonal change rather than the more basic graphic in Aitken et al., 2024.

With respect to the variables and names we have made an effort to define these more clearly through the text.

1. I suggest the authors discuss seasonality in the Younger Dryas. The FLDIL-GLaDs simulation assumes a modern precipitation record, a depressed modern temperature record, and an elevation dependent lapse rate of 7.5 C/km. This climatology is then used to force a hydrology model with a timestep of four days, and a sediment model at an even finer timescale. Thus, I imagine that uncertainties in the starting climatology might influence the results. How sensitive are the results to the prescribed climatic conditions? What is known about seasonality in this region during the Younger Dryas? I encourage the authors here to engage with previous work suggesting that seasonality was different from modern during the Younger Dryas (e.g. Denton et al., 2005; Hall et al., 2008). 

To avoid confusion we must be clear that the GlaDS timesteps were adaptive and typically much less than 4 days, however this is the frequency with which results were stored. In Hepburn et al., 2024, GlaDS was run with an adaptive timestep which was allowed to vary

For the current case-study we restrict commentary to the modelling results of Hepburn et al. (2024). It is noted in the text that this is one realisation and refer readers to the existing discussion in Hepburn et al (2024) on the climate forcing applied.

1. I suggest the authors discuss the time dissonance with 10 years of model time aiming to replicate geomorphic features that take much longer to form. The manuscript presents spatial agreement between the modeled sediment distribution and the observed geomorphic features of Salpausselkä II. However, the authors note (lines 321–331) that the rates of discharge in most of their models would take much longer to form these features than implied by the geology. I encourage the authors to expand on this discussion of the temporal framework of the simulations. The experiments are conducted over 10 years of model time yet the geomorphic features being replicated potentially require substantially longer timescales to form (1000-10000 years). This creates a potential time-scale dissonance that can be expanded on in the discussion. Specifically, what does a 10-year simulation meaningfully inform about features that develop over millennial timescales?

Notwithstanding limitations from the compute resources and the researcher’s patience, our work is driven by the timescales of subglacial hydrology and fluvial sediment transport – which are relatively short – and the time needed to development a dynamic balance between gains and losses. 

For the ‘hard bed’ low-H and high-H scenarios a consistent discharge rate is reached within 5 years. Longer durations, without changing the forcing, would yield more cycles, but no significant change in discharge rate. For the mixed bed scenario, the persistence of supply means that we do not reach a balance, but output volumes are declining only slightly (see Fig 6). Longer durations with this model may see ongoing slow decline until supply is exhausted.

Longer-term evolution can occur related to the sustainability of sediment supply, and processes such as glacial advection are significant on longer timescales. However, all of these are likely to be trumped by the changing ice dynamics and the impact on hydrology (see the next comment)   

This connects to another question that I had in the context of seasonality - even on short timescales the ice margin might advance and retreat, and during the younger dryas we know there were phases of ice sheet growth and decay. It’s fair enough that for this study the authors use a fixed ice sheet geometry, but I would appreciate a deeper discussion of how they think the results might change in the future were a transient geometry to be considered. For example, the discussion outlined in L305-319 regarding implications of the modeled sediment distribution against observations is expanded in a way that discusses such limitations.

We include some further comments in this section with respect to alternative margin evolutions. However, we note that these have not been modelled and we cannot know what the model would show.

We reiterate our tow main points: that for the Salpausselkä II, erosion cannot yield enough sediment in the formation timeframe – erosion rates would need to be ~; and that hydrologic transport can generate this sediment very quickly from pre-existing sediment cover. Therefore, the important factor is the distribution of this cover relative to the margin.

A retreating or oscillating margin location would be able to tap into different sources of sediment and could yield a faster deposition rate.

1. I suggest the authors replot the spatial results in Figures 4-6 (b-c) with a focus closer to the ice margin. The changes occurring in the modeled interior do not appear to be the primary focus of the results. A zoomed-in view of the model domain near the ice margin would likely improve visualization of the spatial patterns and make the color scale easier to interpret. The model extent shown in Figure 7 may provide a useful domain for the subplots in Figures 4-6. 

We now show the zoomed in view in all figures except Figure 1 (now Figure 2). Noting that the zoomed out view is accessible through supplementary material

Specific comments:

L29: A reference to Sommers et al., 2022 (ISSM-SHAKTI application at Helheim Glacier, Greenland) is also appropriate here when referring to subglacial hydrology models coupled with ISSM.

 Indeed, I have included this reference, also noting that GraphSSeT does not specifically require GlaDS hydrology models.

L57-L62 and Figure 1: It may be beneficial to add and label Salpausselka I and II on Figure 1a, whether as an average line or arrow. 

Labels are added 

Figure 1: A label for the model domain in black would be useful to add either to the legend or figure caption.

Added to legend 

L110-114: A figure or schematic showing how GLaDs and GraphSSeT are coupled would be useful here (see point 1 above).

See response above

Figure 3d - a gradient color bar would be suitable rather than a divergent color bar. Additionally, is the percentage calculated here = basal sediment divided by basement sediment? I think a clarification in the figure legend would be helpful.

Agreed and updated

Table 1 - include a column in Default 2 parameters. This seems important in Figure 3, but is not mentioned in this table.  

Default 2 uses the same parameters as default 1 so I edit the caption

L183: Are you referring to Figure 3b or 3d for the variable volumetric discharge? A clarification on which sub-figure is useful here. 

Figure 3b, this is now specified

Figure 9: For the sediment classifications outlined in the legend, I would appreciate it if they can connect back to the legend in Figure 1b. Is metasedimentary all sedimentary classifications in the Figure 1b legend? A connection back to the region can be helpful here.

The bedrock classes link directly Fig 9 to Fig 1 – so we have made the colours and names the same. We also make the basal and eroded sediment more obviously different colours. Unlisted volume elements are not observed in the detrital output.