Based on ocean-circulation modelling, this paper examines circulation in proglacial fjords where the exchange flow is primarily driven by buoyancy forcing from subsurface ice melt and subglacial discharge. More than seventy numerical experiments have been conducted to investigate how the fjord circulation depends on conditions such as sill height, ambient water temperature and stratification, and tidal flow. The paper offers many interesting results, including how high sills can cause a transition to hydraulically-controlled exchange flow, which induces cooling and recirculation in the waters between the sill and the glacier. However, some of the main findings could be presented in a more general and accessible way; particularly making limitations of the results clearer when applying them across seasons and between oceanographically different fjord systems.

• We appreciate the reviewer’s careful and thoughtful consideration of our manuscript and the many suggestions for improvement. Our responses to each of them are detailed below.

Main comments

1. How is the model stratification set?

It is not clearly described how the stratification is set in the model in section 2.1. Are the authors restoring to the idealised Greenlandic profile in the open-ocean part of the model domain? Or are they initialising the model with a stratification that is allowed to evolve in the simulations? This issue is highly important for the interpretation of the model results; it decides if the simulations yield fjord stratifications that are determined by the interplay between melt dynamics, the sill, and the open ocean conditions. (The results in Fig. 7, for example, make me suspect that the stratification is set by the initial conditions.) Explain this clearly.

• The model is initialized with a prescribed stratification that is then allowed to evolve.  Initial temperature and salinity profiles are indeed restored at the open boundaries throughout the simulation. But because there’s a shelf outside the fjord (set to 25 km in the cross- and along-shelf directions), restoring the conditions does not change our results. Prompted by this comment, we ran several simulations with a wider (60 km) shelf and confirmed this change did not impact our results. That said, the initial conditions do influence our analysis as the initial stratification and temperature help set the initial melting rate, the height of the plume, etc. 
  
  
  
   We clarified the initial fjord conditions setup as follows in the revised manuscript:  “Initial fjord conditions are horizontally homogeneous, with temperature and salinity profiles restored only at open boundaries throughout the simulation. Tests expanding the shelf from 25 km to 60 km showed these boundary conditions did not impact our results in the fjord. Initial water temperature is a constant of 2 to 10 °C. The vertical salinity gradient is based on an idealized Greenland fjord profile (Cowton et al., 2015) where the salinity ranges from 32 to 33.8 in the upper 80 m and slowly increases to 34.5 at the bottom.”     

2. The sill height

The authors use the ratio between the sill depth (h s) and the fjord depth (h f) as a measure to distinguish/discuss flow regimes, and also refer to effects of hydraulically-controlled exchange flows. However, in a two-layer description of hydraulic flows (see e.g. Pratt and Whitehead, 2007) it is only the upstream height of the layer interface above the sill crest that matters for the dynamics — the fjord depth does not enter. The authors need to expand on this matter and discuss how the specific T and S profile they use as initial conditions (or restoring open ocean conditions?) relates to the impact of h s on the flow. [L102: here you should describe the idealised Greenlandic salt stratification and estimate an approximate two-layer representation of the vertical density distribution; with an interface depth of say h i : Note that h s − h i and the layer density difference are key variables determining the flow characteristics (see e.g. Pratt and Whitehead, 2007; Schaffer et al., 2020; Nilsson et al., 2022).] Results of for instance Jakobsson et al. (2020), Schaffer et al. (2020), and Nilsson et al. (2022) show that hydraulic control can emerge in North Greenlandic fjords with marine glaciers that have relatively deep sills (h s ≈ 400 m) where h s /h f ≈ 0.5. It could be relevant to mention that the present study to some extent also reveals impacts of sill geometry on the exchange circulation and reflux in such fjord systems.

• We agree with the reviewer that the studies mentioned are very relevant, and we’ve added references and text acknowledging the potential for hydraulic control in deep-silled fjords. Those systems are not the focus of our paper, however, and a complete description of all possible hydraulic control scenarios, which can be rather complex in the presence of tides and multi-layer flows, is beyond the scope of our work. The text reads:
  
  "Hydraulic control can also occur in deep-silled systems (Pratt and Whitehead, 2007), where this process can also play a significant role in modulating deep water properties and the heat supply to marine-terminating glaciers in fjords (Schaffer et al., 2020; Nilsson et al., 2023)..."
  
  
  
  We should note that the depth of the fjord (and the ratio to the sill height) is a relevant parameter in our framework of analysis because the initial meltwater plume height might or might not reach the surface of the fjord or the sill.
  
  The interface depth at the sill/constriction (rather than upstream as in the 11/2 layer model used in the Pratt and Whitehead as applied in the papers referenced above) is our preferred approach here as it wasn’t self-evident to us the upper layer couldn’t also be hydraulically controlled (that is, this is something we did not want to assume), and certainly we couldn’t assume that the upstream interface depth is much larger than the sill depth.
  
  
  
  Here we are following dynamics that are well established in prior literature (e.g., Stigebrandt, 1981; Farmer and Freeland, 1983, Geyer and Ralston, 2011), and thus we use the properties over the sill to define, for example, the relevant Froude number, which is what we use to diagnose hydraulic control. We do not use the ratio of sill depth to fjord depth for this purpose. 
  
  
  
  We thank the reviewer for bringing these recent papers on deep-silled dynamics to our attention. 
  
  
  
  References:
  
  Stigebrandt, A., 1981. A Mechanism Governing the Estuarine Circulation in Deep, Strongly Stratified Fjords. Estuarine, Coastal and Shelf Science 13, 197–211. https://doi.org/10.1016/S0302-3524(81)80076-X.
  
  Farmer, D.M., Freeland, H.J., 1983. The Physical Oceanography of Fjords. Progress in Oceanography 12, 147–220.
  
  Geyer, W.R. and Ralston, D.K., 2011. 2.03—The dynamics of strongly stratified estuaries. Treatise on Estuarine and Coastal Science. Amsterdam: Elsevier, pp.37-52.
  
  
  
  
  
  

3. Submarine melting

If the flow is hydraulically controlled, a taller sill is expected to diminish the ocean heat transport towards the glacier, and hence to reduce the subsurface melt (Schaffer et al., 2020; Nilsson et al., 2022). A puzzling result of Table 3 is that the shallowest sill experiment yields slightly higher subsurface melt that the ”no sill experiments”, with the highest melt found for hs /hf = 0.12. Can you explain this? (Are the experiment transient in character and do not give equilibrated melt rates?) Additionally, could some general information in Tables 2 and 3 be extracted and represented graphically in a figure? I find it difficult to digest the results in the tables.

• The melting rate depends on both the deep fjord temperature and the stratification. Sill processes modify both, but in opposite directions (reducing stratification, which promotes a ‘taller’ freshwater plume and enhanced submarine melting, but also cooling, which reduces melting. We have revised the related sections to provide what we hope is a clearer explanation of how these changes can result in enhanced submarine melting despite the fact that the inflow is cooled.
  
  
  
  In response to the second reviewer’s comments, we also decided to remove the discussion of Table 2 and Table 3 and focus instead on the impact of shallow sills on the variables impacting melting directly, i.e., stratification and deep water temperature. 
  
  
  
  With constant subglacial discharge in a linearly stratified fjord, the dependence of submarine melt on fjord temperature and stratification can be scaled by (T_a-T₀)(N²)^-5/8, where T_a is the ambient fjord temperature, T₀ is the freezing temperature of seawater. As shown in our new plot (Fig. R1, see attachment), the effect from stratification is more pronounced than that from fjord temperature with the presence of a shallow sill, resulting in fjord cooling but increased submarine melting.
  
  
  
  You can see that in several cases the combined cooling and destratification effect largely cancel each other out (i.e. it results in motion along an “isomelt” in the diagram) while in some others (particularly the warm, low stratification cases) the cooling effect dominates, creating significant change in melting. 
  
  
  
  

  References:
  
  Slater, D.A., Goldberg, D.N., Nienow, P.W. and Cowton, T.R., 2016. Scalings for submarine melting at tidewater glaciers from buoyant plume theory. Journal of Physical Oceanography, 46(6), pp.1839-1855.

4. Unsteady flow regimes

The authors present four flow regimes, two of which are unsteady. In regime III, freshening due to subsurface melt will continuously increase the buoyancy of the fjord water below the sill level. Also in regime IV, the subsurface melt increases the buoyancy of the fjord water, and at the same time the exchanges flow transports buoyancy into the fjord. In both cases, the fjord will convectively overturn after sometime, establishing a circulation in the regime I or II. The question is after how long.

On L345, the authors state that their results suggests that Regime IV may persists on seasonal timescales. There seems to be very little support for this statement. Furthermore, the lifespan of transient regimes like III or IV can presumably vary greatly, depending on the particular fjord system and what processes that forced a transition into the transient regime. This needs to be discussed and better quantified.

• We appreciate the reviewer’s suggestion to include additional support in the unsteady regime timescales. In Regime III, the fjord basin may act as a closed control volume or a ‘‘filling box’’ [Baines and Turner, 1969], where the outflowing plume is blocked by the sill and progressively fills the basin downward from the initial level of neutral buoyancy. In a linearly stratified environment such as our cases, the time scale ta of the ascending plume can be expressed as:
  
  t_a=0.12γ^-4/3H₀^-2/3AB^-1/3τ
  
  
  
  where H₀=0.41γ^-1/2B^1/4(N₀²)^-3/8 is a characteristic length scale (Morton et al., 1956) proportional to the initial plume height h_p in Eq. (6). A is the horizontal cross-section area from glacier front to sill, B=g₀^'Q_sg is plume buoyancy flux. Following the ‘stratified filling box’ model developed by (Cardoso & Woods, 1993), the non-dimensional time can be predicted by
  
  τ=2^-7/3[(h_f-h_s/H₀)²-5.76]
  
  
  
  For initial stratification 1N₀²~4N₀², theory gives t_a≅7, 20, 33, 44 days, as compared to the model output (Fig. R2).
  
  
  
  For Regime IV, the sill-level outflow may be still partly blocked, filling the fjord basin progressively as in Regime III. This scenario resembles the summer circulation and stratification regime at LeConte Bay, Alaska (Hager et al., 2022), it may last for several months before the decline of subglacial discharge and associated meltwater as well as collapsed stratification in fall and winter. If the exchange flow is fully reversed above the sill like a reverse estuary (Giddings and MacCready, 2017) the circulation can be impacted by variability (usually with time scales of days to months) outside the fjord, such as density variations originating from along-shore winds or density anomalies advected past the mouth of the fjord (Straneo & Cenedese, 2015).
  
  
  
  Our model does not resolve the shelf variations and anomalies, therefore the transitory time for Regime III or IV is expected to be mostly modulated by the fjord-sill geometry and properties of the turbulent plume. For the 60-day experiments, we didn’t see the plume outflow reaching the fjord surface in the shallowest-sill and strongest-stratification case. And based on the real-world study mentioned above, the unsteady flow regimes could last a few months. This would obviously depend on the bathymetry and other properties of the freshwater inflow.
  
  
  
  References:
  
  Baines, W.D. and Turner, J.S., 1969. Turbulent buoyant convection from a source in a confined region. Journal of Fluid Mechanics, 37(1), pp.51-80.
  
  Cardoso, S.S. and Woods, A.W., 1993. Mixing by a turbulent plume in a confined stratified region. Journal of Fluid Mechanics, 250, pp.277-305.
  
  Giddings, S.N. and MacCready, P., 2017. Reverse estuarine circulation due to local and remote wind forcing, enhanced by the presence of along‐coast estuaries. Journal of Geophysical Research: Oceans, 122(12), pp.10184-10205.
  
  Straneo, F. and Cenedese, C., 2015. The dynamics of Greenland's glacial fjords and their role in climate. Annual review of marine science, 7, pp.89-112.
  
  Hager, A.O., Sutherland, D.A., Amundson, J.M., Jackson, R.H., Kienholz, C., Motyka, R.J. and Nash, J.D., 2022. Subglacial discharge reflux and buoyancy forcing drive seasonality in a silled glacial fjord. Journal of Geophysical Research: Oceans, 127(5), p.e2021JC018355.

Language

Please improve the language: there are grammatical mistakes and some formulations that are a bit unclear

• Checked and improved.

Minor comments

L1: perhaps change ”of glaciers is” to ”of ice sheets are” (since glaciers are not main contributors to sea level).

• We appreciate the reviewer’s attention to detail. To further clarify, we have cross-checked our statement with relevant studies (e.g., Pörtner et al., 2019; Hugonnet et al., 2021). In the early 21st century, the mass loss of glaciers contributes approximately twenty percent of the observed sea-level rise. Particularly, thinning rates of glaciers outside ice sheet peripheries doubled over the past two decades. Glaciers currently lose more mass, and at similar or larger acceleration rates, than the Greenland or Antarctic ice sheets taken separately.
  
  
  
  References:
  
  Pörtner, H.O., Roberts, D.C., Masson-Delmotte, V., Zhai, P., Tignor, M., Poloczanska, E., Mintenbeck, K., Alegría, A., Nicolai, M., Okem, A. and Petzold, J., 2019. IPCC special report on the ocean and cryosphere in a changing climate. IPCC Intergovernmental Panel on Climate Change: Geneva, Switzerland, 1(3), pp.1-755.
  
  Hugonnet, R., McNabb, R., Berthier, E., Menounos, B., Nuth, C., Girod, L., Farinotti, D., Huss, M., Dussaillant, I., Brun, F. and Kääb, A., 2021. Accelerated global glacier mass loss in the early twenty-first century. Nature, 592(7856), pp.726-731.

L8: ”leads to 10% cooling”. Obviously, the cooling rate depends on reflux as well as the temperature difference between inflowing and outflowing waters. Thus, ”leads to 10% cooling” needs to be related to actual temperatures or a specific fjord type; e.g. an Alaskan fjord in summer.

• 1°C decrease in fjord temperature, we modified the wording in the abstract.

L103: Are you restoring to the idealised Greenlandic profile in the open-ocean part of the model domain? Please explain.

• We explained it in the main comments.

L119 and Table 1: The authors are analysing many cases, and it may be helpful if they define a reference case, meant to characterise a particular fjord (or group of fjords). The authors use a Greenlandic salt profile, and mention glaciers in Patagonia and Alaska. I note that in most of the experiments, the temperature is 10 ◦ C. This is much warmer than subsurface Atlantic Water temperatures around Greenland (Straneo et al., 2012). So to present one reference case would be helpful.

• The Greenlandic salt profile is used as a reference profile, but it is then modified to explore the parameter space, including using a weaker temperature forcing that more closely resembles shelf conditions off Greenland (see Fig. R1).

L123: Is subglacial discharge values used here small or large for a 2 km wide fjord? I would expect that a subglacial discharge of 1000 m s −1 into a 2 km wide fjord is a bit extreme; or even unrealistic?

• We agree this is a somewhat extreme example, which we included to have a fuller description of the parameter space. The purpose is to cover a full range of very weak to very strong subglacial discharges, although we mostly focus our analyses on the case using 250 m³/s.

L174: ”driven only by subglacial discharge”; I suppose you mean driven by the temperature forcing and subglacial discharge.

• Correct, the text is modified and it is the buoyancy-driven circulation fed only by subglacial discharge.

Fig. 3: mention the temperature of the experiment.

• T_ini = 10℃, added to the figure caption.

Table 3: I assume that ”h s /h” should be ”h s /h f ”. Also, why do not Q 0 f and Q 0 s balance each other? Is this due to the plume parametrisation?

• Yes. The divergence between Q_i^f and Q_i^s corresponds to the freshwater discharge, or Q_sg based on the conservation of mass.

Figure 7 and Eq. (6): I repeat that I don’t understand how the salinity stratification is set or prescribed in the model (see main comment and L103 above). If the stratification would have been restored in the open ocean, then I don’t see how this could affect the stratification below the sill level in the fjord. Explain what is going on.

• The application of Eq. (6) and results in Fig. 7 are for initial fjord stratification, which is allowed to evolve. We do restore the properties at the boundaries, not everywhere on the shelf. Our sensitivity experiments show that the boundary far-field shelf conditions do not impact our results significantly. 

L420: ”With a sill depth of h s /h f = 0.04, about 70% of the plume-driven outflow is refluxed to depth.” As stated in the major point, h s /h f is essentially irrelevant for whether the flow is hydraulically controlled, and outflowing glacially-modified water is entrain into the inflowing oceanic water.

• We do not diagnose hydraulic control from hs/hf, but rather from the calculations of the layers Froude numbers at the sill, which does show this quantity reaching or exceeding 1 for the shallower sill cases. We should also note that “reflux” quantifies all processes transporting water from the upper to the lower layer, even though as the reviewer points out one expects that transport to be greatly enhanced when there’s strong mixing associated with hydraulic control. 

This manuscript explores fjord circulation in proglacial fjords with an extensive set of numerical simulations. The study focuses on the impact of a sill in driving reflux _ an exchange from the outflowing layer back into the deeper inflow, which impacts the fjord circulation patterns, heat content, stratification, and – in some way – the heat transport to the glacier that drives melt (though as outlined below, this last point is the weakest part of the paper). The authors identify 4 different regimes based on the relative depths of the fjord, sill, and plume, which is an important step towards mapping out the parameter space of different fjords around the globe. Overall, the paper is well-written with clear figures and interesting results. While I have several major comments, I think these are addressable and this paper is well on its way to being a valuable contribution to the literature on ocean-glacier interactions in fjords.

• We thank the reviewer for this very constructive and thorough review. We believe that in addressing them the manuscript has improved significantly. 

Major Comments 

1) Heat transport, heat budgets & melting. In my view, the parts of the paper requiring most attention (or cutting) are the aspects related to heat transport and melting. The paper tries to evaluate parameters that affect the “heat transport to the glacier” and thereby submarine melt rates. But there seem to be multiple, intertwined points of confusion about heat budgets and melt rates.

The primary “heat transport” calculations are reported in Table 2, as the “heat flux” in the upper or lower layers at various sections. These heat flux calculations – or their relationship to the total heat budget – are not described in the Methods section, so it’s hard to follow exactly what has been done, but presumably it’s a simple calculation of velocity x temperature over the layer. There is one sentence at the end of the Methods sections saying that they “construct the volume and heat budgets within the fjord” following Jackson & Straneo (2016). However, the components of a full budget are not presented, and instead the results are reported in terms of heat fluxes over different layers of a partial cross section. (Also, there is some ambiguity with different control volumes mentioned, but a complete heat budget is not addressed for any control volume, so that’s a secondary concern.)

One issue with the results, as reported, is that the heat flux through a transect with a net mass transport doesn’t have a clear meaning – it’s the heat flux divergence between all bounding surfaces of a control volume that has a well-defined meaning (see, e.g. Montgomery 1974 or Schauer & Beszczynska-Moller, 2009). For example, the text says that, near the glacier, the heat flux in each layer H^f_0 and H^s_0 is smaller with shallower sill, then attributes (without explanation) this to cooler temperatures. But given how these quantities are calculated, couldn’t this reduced heat flux be partially or entirely from reduced volume flux? (It’s stated that Q_0 also decreases with reduced sill). A comparison of Table 2 and 3 suggestions that heat flux through a layer is largely proportional to the volume flux. This gets at the fact that the heat flux/transport through only one boundary of a control volume does not have meaning in an absolute sense. These calculations of heat transport across partial sections (i.e. over certain layer) will intertwine the temperature of the transport with the volume transport. Relatedly, the heat flux is referred to as the “heat supply to the glacier” – this seems problematic because the heat transport *towards* the glacier in the lower layer largely just feeds into the heat transport away from the glacier in the upper layer. The net heat flux *to* the glacier is a tiny difference between the inflow and outflow.

If heat transports are to be calculated, I would encourage the authors to put these in the context of a total heat budget for the control volume. If you are trying to make heat transport calculations that are relevant to submarine melting, I’d consider defining a control volume between one of the transects and the glacier. Then, the heat flux divergence (i.e. the difference between inflowing and outflowing heat) will go to [1] changing the T of control volume waters, i.e. storage term, and [2] heat for submarine melting. For a control volume that is contained mid-fjord, heat flux divergence should just balance a heat storage term.

• We appreciate your feedback and attention to the heat flux/transport calculations presented in our study. We agree that a full discussion of the heat budget terms might distract from a more focused presentation of the impact of sill processes on the stratification and temperature, which can readily be related to submarine melting. In the revised manuscript we have removed the heat budget tables and results in favor of the key parameters (subglacial discharge, fjord temperature, and fjord stratification) that determine the modeled glacial melt.

When it is claimed that certain cases reach steady state, can you quantify this? E.g. around L296-298, for the “steady” case, it’s claimed that the heat exchange compensates for heat loss due to mixing and melting… but did you quantify the heat storage? Is the storage term actually small compared to heat for melting? Here and elsewhere, it would be helpful to define what you mean by reaching steady state.

• As we should’ve clarified earlier, the “steady state” for our simulations is defined as the point when exchange flow fluxes converge to near-constant values (adding to L125). For the heat budget of the control volume between the glacier and the sill, we did check the storage term in the base cases. Near the end of the 60-day simulation, the heat storage change term is close to zero (Fig. R1, see attachment), and the storage term is small compared to the advective flux term (Fig. R2).

Also, it would be helpful to be more explicit about what actually sets the melt rate in the model. Heat transport (volume flux x temperature) in the lower layer is not a good metric of heat going to melt. Instead, in this model, submarine melting is calculated from a plume parameterization where melt is going to be a function of 3 basic inputs: the subglacial discharge flux, the near-glacier temperature, and near-glacier stratification N2 (and technically salinity, but let’s ignore that for now). First, it would be helpful to say that based on the model setup, the melt will vary with these three parameters. And we know how the modeled melt rate will vary with these three parameters, based on tons of BPT studies. The novel question addressed in this paper is how do sills and reflux modify the near-glacier T and N2. In this model, any change in the circulation regime or heat transport in a layer can’t affect the melt rate unless it changes the near-glacier T or N. For example, consider a hypothetical scenario where an increase in subglacial discharge enhances the exchange flow but does not change temperature and stratification near the glacier. Then you would observe an increased magnitude of heat transport in each layer (volume flux*temperature); however, the modeled melt rate would be unchanged, i.e. the net heat going to glacier melt would be the same.

• To incorporate these useful suggestions we have added new text to provide an explicit dependence of the modeled submarine melting on the three parameters in the Methodology section. The new text reads: 
  
  "In grid locations where subglacial discharge is specified, the submarine melt rate is calculated based on the plume temperature, salinity, and velocity, as well as the ice-ocean boundary layer temperature and salinity (Holland and Jenkins, 1999). In the grid points along the remainder of the glacier front, the melt rate is obtained using the temperature, salinity, and velocity as determined by the MITgcm."

One suggestion – if you don’t want this paper to get bogged down in the details of a full heat budget– would be to keep the focus more on the results about how reflux/sills/etc. affects the temperature and stratification near the glacier. These are the two concrete quantities that then go into the melt parameterization to affect the melt rate in a way that we understand. The changes in the Q volume fluxes are also meaningful, as a well-defined quantity to see the impact of sill and reflux. But multiplying together the volume flux and the temperature within a layer gives a quantity of unclear meaning – unless you evaluate the whole heat budget for the control volume. And then, even if you do go this more complicated route, in the end the heat going to melt is determined by the near-glacier temperature and stratification (given the model’s melt parameterization), so maybe just focus on reporting those simpler metrics. Related, I would emphasize more that you have an expression to predict the lower layer cooling based on the sill. In the text, this point is somewhat under-emphasized because it jumps right to showing the modeled vs. estimated quantity, which is just a proof that relationship holds, but it doesn’t present the main result of how T changes with the controlling parameters (e.g. a plot of T or some metric of cooling vs sill depth).

• We have removed the heat flux calculations, including Table 2, as suggested above.
  
  We have also added a new figure (Fig. R3) that shows how the stratification and temperature at depth change as a result of the presence of the shallow sill, which we hope highlights the importance of understanding the changes in stratification. This complements the other figures that show the change in reflux (Qr) as a function of sill to fjord depth (hs/hf). 

Around L319, it says that says that reflux “result[s] in less heat supply to the glacier but increased submarine melting”. I think this should be “lower heat content near the glacier.” Heat supply that actually goes to the glacier *is* exactly proportional to the submarine melting. I think this gets at the point of confusion between heat transport versus changes in heat content (i.e. the heat storage term in a heat budget).

• Thanks for catching this, we’ve changed the text around L319.

The conclusion implies that regime I and II are steady because heat lost to melting is replaced, whereas regimes III and IV are unsteady because heat lost to melting. Are you confident that regimes III and IV are cooling because of melting? That does not seem supported by the results presented. I would guess that the fjord is cooling because more and more subglacial discharge, which is colder than shelf water, is being mixed into the deep layer. It does not seem to be about heat lost to melting being replaced or not, right?

• This reading of the text was not what we meant and we agree with the reviewer that the cooling in those regimes is largely caused by the cooling driven by subglacial discharge melting with deep warm ocean water. This was correctly stated in the first paragraph of section 2.1, where we mention that the cooling driven by the melting flux is small. The corresponding text is modified to clarify this point. Thanks for catching this.

Finally, two minor questions about the melt implementation in the model: First, is there a heat sink in the way submarine melt is put into the MITgcm? On L44, it says submarine melting is parameterized as a virtual salt flux… is there no heat flux at the terminus associated with submarine melting? Second, at L70-75, the text is a little muddled about submarine melting behind discharge plume versus across rest of the terminus. Are you calculating melt rates across rest of the terminus based on plume velocities (as indicated) or the MITgcm-resolved velocity field or neither? Presumably the melt rates are calculated in the discharge plume using the plume velocity, but what is done for the melt rates across the rest of the terminus?

• Yes, the impact of melting is taken into account by generating a virtual salinity and heat flux that freshens and cools adjacent grid cells. The plume module is coupled to remove the heat from the ambient water in the MITgcm domain, the corresponding text has been modified to clarify this point. Along the remainder of the ice front where there is no runoff input, the melt rate in these cells is calculated using the same formulation as in the BPT model using the ambient temperature, salinity, and velocity. We have added a sentence (included in previous comments) to make this clear. 

2) Volume flux calculations in TEF and reflux

First, I found the notation used for the TEF and the reflux formalism a bit confusing. In the paragraph introducing TEF, Q_in and Q_out are defined as landward and seaward transports. But then for the notation on reflux, with Q^in_1, Q^out_1 the “in” and “out” refer to in and out of the reflux control volume, not landward or seaward. Also, in this latter case, you can’t know whether in/out is landward/seaward or upper/lower layer without knowing where the sections are relative to the control volume. I found this all a bit confusing, especially given the common practice of having “in” and “out” refer to landward and seaward in an estuary. Also, is there a reason to define yet set of layer notation with U_upper and U_lower, etc on L217? (More minor, but it’s also slightly confusing how sometimes the 1/2/3 section indices are subscripts and sometimes superscripts.)

• We’re in agreement with the comment and have made the necessary revisions to clarify and simplify this notation. We are using in/out for transports towards/away from the glacier. The corresponding notions are modified in Fig. 2, Fig. 5, and the text in Section 2.2.

Second, it seems like the beauty of the reflux calculation is that it “expresses vertical fluxes as volume transports, which is equivalent to the horizontal fluxes in TEF”. But then I’m slightly confused about the fact that the layers seemed to be defined differently in TEF vs. the reflux calculations. For TEF, the layers are determined with the ‘dividing salinity’ method, whereas for reflux calculation the layer interface is determined by zero-crossing in the velocity profile. Is this consistent? Can you explain?

• We use the zero-crossing of the layer interface for two purposes only: to calculate the Froude number, and to calculate the volume-averaged temperature (Equation 4, following MacCready et al 2021). This is because a choice has to be made to define the upper/lower layer volumes throughout the entire control volume, i.e., not only at the sections where the TEF/reflux calculations are being made. In MacCready et al 2021, for example, they chose to assume that a fixed percentage of the total volume is the upper layer and provided some sensitivity estimates to show the results didn’t depend on the details of this choice. Here we decided the zero crossing was a reasonable way to look at the two-layer cases without making a priori assumption. The TEF/reflux calculations are otherwise being made using the dividing salinity method at the sections. We’ve clarified this text. 
  
  
  
  References:
  
  MacCready, P., McCabe, R.M., Siedlecki, S.A., Lorenz, M., Giddings, S.N., Bos, J., Albertson, S., Banas, N.S. and Garnier, S., 2021. Estuarine circulation, mixing, and residence times in the Salish Sea. Journal of Geophysical Research: Oceans, 126(2), p.e2020JC016738.

Also it’s a little unclear to me how the material in the Results section around L228-233 about calculating TEF transports relates to the methods outlined in the Methods section. Is this just repetitive or is something different being said? Should it be here or in Methods?

• Here we defined the calculation of entrainment flux Qr between two cross-sections based on TEF transports, which was not mentioned previously and we’ve moved it to the Methodology section.

For both the TEF and reflux calculation, it seems that a 2-layer system is assumed, but how do these calculations work with the regimes that have a 3 layer circulation pattern?

• Thank you for raising an important point. There are several cases that have a 3-layer circulation when the subglacial discharge is small (<50 m³/s), or the initial stratification is set to be strong (>3N₀²). In these cases, a very weak inflow was observed at the fjord surface, which complicates the estimate of TEF transports and refluxes. For this study, the calculations of TEF/reflux were only applied to two-layer cases. We have clarified this in the revised manuscript. 

In Table 3, the difference in upper and lower layer volume flux is exactly 250 m3/s, the subglacial discharge input, for all columns. This is required by mass conservation, right? Thus maybe you only need to show the volume flux in one layer (can state in caption or text that the other layer is the same flux, plus or minus the FW input). Also, you could consider making this into a plot of volume flux vs. h_s/h_f – might be easier to see the trends.

• Yes, the reviewer is correct. We now rely on the figures showing Qr as a function of hs/hf to show the impact of the sill on circulation.

Finally, I found it somewhat confusing how there’s a discussion in the text of Table 3, about how the sill affects volume flux in terms of Q_0^s and Q_0_f. And then later, there’s another discussion about how the sill affects the volume fluxes of Q^2_in and Q^1_out. Could these be consolidated?

• To clarify, the discussion of Table 3 is based on volume budget analysis of the control volume starting from the fjord head to the sill. It serves as an initial investigation of the water transport along the fjord before applying TEF, and helps identify that the strong downward transport occurs over the sill. We’ve moved away from the discussion of these fluxes to focus on the impacts of deep temperature and stratification on the melting, as suggested by the reviewer.

3) Testing the scaling for terminal depth of plume

In this section where you test the scaling for the terminal depth of the plume, is this just testing the plume parameterization at the glacier boundary of the MITgcm domain? My understanding is that the scaling from Slater et al (2016) is derived from BPT, and then here you’re testing if it works well for the MITgcm plume module, which employs BPT to represent the plume. This seems a bit circular, unless I’m missing something. I understand why it’s helpful to introduce the scaling, but it seems like could just be a sentence that says, “based on a BPT study, terminal height of plume scales with (N^2)^(-3/8), and our simulations follow this scaling since the plume module uses BPT to calculate the depth of injection into the fjord domain…” or something. The interesting point is that the sill/reflux changes N^2, and so that changes the terminal height of the plume, right?

In any event, I would try to make clear if this part is testing anything about what is modeled within the MITgcm domain, or is this just testing the BPT parameterization of the plume, in which case there doesn’t seem to be anything novel in testing a scaling that was derived from BPT.

• We agree with this. This section was just meant to show that this existing scaling (which is helpful to define the regimes) works for our modeling output, but we did not mean to imply it was an original result or unexpected. We have removed the figure in question and simplified the text to make this clearer. 

Minor comments

a) How many grid points in the vertical direction are over the sill? For the shallower sills, is this sufficient to resolve sill-related mixing and reflux dynamics?

• In the vertical, the model has 90 grids with resolutions distributed as 2m*20, 4m*30, and 6m*40. For the shallowest sill (hs=16m), there are 8 grid points above. And we tested with a higher resolution (0.5m), the sill-related dynamics and the corresponding results remained consistent.

b) Abstract and L250: I don’t think a % is the best way to quantify cooling. 10% cooling with degrees Celsius would be a very different % on Fahrenheit or Kelvin temperature scale. Would just drop the percentage and say in degrees, or I guess it could be a percentage reduction in the temperature above freezing.

• Yes, we agree with the comment. The text is changed to “the deep fjord temperature dropped by 1℃ ”.

c) When discussing the result that most of the vertical transport occurs over the sill (L206), isn’t that is somewhat by design since there is no topography anywhere else in the fjord? Realistic bathymetry with bends and shallow sides would promote some enhanced mixing. In this case, the sill is really the only place that the flow interacts with topography since it’s a rectangular fjord. This seems like a relevant caveat to mention when stating this result.

• We agree with this. In L207 we mentioned that with our simplified fjord bathymetry mixing processes are expected over the sill and near the glacier. In realistic fjord systems, however, enhanced mixing could also be promoted by bends and shallow regions.

d) In Regime III, when the outflow hits the sill and feeds into the deep layer which is then drawn back towards the glacier, why is this not reflux? Seems like the outflowing layer feeds the deep inflow layer, right? L323 states that “reflux barely generated when plume-driven outflow progressive fills…” – but why isn’t it reflux?

• This is an interesting point and while we were using ‘reflux’ to that occurring near the sill, we agree this broader use is useful and we have modified the text accordingly. 

e) I found the section on tides a bit confusing. One place is says melt increases by 6-30% with tides, then a few sentences later says melt “stayed unchanged with higher tidal amplitudes”. Also, I don’t quite follow the explanation of why reflux decreases with introduction of tides – could you try explain this a bit more clearly?

• At ~L345, we intended to say “within a certain range (Ut/Ue=0.7~1.9), increasing tidal amplitudes has little impact on the submarine melting”.
  
  We understand the importance of providing a clear rationale for the tidally-modified reflux coefficients. An explanation is proposed here: the reflux volume is determined by the outflow flux and the reflux coefficient. The outflow flux increases with stronger tidal forcing because it enhances the exchange flow; the downward reflux fraction was modeled to decrease with stronger tides, a result consistent with the simulation by Hager et al., (2022) (we also mentioned in the discussion Section 4.1). A preliminary understanding is that tides enhance vertical mixing in both directions, suggesting the consequent upward transport may counter the effect of downward reflux. We found the upward reflux fraction increased with stronger tidal forcing (Fig. R4), but further investigation is needed to understand this more generally.
  
  
  
  References:
  
  Hager, A.O., Sutherland, D.A., Amundson, J.M., Jackson, R.H., Kienholz, C., Motyka, R.J. and Nash, J.D., 2022. Subglacial discharge reflux and buoyancy forcing drive seasonality in a silled glacial fjord. Journal of Geophysical Research: Oceans, 127(5), p.e2021JC018355.

Line by Line

L19: citation for first half of sentence?

• We’ve changed the sentence to “Increased submarine melting of glaciers terminating in fjords can be a significant contributor to glacier retreat”.

L25, missing comma: “freshwater, leaving”

• Fixed.

L32: typo? “no deep or no sill” 

• Fixed.

L47: “besides” – word choice 

• We changed it to “In addition”.

L94 shouldn’t “h_s” be “h_f”?

• We’ve changed the text to “... a shallowest depth of hs”.

L95: what do you mean the “fjord domain is set to 2 km”? what dimension is that? The previous sentence says the width is 4 km, and the length seems to be 20 km in Fig 1.

• “The fjord width is set to 2 km”

L96: typo, “this *does* not prevent” 

• Fixed.

L97: says h_f is 400 m in all simulations, but doesn’t Table 1 say that h_f is 200 m in one set? Or what is h (second column) in Table 1?

• Right, hf is 400 m in most cases with an option of 200m for comparison. It should be “hf” in the second column of Table 1. We have corrected the text. 

L110: missing units on Coriolis parameter 

• Fixed.

L113: here and elsewhere, the “~” symbol should be replaced with a dash 

• Fixed.

Table 1: shouldn’t “h” be “h_f”?

• The table has been removed.

L135: missing word: of “of *whether* they are…”

• Fixed.

L179 typo “with with” 

• Fixed.

Figure 3: label glacier location on these (optional, just a suggestion) 

• An arrow was added to indicate the location of the glacier.

L 181: “As the sill depth increases” ambiguous wording, would say as the sill becomes shallower.

• We agree to change the wording.

Table 2 & Table 3: specify for which set of experiments? Looks like this is the first set from Table 1 where Qsg is held constant at 250 m3/s, right?

• Yes, they are from the Base Case.

Figure 6: put a box around the upper left corner point – I assume this is a legend, but it looks like a data point 

• We’ve adjusted the plot as suggested.

L241. Cite Figure 5 here?

• Citation added.

L249” cite Figure in this sentence about cooling?

• We believe the reviewer might mean Figure 6 here? We’ve added a cross-reference.

L256-258: reword this sentence and how it related to the one before 

• The new sentence reads: “A strong stratification constrains the plume terminal height and thus reduces the distance from the plume detachment location at the glacier, it also impacts the overall entrainment of warm ambient water and resulting submarine melting.”

L275: typo? how can hp/hf be greater than 1? Wouldn’t hp/hf = 1 be the plume at surface, and hp/hf<1 be subsurface plume… so hp/hf>1 would be plume above the surface?

• The estimated value of hp from theory can be greater than hf, but we’ve changed the text to say that hp/hf=1 when the plume reaches fjord surface.

Figure 11: the U and T profiles sketched on panel c does not match up with the arrows… should the blue outflow arrow be roughly the same height as the peak of positive U velocity and negative T?

• We’ve adjusted the plot as suggested.

L326-326: run-on sentence

• We’ve modified the sentence and included a new discussion of the regime timescale.

Increasing number of typos towards end of the paper.

• Checked and improved.

L386: should this be H_f^0 or H_s_0? Also, this doesn’t make sense to me 

• We’ve removed the discussion of heat fluxes.

L398: “where background melting dominates the freshwater output”… Background melt is shown to be a significant portion of the total submarine melt, but not the total freshwater input (subglacial discharge is still much larger) 

• Thank you for pointing it out, we’ve modified the wording as “... dominates the total submarine melting”.

Label/refer to regimes in text, if going to label them in figures

• We’ve added additional references to the regime in the text.