Local forcing mechanisms challenge parameterizations of ocean thermal forcing for Greenland tidewater glaciers

Frontal ablation has caused 32 %–66 % of Greenland Ice Sheet mass loss since 1972, and despite its importance in driving terminus change, ocean thermal forcing remains crudely incorporated into large-scale ice sheet models. In Greenland, local fjord-scale processes modify the magnitude of thermal forcing at the ice–ocean boundary but are too small scale to be resolved in current global climate models. For example, simulations used in the Ice Sheet Intercomparison Project for CMIP6 (ISMIP6) to predict future ice sheet change rely on the extrapolation of regional ocean water properties into fjords to drive terminus ablation. However, the accuracy of this approach has not previously been tested due to the scarcity of observations in Greenland fjords, as well as the inability of fjord-scale models to realistically incorporate icebergs. By employing the recently developed IceBerg package within the Massachusetts Institute of Technology general circulation model (MITgcm), we here evaluate the ability of ocean thermal forcing parameterizations to predict thermal forcing at tidewater glacier termini. This is accomplished through sensitivity experiments using a set of idealized Greenland fjords, each forced with equivalent ocean boundary conditions but with varying tidal amplitudes, subglacial discharge, iceberg coverage, and bathymetry. Our results indicate that the bathymetric obstruction of external water is the primary control on near-glacier thermal forcing, followed by iceberg submarine melting. Despite identical ocean boundary conditions, we find that the simulated fjord processes can modify grounding line thermal forcing by as much as 3 °C, the magnitude of which is largely controlled by the relative depth of bathymetric sills to the Polar Water–Atlantic Water thermocline. However, using a common adjustment for fjord bathymetry we can still predict grounding line thermal forcing within 0.2 °C in our simulations. Finally, we introduce new parameterizations that additionally account for iceberg-driven cooling that can accurately predict interior fjord thermal forcing profiles both in iceberg-laden simulations and in observations from Kangiata Sullua (Ilulissat Icefjord).

1 Introduction

Mass loss from the Greenland Ice Sheet (GrIS) contributed 10.8 ± 0.9 mm to mean sea level rise from 1992 to 2018 (The IMBIE Team, 2019) and is projected to raise sea level by 90–180 mm by 2100 (Fox-Kemper et al., 2021). This mass loss has, in part, been triggered by the tidewater glacier response to warming ocean temperatures (e.g., Nick et al., 2009; Holland et al., 2008; Murray et al., 2010; Motyka et al., 2011; Straneo and Heimbach, 2013; Wood et al., 2018), with frontal ablation accounting for 32 %–66 % of GrIS mass loss since 1972 (Enderlin et al., 2014; Van den Broeke et al., 2016; Mouginot et al., 2019). In Greenland, fjords are the principal pathways connecting tidewater glacier termini to the coastal ocean, in which local processes relating to sill-driven mixing and silled obstruction of external water (Mortensen et al., 2011, 2013, 2014; Moffat et al., 2018; Jakobsson et al., 2020; Hager et al., 2022), submarine melting of icebergs and glacier termini (Davison et al., 2020; Jackson et al., 2020; Magorrian and Wells, 2016; Moon et al., 2018), and subglacial discharge (Carroll et al., 2015; Jenkins, 2011) modulate the magnitude of ocean forcing at the ice–ocean boundary, often on a seasonal basis (e.g., Moffat et al., 2018; Mortensen et al., 2013; Hager et al., 2022). However, such processes are too small scale (∼ 1 m to ∼ 1 km length scales at hourly to seasonal timescales) to be resolved in global climate models (grid resolutions of ∼ 30–60 km at annual timescales; e.g., Watanabe et al., 2010; Golaz et al., 2019). To date, sea level rise projections have instead relied on poorly validated simplified parameterizations of oceanic boundary conditions in ice sheet models – such as those developed in Morlighem et al. (2019), Jourdain et al. (2020), and Slater et al. (2020) – that are large sources of uncertainty when predicting future mean sea level (Goelzer et al., 2020; Seroussi et al., 2020). This paper focuses on the ocean thermal forcing of GrIS outlet glaciers, yet Antarctic ice sheet models face similar challenges when prescribing ocean boundary conditions beneath ice shelves (e.g., Seroussi et al., 2020; Jourdain et al., 2020; Burgard et al., 2022).

Recent studies have shown that multi-decadal retreat across a population of tidewater glaciers can be reasonably approximated as a linear function of the climate forcing they experience (Cowton et al., 2018; Slater et al., 2019; Fahrner et al., 2021; Black and Joughin, 2022). For many Greenland tidewater glaciers, change in terminus position is specifically thought to be the result of enhanced submarine melting of the terminus and subsequent changes to ice dynamics (e.g., Holland et al., 2008; Murray et al., 2010; Straneo and Heimbach, 2013; Luckman et al., 2015). Taken together, these findings have prompted the development of parameterizations that use submarine melting to drive frontal ablation in ice sheet models. In particular, the Ice Sheet Intercomparison Project for CMIP6 (Coupled Model Intercomparison Project; ISMIP6) (Nowicki et al., 2016, 2020), which produced sea level contribution projections for Greenland in the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC, 2023), relies on two such parameterizations. The first parameterization (called the retreat implementation) is the simplest – being designed to be implementable in all participating ISMIP6 models – and is used to determine changes in the glacier terminus position over a given time (Slater et al., 2019, 2020):

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }L=\mathit{\kappa }\left(Q_{t_{\mathrm{2}}}^{\mathrm{0.4}}{\mathrm{\Theta }}_{t_{\mathrm{2}}}-Q_{t_{\mathrm{1}}}^{\mathrm{0.4}}{\mathrm{\Theta }}_{t_{\mathrm{1}}}\right),\end{array}$$

where ΔL is the retreat/advance distance (km) between times t₁ and t₂, Q is the mean summer subglacial discharge, Θ is the ocean thermal forcing (°C above freezing temperature), and κ is a coefficient tuned to fit the observed terminus positions of almost 200 Greenland tidewater glaciers between 1960–2018 (Slater et al., 2019). Here, submarine melting is assumed to scale proportionally to Q^0.4Θ.

The second parameterization, called the submarine melt implementation, encompasses only submarine melt and leaves the subsequent glacier response (as given by the relationship between ice flux, submarine melt, and calving) to be calculated by the ice sheet model (e.g., Morlighem et al., 2016). Here, submarine melt rate ($\dot{m}$) is (Rignot et al., 2016)

$$\begin{array}{}\text{(2)}& \dot{m}=(\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}h{q}^{\mathrm{0.39}}+\mathrm{0.15}){\mathrm{\Theta }}^{\mathrm{1.18}},\end{array}$$

where h is grounding line depth and q is the annual mean subglacial discharge normalized by the calving front area. In both implementations, ice sheet models need a method for prescribing Θ based on offshore ocean conditions. In ISMIP6, this was done by first taking a spatial average of annual mean ocean conditions within seven large regional zones surrounding Greenland (Slater et al., 2020). For the retreat implementation (Eq. 1), glaciers are forced with a depth-averaged Θ so that all glaciers within a region experience the same thermal forcing. In the submarine melt implementation (Eq. 2), an adjustment is made accounting for fjord bathymetry preventing deep currents from reaching the glacier face (Sect. 2.2; Morlighem et al., 2019). However, neither parameterization explicitly incorporates water transformation between the coast and glacier termini (e.g., Gladish et al., 2015; Straneo et al., 2012), which can vary greatly even between neighboring fjords (Bartholomaus et al., 2016). Furthermore, accelerated mass loss from the Greenland Ice Sheet can be largely attributed to the dynamics of only a small number of individual glaciers (Enderlin et al., 2014; Fahrner et al., 2021), which can dominate uncertainty in regional retreat projections (Goelzer et al., 2020). There is thus an urgent need for improved parameterizations that incorporate local water transformation and that are validated by high-resolution models or extensive observations.

The large-scale and long-term observations necessary to validate such parameterizations are logistically difficult in Greenland, suggesting a modeling approach is warranted. However, until recently, general circulation models lacked the ability to realistically incorporate iceberg melting, which can be the primary freshwater source in Greenland fjords (Enderlin et al., 2016; Moon et al., 2018). Here, we employ the newly developed IceBerg package (Davison et al., 2020) within the Massachusetts Institute of Technology general circulation model (MITgcm) (Marshall et al., 1997) that enables the inclusion of icebergs to test the accuracy of both ISMIP6 thermal forcing parameterizations across a variety of local forcing conditions. We create a suite of idealized model simulations, each forced with different combinations of subglacial discharge, iceberg prevalence, tidal forcing, and sill geometry but all experiencing the same offshore temperature and salinity conditions at the open boundary. In doing so, our objective is to simulate the diverse array of neighboring fjord conditions that can result from the same regional ocean forcing when local factors are accounted for. We then quantify the error of each ISMIP6 thermal forcing parameterization for all model runs and determine the primary contributors to local water transformation and uncertainty within each formulation. Based on our results, we recommend simple improvements to current thermal forcing parameterizations that substantially improve their accuracy. This paper focuses solely on the accuracy of thermal forcing parameterizations, while assessing the validity of Eqs. (1) and (2), which are provided here only for context, is left to other studies.

Figure 1 (a) Model domain with the location of the sponge layer along open boundaries. (b) Enlargement of the brown box in (a) depicting the locations of the entrance sill, total exchange flow transects (dashed lines), glacier face (blue line), subglacial discharge plume center point (purple dot), and the distribution of iceberg concentrations. (c) Initial and open boundary conditions in relation to the depths of each sill. (d) Along-fjord cross-section of (b) depicting the vertical distribution of iceberg keel depths (binned every 25 m) for each iceberg scenario (labeled by the maximum coverage of grid cell surface area, SA), plotted with the depths of each sill.

Table 1 Forcing conditions and primary diagnostics for each simulation. For subglacial discharge, HC indicates a half-conical plume and L denotes a line plume.

2 Methods

2.1 Model setup

MITgcm fjord geometries were typical of Greenland fjords (e.g., Straneo and Cenedese, 2015) and were 800 m deep, 5 km wide, and 60 km long and had a laterally uniform, Gaussian-shaped sill near the mouth of the fjord with a minimum depth (Z_sill) of either 100, 250, or 400 m (hereafter distinguished as S100, S250, and S400 runs; Fig. 1, Hager, 2023). We varied the sill depth because sills are dominant mechanisms for local water transformation within fjords (e.g., Ebbesmeyer and Barnes, 1980; Cokelet and Stewart, 1985; Hager et al., 2022; Bao and Moffat, 2024), while other geometric constraints, such as length and width, influence fjord circulation (e.g., Carroll et al., 2017) but are not expected to greatly alter fjord water properties. We did not include runs without sills because preliminary simulations showed no difference from S400 runs. The vertical resolution was 10 m in the upper 100 m, 20 m between 100–500 m depth, and 50 m below 500 m. The majority of runs had horizontal resolutions (H_res) of 200 m; however, a few high subglacial discharge and iceberg meltwater flux runs were conducted at 500 m resolution to avoid running at an impractical time step (Table 1). An 800 m deep coastal zone was constructed outside the fjord with an additional 30 cells to the west and 20 cells to the north and south. H_res in the coastal zone linearly telescoped to 2 km at the northern, western, and southern open boundaries. A 10-cell sponge layer was imposed at each open boundary to inhibit internal waves from reflecting back into the domain.

All simulations were run in a hydrostatic configuration with a nonlinear free surface and a Coriolis frequency of 1.3752 × 10⁻⁴ s⁻¹. High- and low-resolution simulations were run at time steps of 25–30 and 60 s, respectively. The horizontal viscosity was prescribed using a Smagorinsky scheme (Smagorinsky, 1963) and a Smagorinsky constant of 2.2, while horizontal diffusivities were set to 0 (though numerical diffusion will still exist). We used the K-profile parameterization (KPP) (Large et al., 1994) for vertical mixing, setting the background and maximum viscosity to 5 × 10⁻⁴ and 5 × 10⁻³ s⁻¹, respectively, and background and maximum diffusivities to 0 and 5 × 10⁻⁵ s⁻¹, respectively. Simulations were run until all fjord water below the sill depth had been flushed and water properties stopped evolving (200–1000 d, depending on fjord flushing time). This step was necessary to ensure simulated fjords had time to fully respond to their unique forcing conditions and to remain consistent with the tacit steady-state assumption made in ISMIP6 thermal forcing parameterizations. The output was averaged over the last 10 d of model time to remove the influence of tides or internal waves, and all “near-glacier” model output was averaged over the two closest grid cells to the glacier face. “Grounding line” water properties are an average of the bottommost near-glacier cells. In total, we ran 27 simulations, each with differing combinations of sill depths, tidal velocities, subglacial discharge plume magnitudes and geometries, iceberg concentrations, and iceberg keel depths (Table 1).

Simulations were initialized from temperature and salinity data profiles observed in 2013–2015 outside the Uummannaq Fjord system, West Greenland (Bartholomaus et al., 2016), which shares a similar vertical structure to summer coastal properties around Greenland: a warm, fresh summer surface layer underlain by cold Polar Water and warm, salty Atlantic Water at depth (Fig. 1c; Straneo et al., 2012; Straneo and Cenedese, 2015). The same profiles were used as boundary conditions along the open boundaries (Fig. 1a). M2 frequency tidal velocities (U_T) of 5 × 10⁻³ m s⁻¹ were imposed along the western boundary, creating tidal amplitudes of ∼ 1.5 m within the fjord, typical of tides throughout East and West Greenland (Howard and Padman, 2021). In shallow-silled fjords (S100 and S250 simulations), where significant tidal mixing was plausible (e.g., Hager et al., 2022; Bao and Moffat, 2024), we tested additional high and low tidal forcing scenarios with U_T values of 7 × 10⁻³ and 3 × 10⁻³ m s⁻¹, creating tides of 2.06 and 0.88 m, respectively (Table 1).

Subglacial discharge (Q_sg) and glacier submarine melting (Q_sm) were parameterized with the IcePlume package (Cowton et al., 2015) using a straight glacier face along the eastern boundary (Fig. 1b, d). Summer high-resolution runs were forced with a subglacial discharge of 300 m³ s⁻¹, which is typical of summer values from Kangilliup Sermia (Rink Isbræ) (Bartholomaus et al., 2016; Carroll et al., 2016). Summer low-resolution runs were designed to resemble the largest Greenland glacial fjords and were forced with a subglacial discharge of 1000 m³ s⁻¹, characteristic of glacier runoff entering Sermilik Fjord and Kangiata Sullua (hereafter referring to Ilulissat Icefjord) (Echelmeyer and Harrison, 1990; Gladish et al., 2015; Carroll et al., 2016; Moon et al., 2018). Subglacial discharge plumes are parameterized to have a half-conical geometry in most runs; however, we test the influence of plume geometry (and thus near-terminus subglacial hydrology) by repeating five runs with subglacial discharge spread out across a 1 km wide line plume (Table 1; e.g., Jenkins, 2011), which may be more realistic for some fjord systems (Jackson et al., 2017; Hager et al., 2022; Kajanto et al., 2023). Winter scenarios were reinitialized from the steady state of summer runs with the same tidal, iceberg, and geometric constraints but were forced by a 10 m³ s⁻¹ line plume across the entire glacier width (Table 1) to account for a switch to basal-friction-generated, distributed subglacial drainage in the winter (e.g., Cook et al., 2020). To be consistent with ISMIP6 methodology, we do not account for seasonal differences in offshore waters; thus, our seasonal sensitivity runs only test seasonal variation in subglacial discharge. In all runs, a background velocity of 0.1 m s⁻¹ was applied along the glacier face to facilitate ambient submarine melting.

Icebergs were parameterized using the IceBerg package (Davison et al., 2020), which treats icebergs as stationary barriers to flow and adjusts surrounding fjord water properties according to calculated iceberg meltwater fluxes (Q_berg). Iceberg depths were set using an inverse-power-law size frequency distribution with an exponent of −1.9, similar to those observed in Sermilik Fjord and Kangilleq – the fjord where Kangilliup Sermia terminates (Sulak et al., 2017). Consistent with observations (e.g., Enderlin et al., 2016; Sulak et al., 2017; Moon et al., 2018; Schild et al., 2021), we prescribed a maximum iceberg depth ($Z_{\text{berg}}^{max}$) of 300 and 400 m in high- and low-resolution runs, respectively (Fig. 1d). Icebergs were concentrated at the fjord head, covering a maximum surface area (${\mathrm{SA}}_{\text{berg}}^{max}$) of either 25 % or 75 % within 10 km of the glacier, which linearly decreased to a minimum surface area (${\mathrm{SA}}_{\text{berg}}^{min}$) of 5 % just inside the entrance sill (Fig. 1c). These concentrations approximate those observed at Kangilleq and Sermilik fjords (Sulak et al., 2017). Additional S250 simulations targeting the ice-choked conditions of Kangiata Sullua were conducted at low resolution using an iceberg concentration of 90 % throughout the fjord. Meltwater plumes resulting from iceberg submarine melt were parameterized by imposing a background velocity of 0.06 m s⁻¹ along the iceberg face. All forcing and geometric conditions (except tidal sensitivity experiments) were repeated with and without icebergs. We did not include sea ice formation in our simulations and do not expect the neglect of associated latent heating and brine rejection to significantly affect our results, particularly in deep fjords. However, it is possible sea ice formation could influence thermal forcing in some shallow fjords.

2.2 Testing of ISMIP6 thermal forcing parameterizations

For comparison to our simulations, we calculate the thermal forcing that would have been imposed in ISMIP6 experiments assuming regional water properties equal to our open boundary conditions. In the first thermal forcing parameterization (ISMIP6retreat; Table 2) used with Eq. (1), thermal forcing is determined by

$$\begin{array}{}\text{(3)}& \mathrm{\Theta }\left(z\right)=\mathit{\theta }\left(z\right)-{\mathit{\theta }}_{\mathrm{f}}\left(z\right)=\mathit{\theta }\left(z\right)-\left[{\mathit{\lambda }}_{\mathrm{1}}S\right(z)+{\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}z],\end{array}$$

where θ and S are the profiles of prescribed potential temperature and practical salinity at the open boundaries; θ_f is the local freezing temperature at depth z; and λ₁ = −5.73 × 10⁻² °C psu⁻¹ (practical salinity unit), λ₂ = 8.32 × 10⁻² °C, and λ₃ = 7.61 × 10⁻⁴ °C m⁻¹ (Jenkins, 2011). Profiles of Θ(z) are then depth-averaged between 200–500 m depth (Slater et al., 2020) to provide a singular value, ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{z}}$ = 4.7 °C (Table 2), across all simulations for ISMIP6retreat. This range of depth was chosen to encompass most Greenland tidewater glacier grounding lines (Slater et al., 2019).

Figure 2 Near-glacier (a–c) salinity and (d–f) thermal forcing profiles for all iceberg (light blue) and non-iceberg (gray) runs (a and d are S400 runs, b and e are S250 runs, c and f are S100 runs). Also plotted are the profiles used to calculate ISMIP6retreat and AMretreat (orange), ISMIP6melt and AMmelt (purple), AMberg (blue; Sect. 4.2), and AMconst (red; Sect. 4.2). The orange profile is also equivalent to the boundary conditions. The gray background in all plots illustrates the range across all runs, and the horizontal dashed line depicts the sill depth. Triangles in (a–c) represent the terminal plume depth for line plumes (white) and half-conical plumes (black). The vertical distributions of iceberg freshwater fluxes (Q_berg) and heat fluxes (H_berg) are shown in (a)–(c) and (d)–(f), respectively, to depict the depth of iceberg melt relative to the sill depth and profile variability.

In contrast to ISMIP6retreat, the submarine melt thermal forcing parameterization (ISMIP6melt) accounts for bathymetry preventing external water from entering the fjord below the sill depth (Table 2; Morlighem et al., 2019). This is accomplished by first defining an effective depth as the deepest part of the near-glacier water column in direct contact with the open ocean (here equal to the sill depth). Fjord water properties above the effective depth are directly extrapolated to the glacier terminus, while water properties below the sill depth are made equal to those at the effective depth (e.g., Morlighem et al., 2019). Extrapolated potential temperature and practical salinity are then converted to in situ temperature and absolute salinity before calculating thermal forcing across the glacier face: Θ(z) = T(z)−T_f(z). Here, T and T_f are the in situ temperature and freezing temperature, which together with the absolute salinity are calculated using the non-linear TEOS-10 toolbox (McDougall and Barker, 2011). As in Slater et al. (2020), the final Θ value used in Eq. (2) is taken from the grounding line depth. The ISMIP6melt formulation therefore predicts the same grounding line thermal forcing (${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{\text{gl}}}$) for all runs within each of the S100 (2.9 °C), S250 (4.6 °C), and S400 (5.5 °C) groups (Fig. 2).

Although methods differ, both ISMIP6 parameterizations are designed to at least crudely target thermal forcing at glacier grounding lines for two reasons: (1) buoyant upwelling along the glacier terminus may homogenize the near-glacier water column (e.g., Mankoff et al., 2016), and (2) grounding line thermal forcing may promote undercutting, and thus calving, of the glacier face. To remain consistent with their original formulations, we test each parameterization against our simulations by comparing it to its own intended metric of thermal forcing near the grounding line. ISMIP6retreat is thus compared to ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{z}}$, the modeled near-glacier thermal forcing averaged between 200–500 m depth, while ISMIP6melt is compared to ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{\text{gl}}}$, the modeled near-glacier thermal forcing at the grounding line. Later, a third thermal forcing metric, ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$, the modeled area-mean (vertically and horizontally averaged) near-glacier thermal forcing, is used to test additional parameterizations (Table 2).

We use root mean square errors (RMSEs) to compare like values between parameterizations and model runs. However, the best depth at which to prescribe thermal forcing at a calving face is still a topic of ongoing debate (see Sect. 4.1). We therefore also employ a Willmott skill score (Willmott, 1981) to assess each parameterization's ability to predict thermal forcing throughout the entire water column. The Willmott skill score (SS) is defined as

$$\begin{array}{}\text{(4)}& \mathrm{SS}=\mathrm{1}-{\displaystyle \frac{\frac{\mathrm{1}}{N}{\mathrm{\Sigma }}_{i=\mathrm{1}}^{i=N}(p_i-m_i{)}^{\mathrm{2}}}{\frac{\mathrm{1}}{N}{\mathrm{\Sigma }}_{i=\mathrm{1}}^{i=N}(\mid p_i-\stackrel{\mathrm{‾}}{m}\mid +\mid m_i-\stackrel{\mathrm{‾}}{m}\mid {)}^{\mathrm{2}}}},\end{array}$$

where m_i is the MITgcm or observed near-glacier thermal forcing at z coordinate i, $\stackrel{\mathrm{‾}}{m}$ is the mean of all m_i values, p_i is the corresponding value predicted by each parameterization, and N is the number of z coordinates in a profile. An SS value of 1 indicates perfect agreement, while a value of 0 indicates no agreement between profiles.

Table 2 Descriptions of ISMIP6 thermal forcing parameterizations, new thermal forcing parameterizations presented in this paper (see Sect. 4.2), and thermal forcing metrics extracted from our simulations.

2.3 Quantification of sill-driven mixing

Sill-driven mixing is a primary control on circulation and water properties in shallow-silled glacial fjords in Alaska and Patagonia (Hager et al., 2022; Bao and Moffat, 2024) and may additionally influence near-glacier thermal forcing in Greenland. Following MacCready et al. (2021), Hager et al. (2022), and Bao and Moffat (2024) we quantify the net effect of sill-driven vertical mixing by pairing the estuarine total exchange flow (TEF) framework (MacCready, 2011) with efflux–reflux theory (Cokelet and Stewart, 1985). We use this approach because it provides bulk mixing transports that are easily relatable to other forcing processes. TEF utilizes isohaline coordinates to identify inflowing and outflowing transports that satisfy the Knudsen relations and account for both tidal and subtidal fluxes (Knudsen, 1900; MacCready, 2011; Burchard et al., 2018). We use 1000 salinity classes to bin salt and volume fluxes across each transect and employ the dividing salinity method (MacCready et al., 2018; Lorenz et al., 2019) to calculate inward and outward transports, allowing for the potential for multiple layers of each to exist. Inflowing and outflowing transport-weighted salinities are given by $S_{\text{in},\text{out}}=F_{\text{in},\text{out}}^s/Q_{\text{in},\text{out}}$, where $F_{\text{in},\text{out}}^s$ and Q_in,out are the sums of all inflowing and outflowing salt and volume fluxes, respectively. We treat temperature as a tracer corresponding to each salt class so that the transport-weighted inflowing and outflowing temperatures are calculated as $T_{\text{in},\text{out}}=F_{\text{in},\text{out}}^{\mathrm{t}}/Q_{\text{in},\text{out}}$ (e.g., Lorenz et al., 2019).

Efflux–reflux theory assumes an estuarine system where mixing is concentrated at constrictions (such as sills) separated by deep basins (dubbed reaches) where mixing is minimal (Cokelet and Stewart, 1985). At each mixing zone, some portion of inflowing or outflowing water is vertically mixed and recirculated, or refluxed, into the opposing layer and back into its original reach, while the remainder, the efflux, is transported across the mixing zone to the next reach (Fig. A1). Using mass and volume conservation, the percentage of inflowing or outflowing water that is refluxed or effluxed can be written in terms of TEF variables (MacCready et al., 2021), but for the purposes of this paper, we are primarily concerned with ${\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}$, which represents the percent of the outflowing fjord water that is refluxed at the entrance sill:

$$\begin{array}{}\text{(5)}& {\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}={\displaystyle \frac{Q_{\text{in}}^{\mathrm{g}}}{Q_{\text{out}}^{\mathrm{g}}}}{\displaystyle \frac{S_{\text{in}}^{\mathrm{o}}-S_{\text{in}}^{\mathrm{g}}}{S_{\text{in}}^{\mathrm{o}}-S_{\text{out}}^{\mathrm{g}}}},\end{array}$$

where superscripts o and g denote the TEF transports on the oceanward and glacierward sides of the mixing zone, respectively (Fig. A1). We calculated efflux–reflux budgets between two TEF transects on either side of the entrance sill and avoid prescribing icebergs within the sill region to ensure temperature and salt are conserved across the mixing zone. More information about using TEF with efflux–reflux theory can be found in MacCready et al. (2021), Hager et al. (2022), and Appendix A.

Figure 3 (a) Near-glacier thermal forcing profiles from all iceberg (light blue) and non-iceberg (gray) runs after removing the depth average of each sill group. (b) The three dominant EOF modes with the percentage of variance they contribute, as calculated from the profiles in (a). Gray triangles indicate terminal plume depths of all runs. The teal line represents variance from iceberg melting, the purple line indicates variance stemming from the boundary conditions, and the gray line signifies variance from subglacial discharge plumes.

2.4 Calculation of local heat fluxes

Quantifying the heat fluxes associated with each local forcing mechanism is important when determining the primary causes of local water transformation. The heat flux resulting from the submarine melting of ice (H_melt) is calculated by

$$\begin{array}{}\text{(6)}& H_{\text{melt}}=-{\mathit{\rho }}_{\text{mw}}{Q}_{\text{mw}}[L+c_i({\mathit{\theta }}_{\mathrm{f}}-{\mathit{\theta }}_i\left)\right],\end{array}$$

where ρ_mw is the meltwater density (1000 kg m⁻³), Q_mw is the total meltwater flux as determined from IceBerg or IcePlume, L is the latent heat of fusion, c_i is the heat capacity of ice, θ_f is the potential freezing temperature, and θ_i is the potential temperature of ice. In our experiments, ice is set to its melting temperature so that Eq. (6) collapses to

$$\begin{array}{}\text{(7)}& H_{\text{berg}}=-{\mathit{\rho }}_{\text{mw}}{Q}_{\text{berg}}L\end{array}$$

and

$$\begin{array}{}\text{(8)}& H_{\text{sm}}=-{\mathit{\rho }}_{\text{mw}}{Q}_{\text{sm}}L\end{array}$$

for the iceberg and glacier submarine melt heat fluxes, respectively.

Advective heat fluxes (H_adv) arising from sill-driven reflux and subglacial discharge are given by

$$\begin{array}{}\text{(9)}& H_{\text{adv}}=\mathit{\rho }{c}_{p}Q({\mathit{\theta }}_{\text{adv}}-{\mathit{\theta }}_{\mathrm{r}}),\end{array}$$

where ρ is the water density, c_p is the heat capacity of water, Q is the advective volume flux, θ_adv is the potential temperature of the advected fluid, and θ_r is a reference temperature. To calculate the recirculatory heat flux caused by sill-driven mixing (here called heat reflux), we substitute TEF quantities into Eq. (9) so that

$$\begin{array}{}\text{(10)}& H_{\text{reflux}}=\mathit{\rho }{c}_p{\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}{Q}_{\text{out}}^{\mathrm{g}}\left(T_{\text{out}}^{\mathrm{g}}-T_{\text{in}}^{\mathrm{g}}\right),\end{array}$$

where ${\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}{Q}_{\text{out}}^{\mathrm{g}}$ is the reflux volume; $T_{\text{out}}^{\mathrm{g}}-T_{\text{in}}^{\mathrm{g}}$ is the temperature difference between the outflowing and inflowing layers on the sill's glacierward side; and ρ is the refluxed water density as determined from $T_{\text{out}}^{\mathrm{g}}$, $S_{\text{out}}^{\mathrm{g}}$, and the mid-column water pressure at 400 m depth. Here, H_reflux refers to the heat transfer (positive or negative) from the outflowing layers to inflowing layers on the sill's glacierward side as a result of sill-driven mixing. For the subglacial discharge heat flux, Eq. (9) is specified as

$$\begin{array}{}\text{(11)}& H_{\text{sg}}={\mathit{\rho }}_{\text{sg}}{c}_p{Q}_{\text{sg}}\left(T_{\text{sg}}-T_{\text{in}}^{\mathrm{g}}\right),\end{array}$$

where Q_sg is the total subglacial discharge; ρ_sg is 1000 kg m⁻³; T_sg is 0 °C; and the reference temperature, $T_{\text{in}}^{\mathrm{g}}$, is chosen to be consistent with Eq. (10). In practice, $T_{\text{in}}^{\mathrm{g}}$ works well as a reference temperature in both Eqs. (10) and (11) as inflowing water properties remain largely unaltered between the sill and glacier face.

2.5 Empirical orthogonal functions of near-glacier variability

Empirical orthogonal function (EOF) analysis was conducted on near-glacier thermal forcing profiles to determine the dominant modes of variability between runs. Near-glacier Θ(z) profiles were horizontally averaged across the calving face before removing the mean Θ within each sill group (Fig. 3a). This second step was necessary to account for the dependence of mean fjord temperatures on the sill depth (Fig. 2). EOFs were then calculated on the resultant profiles.

3 Results

3.1 Near-glacier water properties

Despite identical temperature and salinity forcing at the open boundaries, ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$ varied by 2.7 °C across all runs (Table 1). ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{\text{gl}}}$ varied by 2.9 °C between runs, while grounding line salinities differed by 1.4 psu (Fig. 2, Table 1). In S400 runs, near-glacier water properties were largely unmodified from the open boundary conditions, particularly below 200 m, where the influence of subglacial discharge and iceberg melt was negligible (Fig. 2). However, water properties below the sill depth progressively freshened and cooled as the sill depth shoaled, allowing for ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$ to be neatly grouped by depth of sill: 4.2–4.7 °C for S400 runs, 3.6–4.1 °C for S250 runs, and 2.0–2.9 °C for S100 runs (Fig. 2, Table 1).

Water properties below the sill depth were nearly homogeneous in all runs, with only minor variability occurring when iceberg keels extended below the sill depth (see Q_berg and H_berg in Fig. 2) or subglacial discharge plumes reached neutral buoyancy below the sill depth (this most often occurs with line plumes; see black and white triangles in Fig. 2a-c). The greatest variability existed above 175 m, where most iceberg melting occurred and where most summer subglacial discharge plumes reached neutral buoyancy (Fig. 2). Horizontal gradients in water properties were concentrated over the sill region and were negligible within the fjord (Hager, 2023).

Unsurprisingly, iceberg runs were always cooler than their non-iceberg counterparts; however, the difference between these groups was larger for runs with shallower sills (Table 1, Fig. 2). On average, the difference in ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$ between iceberg and non-iceberg runs diminished from 0.7 °C for S100 runs to 0.3 and 0.2 °C for S250 and S400 runs, respectively. Iceberg melt had the greatest impact on water properties in the upper 175 m, contributing to a temperature range of 5.1 °C at the surface, independent of the sill depth (Fig. 2). However, where iceberg keel depth exceeded the sill depth, iceberg melting cooled the entire water column to the grounding line (Fig. 2), indicating some volume of iceberg meltwater was mixed and refluxed at the silled mixing zone. Such cooling is most pronounced in S100 runs, where the difference in grounding line thermal forcing was on average 0.5 °C between iceberg and non-iceberg runs (Fig. 2). ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$ showed no significant difference between tidal sensitivity runs, nor was there an appreciable distinction between runs with line and half-conical plumes and otherwise equivalent forcing. Winter runs were generally cooler than their summer counterparts in the upper water column, with ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$ varying by ≤ 0.3 °C between winter and summer discharge scenarios with otherwise equivalent forcing (Table 1).

After removing the dominant influence of sills (Fig. 3a), EOF analysis indicates the presence or absence of icebergs accounts for 84 % of the remaining near-glacier thermal forcing variability between runs (Fig. 3b). In general, this first EOF mode reflects the same pattern of cooling in the upper 300 m present in all iceberg runs, but its amplitude changes sign for non-iceberg runs, thus imitating the warm surface water that exists when icebergs are absent (Fig. 3). The second EOF mode makes up 8 % of the variance between runs and has a spatial structure identical to the open boundary temperature conditions (Fig. 3b). Therefore, the second mode can be interpreted as the influence of regional ocean temperatures on near-glacier thermal forcing, in large part accounting for the minimal fjord water mass transformation in S400 non-iceberg runs. A physical interpretation of the third EOF mode, contributing 5 % of the variance, is less certain; however, this mode depicts temperature variability coincident with the terminal depths of subglacial discharge plumes (Fig. 3b) that is not easily relatable to any other forcing mechanism. We thus interpret the third EOF mode to represent variable outflowing plume conditions. As reflux primarily affects the water column below the sill depth, which is homogeneous and similar to depth-averaged water properties, variability resulting from reflux was likely removed with the depth average when computing EOFs. It is therefore possible that reflux variability is incorporated into any of the three dominant modes, part of the remaining 3 % of the variance without clear physical corollaries was removed with the mean during EOF computation.

3.2 Internal freshwater sources and heat fluxes

Subglacial discharge and iceberg meltwater had similar contributions (∼ 30 %–65 %) to the total freshwater input into summer iceberg runs, with glacial meltwater contributing less than 4 % (Table 1). In contrast, iceberg melt flux was the dominant freshwater source (∼ 95 % of all freshwater fluxes) in winter iceberg runs. In non-iceberg runs, subglacial discharge made up ≥ 90 % of freshwater input in the summer, while glacier submarine melt was the dominant freshwater source in the winter (53 %–65 % of all freshwater fluxes; Table 1).

Figure 4 (a) Box plots of heat fluxes associated with reflux (H_reflux), subglacial discharge (H_sg), iceberg submarine melting (H_berg), and glacier submarine melting (H_sm). Box plots depict the maximum, minimum, mean, and standard deviation across all runs. Horizontal blue and orange rectangles depict ranges of estimated surface heat fluxes (Sect. 4.1) for winter and summer runs, respectively. Note that summer surface heat fluxes are reflected across the x axis for illustrative purposes but are actually positive and represent surface warming. (b) Reflux fraction (${\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}$) as a function of $Q_{\text{fw}}^{\prime }/Q_{\text{fw}}$, which is the portion of freshwater input released below the sill depth. Stars and circles differentiate between winter and summer runs, respectively. Marker sizes vary by the sill depth. Orange line is Eq. (12). (c–d) Depth-averaged near-glacier potential temperatures below the sill depth (${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\text{bs}}$) as a function of heat reflux for the S250 and S100 runs. Marker shape differentiates between winter (stars), summer high-resolution (circles), and summer low-resolution (triangles) runs. In (b)–(d), light-blue and gray markers represent iceberg and non-iceberg runs, respectively, and white edges depict runs forced with a line plume.

Despite comparable freshwater fluxes in summer runs, the heat flux from iceberg submarine melting surpassed that from subglacial discharge by multiple orders of magnitude, regardless of iceberg concentration or subglacial discharge (Fig. 4a). Iceberg melting removed heat from surrounding waters at rates of −4.9 × 10⁷ to −3.8 × 10⁸ kW. In contrast, subglacial discharge heat fluxes were −2.2 × 10³ to −7.2 × 10⁶ kW and glacier submarine melt heat fluxes were −29 to −150 kW. Heat reflux spanned 5 orders of magnitude from −8.3 × 10³ to −7.4 × 10⁸ kW, at its maximum exceeding the magnitudes of even the largest iceberg heat fluxes, while at its minimum falling near the lower limits of subglacial discharge heat flux.

As opposed to heat fluxes from freshwater sources, which principally cool the upper water column, heat reflux can directly facilitate the cooling of deep water. Our experiments show a pronounced cooling of deep-water temperatures with increasingly negative heat reflux for both S250 and S100 runs, resulting in a decrease of over ∼ 0.3 and ∼ 0.6 °C, respectively (Fig. 4c–d). In general, the runs with the highest heat reflux contained either icebergs, line plumes, or both; however, there is no clear relationship between heat reflux and any specific local forcing processes (Fig. 4a–b). Nevertheless, there is a highly significant (r=0.91, p = 4.9 × 10⁻¹¹) linear relationship between the portion of freshwater input released below the sill depth and the percent of outflowing water refluxed at the entrance sill, ${\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}$ (Fig. 4b). In our experiments, ${\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}$ can be estimated by

$$\begin{array}{}\text{(12)}& {\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}=\mathrm{0.74}{\displaystyle \frac{Q_{\text{fw}}^{\prime }}{Q_{\text{fw}}}}+\mathrm{0.046},\end{array}$$

where $Q_{\text{fw}}=Q_{\text{sg}}+Q_{\text{sm}}+Q_{\text{berg}}$ is the total freshwater input and the prime denotes freshwater entering the domain below the sill depth (subglacial discharge is included in $Q_{\text{fw}}^{\prime }$ if the plume reaches neutral buoyancy below the sill depth).

All S100 runs had significant sill-driven reflux (${\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}\ge$ 37 %). ${\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}$ in S250 runs ranges from 0 %–66 % but is highest in runs with substantial iceberg freshwater flux below the sill depth or in runs with line plumes (Fig. 4b). ${\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}$ was negligible in all summer S400 runs but became significant in winter runs where weak subglacial discharge plumes still intersected the sill at depth (Table 1). Despite an equivalent ${\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}$ between S100 tidal sensitivity runs (∼ 37 %), tidal forcing does have a minor effect on S250 runs and is responsible for a range of 0.2 %–16 % in ${\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}$ across S250 tidal scenarios (Table 1).

Table 3 RMSE and skill score of thermal forcing parameterizations for different groups of model runs. The RMSE of ISMIP6melt is calculated relative to ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{\text{gl}}}$, while the RMSE of ISMIP6retreat is calculated relative to ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{z}}$. The RMSE of all other parameterizations is relative to ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$. Bold values indicate the most accurate parameterization for each group of runs.

^a RMSE is calculated relative to ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{\text{gl}}}$. ^b RMSE is calculated relative to ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{z}}$.

3.3 Thermal forcing parameterizations

Overall, ISMIP6melt accurately predicted ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{\text{gl}}}$ within an RMSE of 0.2 °C, an error that varied by ± 0.08 °C regardless of the sill depth or iceberg prevalence (Table 3). ISMIP6retreat predicted mean near-glacier thermal forcing between 200–500 m depth (${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{z}}$) within an RMSE of 1.35 °C (Table 3). This error was minimal (RMSE = 0.10 °C) for S400 runs, in which fjord water was of similar composition to shelf water but increased substantially with successively shallower sills (RMSE = 0.71 °C for S250 runs and RMSE = 2.17 °C for S100 runs; Table 3). Both profiles used to compute ISMIP6melt and ISMIP6retreat had moderate success at parameterizing the near-glacier thermal forcing profile, with average skill scores of 0.5–0.6 across all runs (Table 3). Skill scores for iceberg runs were relatively poor (SS = 0.44–0.49) compared to non-iceberg runs (SS = 0.57–0.73). Skill scores also steadily decreased with the sill depth from ∼ 0.70 in S400 runs to 0.31–0.46 in S100 runs (Table 3).

4 Discussion

4.1 Uncertainty in thermal forcing parameterizations

Our results demonstrate the wide range of fjord conditions that can result from equivalent regional ocean forcing and emphasize the need for local processes to be incorporated into the coupling of global climate and ice sheet models. We assign importance to each local forcing mechanism based on its impact on near-glacier thermal forcing, which we corroborate with additional evidence where possible. In order of importance, we identify the dominant local forcing mechanisms as the following.

1. Bathymetric obstruction of external water. The 2.7 °C range in ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$ in our runs is strongly dependent on the depth of the entrance sill, which preferentially blocks warm water from entering the fjord when the sill is shallow. The prominent thermocline between 100–400 m depth in our boundary conditions (Fig. 1c) has been observed in fjords throughout Greenland and marks the transition between Polar Water and Atlantic Water (Straneo et al., 2012). In our experiments, external temperatures range from 0.46 °C at 100 m to 3 °C at 400 m depth (Fig. 1c), nearly the exact range in ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$ between S100 runs and S400. Furthermore, we see no overlap in ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$ between the three sill depths. Therefore, we conclude that the depth of the entrance sill in relation to the 100–400 m thermocline plays a first-order role in determining internal fjord temperatures, but we do not expect the sill depth to be a strong control when sills lie below the Polar Water–Atlantic Water thermocline.

2. Presence or absence of icebergs. After adjusting for the silled obstruction of external water, cooling from iceberg meltwater (or lack thereof) is responsible for 84 % of the remaining variability between runs (Fig. 3), as well as a temperature difference of 5.1 °C at the surface. Iceberg-driven cooling primarily occurs in the upper 175 m of the water column; however, where iceberg keels extend below the sill depth, cooling affects the entire water column as a result of sill-driven reflux. Similar magnitudes of iceberg-driven cooling were modeled in Davison et al. (2020), Davison et al. (2022), and Kajanto et al. (2023). Iceberg meltwater fluxes are comparable to subglacial discharge in summer runs and dominate the freshwater budget in the winter, which is in agreement with prior estimates of iceberg meltwater fluxes in Greenland fjords (Enderlin et al., 2016; Moon et al., 2018). However, the additional energy required to melt this volume of water enables icebergs to disproportionately cool the surrounding water column when compared to a similar volume flux of subglacial discharge.

3. Refluxed outflowing water. Heat reflux can rival the heat flux of iceberg melting, leading to a substantial cooling of deep-water temperatures. Where heat reflux is important, the ≲ 0.6 °C cooling occurs throughout the water column from the sill depth to the grounding line, often affecting a much larger portion of the water column than the melting of icebergs. While such an effect is hard to identify with EOF analysis, the decrease in deep-water temperatures with heat reflux (Fig. 4c and d) indicates this process has the potential to significantly impact near-glacier thermal forcing in certain fjords.

4. Subglacial discharge. While subglacial discharge certainly affects near-glacier thermal forcing, its contribution to variability in near-glacier thermal forcing (both inter- and intra-seasonal) is overshadowed by the influence of icebergs. EOF analysis suggests subglacial discharge is responsible for only 5 % of the near-glacier temperature variability. However, subglacial discharge remains a critical driver of submarine melting through its influence on glacial fjord circulation (Carroll et al., 2015; Sciascia et al., 2013; Straneo et al., 2011; Xu et al., 2012), deep-water renewal (Carroll et al., 2017), turbulent upwelling (Slater et al., 2015), near-glacier horizontal circulation (Slater et al., 2018), and enhanced iceberg melt (Kajanto et al., 2023). Our results further indicate the terminal depth of subglacial discharge plumes can directly affect reflux at silled mixing zones (Fig. 4b) and thus influence deep-water temperatures.

5. Surface heat flux. Although we neglect surface heat fluxes in our model, we approximate their magnitude based on previous estimates from Sermilik Fjord. Surface heat fluxes in Sermilik Fjord are thought to be ∼ 80 W m⁻² in the summer and ∼ −100 W m⁻² in the winter (Jackson and Straneo, 2016; Hasholt et al., 2004). Applying these values to the exposed surface area (not covered by icebergs) in our simulations, we estimate total surface heat fluxes in our simulations would be 3.4 × 10⁶ to 2.0 × 10⁷ kW in the summer and −4.3 × 10⁶ to −2.5 × 10⁷ kW in the winter (Fig. 4a). Thus, surface heat fluxes could often exceed those from subglacial discharge but are an order of magnitude less than iceberg melt heat fluxes. However, surface heating only affects the uppermost water column, and thus we do not expect it to substantially affect ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$ in Greenland fjords. Nevertheless, surface heating may significantly affect heat reflux in Alaska, Svalbard, and Patagonia fjords where shallow sills protrude into the surface layer (e.g., Hager et al., 2022; Bao and Moffat, 2024).

6. Glacier submarine melting. While our modeling suggests glacier submarine melting has little effect on near-glacier thermal forcing, we cannot discount it as an important variable. Our model resolution is too coarse to resolve the complexities of the ice–ocean boundary, and recent observations indicate the IcePlume package may substantially underestimate ambient melting (Jackson et al., 2020).

Observational studies support these model-identified drivers of contrasting water properties between nearby fjords. In northwestern Greenland, slight differences in sill geometry allow for warm Atlantic Water to flow unimpeded into Petermann Fjord, while the inflow of Atlantic Water into the inner basin of neighboring Sherard Osborn Fjord is restricted to a cross-sectional area ∼ 7.5× smaller (Jakobsson et al., 2020). When paired with enhanced sill-driven reflux of glacially modified water, this restricted heat inflow is responsible for a 0.2 °C difference in near-glacier water between these otherwise similar fjords (Jakobsson et al., 2020). Furthermore, Kangilliup Sermia, Kangerlussuup Sermia, and Sermeq Kujalleq (hereafter referring to Jakobshavn Isbræ) are all located within the same ISMIP6 region and are thus subject to the equivalent ocean thermal forcing in ISMIP6 projections (except for bathymetric adjustments used for ISMIP6melt) (Slater et al., 2020). Yet, mean fjord temperatures in the upper 100 m differed by up to 2.5 °C in summer 2014, seemingly due to the large iceberg meltwater flux into Kangiata Sullua, where Sermeq Kujalleq terminates (Bartholomaus et al., 2016; Mojica et al., 2021; Kajanto et al., 2023).

In addition to local forcing processes, our results point to a further source of error in the commonly used thermal forcing parameterizations, namely that current parameterizations only define thermal forcing at specific depths (e.g., at the grounding line or between 200–500 m). For ISMIP6retreat, a depth average of regional ocean thermal forcing from 200–500 m was chosen to encompass most Greenland tidewater glacier grounding lines (Slater et al., 2020), yet this depth range also bounds the largest vertical temperature gradient along the Greenland coast. Thus, any sill that exists within or above this range will greatly affect the ability of ISMIP6retreat to accurately predict near-glacier thermal forcing. Indeed, ISMIP6retreat was quite accurate for S400 runs, in which the difference in 200–500 m temperatures between external and internal water was minimal, but its accuracy decreased progressively with each shallower sill as the temperature difference between external and internal water increased to > 1.5 °C.

The current state of frontal ablation laws used in large-scale ice sheet models requires a single scalar to represent thermal forcing across the entire glacier terminus (e.g., Aschwanden et al., 2019; Choi et al., 2021). Unfortunately, this makes it difficult to define the best metric for prescribing near-glacier thermal forcing. Current methods often rely on grounding line conditions, motivated by the fact that grounding line submarine melt may promote undercutting and calving (e.g., Rignot et al., 2016; O'Leary and Christoffersen, 2013); however, glacier termini can actually exhibit a wide range of geometries resulting from spatially variable melt rates (Fried et al., 2019), which can either increase or decrease calving depending on the magnitude and shape of the entire melt profile (Ma and Bassis, 2019). Thus, an overemphasis on undercutting processes may lead to spurious frontal ablation laws. Furthermore, the choice to prioritize grounding line thermal forcing in both ISMIP6melt and ISMIP6retreat negates the iceberg-driven cooling of the upper water column, which our modeling indicates is the secondary control on near-glacier thermal forcing. While buoyant upwelling may partially homogenize near-terminus waters (at a finer scale than our model resolution), particularly within the subglacial discharge plume (e.g., Mankoff et al., 2016), it seems likely that iceberg meltwater remains an important control on thermal forcing at the ice–ocean boundary. Although it remains unclear at what depth glacier dynamics are most sensitive to ocean thermal forcing, from an ocean forcing perspective, it may be prudent to define thermal forcing using a metric that captures all processes relevant to near-glacier conditions, such as iceberg melting. Such an approach would also pave the way for coupled mélange–ice dynamics models, which would require an accurate portrayal of the upper water column to appropriately model mélange dynamics and prescribe back stress to the glacier. Therefore, we here test alternative methods for parameterizing a full profile of ocean thermal forcing, as well as an area mean, which allow for the inclusion of all dominant fjord-scale processes, while remaining relevant to calving processes.

4.2 Alternative thermal forcing parameterizations

We first test revisions to ISMIP6melt and ISMIP6retreat that are independent of specific depths, here called area-mean melt (AMmelt) and area-mean retreat (AMretreat), respectively (Table 2). Thermal forcing profiles are calculated as was done for ISMIP6melt and ISMIP6retreat, but they are extrapolated across the glacier terminus, with an area-mean thermal forcing being calculated instead of pulling from a set depth (Fig. C1). A third parameterization (AMberg) accounts for bathymetric obstructions in the same way as ISMIP6melt, but the surface 175 m follows the Gade slope submarine melt mixing line (Gade, 1979; Straneo and Cenedese, 2015) to approximate the influence of iceberg submarine melting on the upper water column (Fig. 2). For ice at the freezing temperature, the Gade slope is

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}S}}={\displaystyle \frac{\mathrm{1}}{S\left(z\right)}}\left({\displaystyle \frac{L}{c_i}}-\left(T_{\mathrm{f}}\right(z)-T(z)\right).\end{array}$$

As surface salinity in the fjord remains relatively similar to external water (Fig. 2), we leave salinity unchanged from ISMIP6melt and adjust temperatures accordingly to fit the Gade slope (Table 2). A final modification is made in which waters below freezing temperature are set to the in situ freezing point (Fig. C1). In practice, AMberg approximates a lower limit of cooling that occurs through iceberg melting, while AMmelt sets the upper bound on surface temperatures attainable for non-iceberg runs (Fig. 2). Therefore, we test a fourth, middle-ground parameterization (AMconst) in which the ISMIP6melt temperature and salinity at 175 m is extrapolated to the surface before calculating the thermal forcing profile (Fig. 2; Table 2). For both AMconst and AMberg, the final thermal forcing value used in Eqs. (1) and (2) is again an area mean across the glacier face.

The profiles used to calculate AMmelt and AMretreat each had a skill score of ∼ 0.55 and predicted ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$ within an RMSE of 0.30 and 1.24 °C, respectively; however, the RMSE of AMretreat remained heavily dependent on the sill depth (Table 3). As expected, AMmelt was a strong predictor of near-glacier thermal forcing in non-iceberg runs (RMSE = 0.13 °C; SS = 0.73), while AMberg performed well in iceberg runs (RMSE = 0.19 °C; SS = 0.91). AMconst was an adequate compromise between AMberg and AMmelt, with an RMSE of 0.31 °C and a skill score of 0.69 across the entire model ensemble (Table 3). Skill scores for all parameterizations decreased considerably with successively shallower sills (Table 3).

Figure 5 Thermal forcing profiles (gray) observed within Kangiata Sullua (KS) via XCTD casts in (a) August 2014 and (b) August 2019, plotted with the profiles used to calculate ISMIP6retreat (orange) and AMberg (dark blue). The ISMIP6retreat profile is equivalent to external forcing conditions in Disko Bay. Skill score for each parameterization profile is also provided. Horizontal dashed line signifies the maximum depth (245 m) of the Kangiata Sullua entrance sill. (c) Bathymetric map rendering the locations of XCTD casts within Kangiata Sullua (circles) used in (a) and (b), as well as XCTD and CTD casts in Disko Bay (triangles) used to provide external forcing conditions. Solid rose line depicts the Sermeq Kujalleq terminus position in July 2014 (Goliber et al., 2022), and the dashed black line outlines the location of the entrance sill.

Encouraged by our model results, we test the efficacy of AMberg with observations using conductivity, temperature, and depth profiles collected by ship (CTD) and helicopter (XCTD, expandable) within Kangiata Sullua and the adjacent Disko Bay in August 2014 (Beaird et al., 2017) and August 2019 (Fig. 5; Appendix B). Temperature and salinity data from Disko Bay were used to calculate an average external thermal forcing profile analogous to ISMIP6retreat. Following the steps outlined above, a profile of AMberg was then constructed based on average Disko Bay temperature and salinity profiles. In both years, AMberg reproduced XCTD casts within Kangiata Sullua with a mean skill score of 0.95 (Eq. 4), compared to mean ISMIP6retreat and ISMIP6melt skills scores of 0.43–0.62 and 0.19–0.27, respectively. Agreement between AMberg and observations was particularly pronounced in the upper water column where vast volumes of meltwater are released from the fjord's dense and extensive mélange (Enderlin et al., 2014). However, observed thermal forcing below the sill depth was higher than predicted by AMberg and ISMIP6melt (Fig. 5), which is in contrast to our Kangiata Sullua-style simulations (Fig. 2e). We interpret this discrepancy to be the result of temporally varying conditions in Kangiata Sullua that reflect previously warmer conditions in Disko Bay and not due to steady-state along-fjord temperature gradients, which were negligible in our runs (Hager, 2023).

While AMberg demonstrates a marked improvement from other parameterizations in capturing the cooling effects of icebergs, it does so by sacrificing accuracy in non-iceberg runs. Therefore, the true merit of AMberg would be its tandem use with AMmelt and some a priori knowledge or likelihood estimate of iceberg presence, whereby AMmelt could be used in fjords with few icebergs and AMberg could be used where icebergs are prevalent (here called AMfit). Such an approach could predict near-glacier thermal forcing within an error of 0.15 °C and a skill score of 0.82 across all model runs (Table 3) but is particularly effective in S400 (RMSE = 0.10 °C; SS = 0.91) and S250 runs (RMSE = 0.09 °C; SS = 0.83).

4.3 Impact of thermal forcing parameterizations on ISMIP6 submarine melt rates

The substantial range between thermal forcing parameterizations could significantly impact modeled submarine melt rates. However, the submarine melt parameterizations used in MITgcm and ISMIP6 do not currently resemble observed melt rates at tidewater glacier termini (e.g., Jackson et al., 2020; Sutherland et al., 2019), and thus comparison of absolute melt rates may be misleading. Instead, we here compare only the relative range of submarine melt rates provided by Eq. (2) between the various thermal forcing implementations. (Note this requires using Eq. (2) slightly outside its intended purpose, which was to prescribe submarine melting based solely on thermal forcing defined between 200 m depth and the grounding line.)

Figure 6 The range of ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{z}}$, ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$, and ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{\text{gl}}}$ for each simulation (ΔΘ), plotted with the resulting range of parameterized submarine melt rates ($\mathrm{\Delta }\dot{m}$), as defined by Eq. (2). A ΔΘ of 0 indicates perfect agreement between ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{z}}$, ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$, and ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{\text{gl}}}$ for a single run, while a large ΔΘ indicates these thermal forcing metrics vary substantially. Iceberg runs are light blue, and non-iceberg runs are gray. White-edged markers depict runs with line plumes, stars represent winter (Q_sg = 10 m³ s⁻¹) simulations, and triangles indicate low-resolution (Q_sg = 1000 m³ s⁻¹) runs. Marker size varies by the sill depth. Green markers depict the range of all seven thermal forcing parameterizations for each sill depth and subglacial discharge group, again plotted against the resultant range of parameterized submarine melt rates. For green markers, a ΔΘ of 0 indicates perfect agreement between all seven thermal forcing parameterizations, while a large ΔΘ indicates substantial variation between parameterizations.

Our modeling indicates that the magnitude of thermal forcing used in ice sheet models, and thus parameterized submarine melt rates, is highly sensitive to the depth range where thermal forcing is defined. In a given simulation, ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{z}}$, ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$, and ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{\text{gl}}}$ differed by 0.4–1.1 °C (Table 1, Fig. 6), which can translate to a difference of > 200 m yr⁻¹ in submarine melt rates parameterized by Eq. (2) (Fig. 6). An even greater thermal forcing range existed between the seven parameterizations tested in this paper (Table 2). Within a given run, thermal forcing parameterizations ranged from 0.90–2.26 °C (Fig. 6). Depending on the thermal forcing method used, submarine melt rates parameterized by Eq. (2) could therefore differ by 91–617 m yr⁻¹ for a single fjord (Fig. 6). Melt rate parameterization uncertainty of this magnitude (≤ 1.7 m d⁻¹) could significantly impact long-term projections of submarine melt for some Greenland glaciers (compare to projected melt rates in Fig. 10 of Slater et al., 2020). The range of parameterizations was greatest in S100 summer runs, largely due to the pronounced bathymetric effects that ISMIP6retreat and AMretreat do not account for. Also notable is the difference between thermal forcing parameterizations in our Kangiata Sullua-style low-resolution runs (Fig. 6), which due to the high subglacial discharge of Sermeq Kujalleq would contribute an uncertainty in parameterized submarine melt rates of > 300 m yr⁻¹. Sermeq Kujalleq single-handedly accounts for 10 % of GrIS mass loss (Smith et al., 2020) and dominates uncertainty in mass balance projections of central West Greenland (Goelzer et al., 2020). Therefore, this wide range between thermal forcing parameterizations could signify substantial uncertainty in the projected terminus position of Sermeq Kujalleq and thus mass loss projections for the GrIS as a whole.

4.4 Remaining uncertainty associated with sill-driven mixing

While AMberg effectively parameterizes the two largest sources of uncertainty in predicting near-glacier thermal forcing – bathymetric obstruction and iceberg meltwater – it does not account for the modification of deep water through sill-driven reflux, which was found to be significant in this study and in previous work (Hager et al., 2022; Davison et al., 2022; Kajanto et al., 2023; Bao and Moffat, 2024). As discussed in Sect. 3.1, the influence of reflux is difficult to discern through EOF analysis, although multiple lines of evidence highlight its importance in shallow-silled fjords. First, heat reflux has the potential to exceed even the greatest iceberg heat fluxes and is responsible for the cooling of deep water (88 % of the water column in S100 runs) by up to ∼ 0.6 °C (Fig. 4). Second, the RMSE of each parameterization is greatest in S100 runs, in which more than 37 % of outflow water is refluxed. This is true despite correcting for the bathymetric obstruction of external water in most parameterizations. Thus, there must be an additional source of error related to the sill depth, which can only be explained through the cooling of deep water through the reflux of iceberg meltwater and subglacial discharge.

Deep-water cooling from sill-driven mixing is not expected to be important in fjords with sills deeper than iceberg keels or the summer terminal plume depth; however, reflux is likely influential in a number of critical glacial fjord systems. Only 10 individual glaciers are responsible for the majority of GrIS mass loss (Fahrner et al., 2021), 3 of which, Kakiffaat Sermiat, Sverdrup Gletsjer, and Sermeq Kujalleq, terminate in fjords with sills ≲ 250 m deep (as compared to grounding line depths of 400–800 m) (Morlighem et al., 2017). It is in these fjords that we expect heat reflux to significantly influence near-glacier temperatures. Indeed, recent modeling of Kangiata Sullua indicates the sill-driven reflux of iceberg meltwater cools near-glacier water by 0.2 °C (Kajanto et al., 2023), a result shared by our Kangiata Sullua-style low-resolution runs.

Although the updated parameterizations presented in this paper greatly reduce thermal forcing error compared to existing ISMIP6 methods, incorporation of sill-driven mixing could further reduce error in shallow-silled fjords, such as the Kangiata Sullua–Sermeq Kujalleq system. We anticipate such improvements would require the use of box models that contain representations of iceberg melting, subglacial discharge, and reflux. The strong linear dependence of ${\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}$ on $Q_{\text{fw}}^{\prime }/Q_{\text{fw}}$ (Fig. 4b) indicates reflux fractions can be accurately estimated from the vertical distribution of freshwater fluxes in the water column, without any knowledge of tidal forcing. Thus, Eq. (12) could be used within a box model to predict reflux, assuming the model can approximate freshwater fluxes throughout the water column.

5 Conclusions

In summary, we have tested the accuracy of common ocean thermal forcing parameterizations across a wide range of local forcing scenarios and fjord geometries and identified fjord bathymetry and iceberg-melt-driven cooling as the two greatest sources of error when translating regional water properties to tidewater glacier termini. Even a simple adjustment for fjord bathymetry, as done in the ISMIP6 submarine melt implementation, can predict grounding line thermal forcing within an average error of 0.2 °C; however, without accounting for bathymetry, parameterizations overpredict thermal forcing by at least 2 °C in our idealized shallow-silled fjords.

Neither ISMIP6 parameterization could reliably reproduce an entire profile of near-glacier thermal forcing in our simulations due to iceberg-driven cooling of the upper water column. To address this concern, we made a simple correction to the ISMIP6 submarine melt implementation that can predict the near-glacier thermal forcing profiles of our iceberg-laden runs, as well as observed thermal forcing within Kangiata Sullua, to a skill score of > 0.90. While promising, this approach sacrifices accuracy in fjords with few icebergs and thus would be best utilized in conjunction with an iceberg prevalence prediction method, which could apply the iceberg parameterization only to fjords where iceberg melt is significant.

Our modeling also highlights the sensitivity of thermal forcing – and therefore the parameterized submarine melt rate – to the depth range that is chosen. The uncertainty this contributes to thermal forcing in ice sheet models could be substantial, and determining where best to define thermal forcing at a calving face should be a topic of high importance. For example, if a depth range that incorporates surface layers is chosen, a calving front area-mean thermal forcing, then iceberg-driven cooling becomes important and can reduce the prescribed thermal forcing by 1 °C in many Greenland fjords. Regardless, the AMfit parameterization presented here can predict entire near-glacier thermal forcing profiles with an average skill score of 0.82 across all bathymetries and iceberg concentrations tested in our simulations. As such, AMfit could be used to accurately parameterize ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{z}}$, ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{A}}$, ${\mathrm{\Theta }}_{\stackrel{\mathrm{‾}}{\text{gl}}}$, or an entire Θ(z) profile, depending on the requirements and capabilities of the ice sheet model. This result highlights the need for an iceberg prevalence prediction method to be developed and implemented in the next generation of ice sheet models.

Additional improvements to the parameterizations presented here could take the form of a box model that can effectively represent sill-driven reflux. While such an approach is outside the scope of this paper, we have identified a strong linear relationship between reflux and the ratio of freshwater released below the sill depth (Eq. 12), which we suggest could be used to efficiently estimate reflux in box models of Greenland fjords. We emphasize that such models also need to accurately parameterize iceberg melting and subglacial discharge plumes, as well as incorporate local fjord bathymetry. Still, the revised thermal forcing parameterizations presented in this paper are an improvement to existing methods, while their simplicity makes them relatively straightforward to implement. Limited observational evidence throughout Greenland supports the efficacy of these thermal forcing parameterizations, yet robust validation is still needed in realistic settings.

Appendix A: Total exchange flow (TEF) and efflux–reflux theory

Efflux–reflux theory quantifies the net effect of mixing without the need to resolve individual mixing processes (Cokelet and Stewart, 1985). Effectively, efflux–reflux transports can be thought of as the vertical equivalent of the horizontal TEF budget (MacCready et al., 2021). In TEF terms, mass and volume conservation require

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}{\mathit{\alpha }}_{\text{in}}^{\mathrm{r}}{Q}_{\text{in}}^{\mathrm{o}}+{\mathit{\alpha }}_{\text{out}}^{\mathrm{e}}{Q}_{\text{out}}^{\mathrm{g}}& =Q_{\text{out}}^{\mathrm{o}},\\ {\mathit{\alpha }}_{\text{in}}^{\mathrm{r}}{S}_{\text{in}}^{\mathrm{o}}{Q}_{\text{in}}^{\mathrm{o}}+{\mathit{\alpha }}_{\text{out}}^{\mathrm{e}}{S}_{\text{out}}^{\mathrm{g}}{Q}_{\text{out}}^{\mathrm{g}}& =S_{\text{out}}^{\mathrm{o}}{Q}_{\text{out}}^{\mathrm{o}},\\ {\mathit{\alpha }}_{\text{in}}^{\mathrm{e}}{Q}_{\text{in}}^{\mathrm{o}}+{\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}{Q}_{\text{out}}^{\mathrm{g}}& =Q_{\text{in}}^{\mathrm{g}},\\ {\mathit{\alpha }}_{\text{in}}^{\mathrm{e}}{S}_{\text{in}}^{\mathrm{o}}{Q}_{\text{in}}^{\mathrm{o}}+{\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}{S}_{\text{out}}^{\mathrm{g}}{Q}_{\text{out}}^{\mathrm{g}}& =S_{\text{in}}^{\mathrm{g}}{Q}_{\text{in}}^{\mathrm{g}},\end{array}\end{array}$$

where superscripts o and g denote the TEF transports on the oceanward and glacierward sides of the mixing zone, respectively, and subscripts indicate whether the transport is inflowing (glacierward) or outflowing (oceanward). Superscripts on α signify the percent of the inflowing or outflowing layer that is refluxed (${\mathit{\alpha }}_{\text{in},\text{out}}^{\mathrm{r}}$) or effluxed (${\mathit{\alpha }}_{\text{in},\text{out}}^{\mathrm{e}}$) at the sill (Fig. A1). The solutions to Eq. (A1) are

$$\begin{array}{}\text{(A2)}& {\displaystyle }\begin{array}{rl}{\mathit{\alpha }}_{\text{in}}^{\mathrm{r}}={\displaystyle \frac{Q_{\text{out}}^{\mathrm{o}}}{Q_{\text{in}}^{\mathrm{o}}}}{\displaystyle \frac{S_{\text{out}}^{\mathrm{o}}-S_{\text{out}}^{\mathrm{g}}}{S_{\text{in}}^{\mathrm{o}}-S_{\text{out}}^{\mathrm{g}}}},& {\mathit{\alpha }}_{\text{out}}^{\mathrm{e}}={\displaystyle \frac{Q_{\text{out}}^{\mathrm{o}}}{Q_{\text{out}}^{\mathrm{g}}}}{\displaystyle \frac{S_{\text{in}}^{\mathrm{o}}-S_{\text{out}}^{\mathrm{o}}}{S_{\text{in}}^{\mathrm{o}}-S_{\text{out}}^{\mathrm{g}}}},\\ {\mathit{\alpha }}_{\text{in}}^{\mathrm{e}}={\displaystyle \frac{Q_{\text{out}}^{\mathrm{g}}}{Q_{\text{out}}^{\mathrm{o}}}}{\displaystyle \frac{S_{\text{in}}^{\mathrm{g}}-S_{\text{out}}^{\mathrm{g}}}{S_{\text{in}}^{\mathrm{o}}-S_{\text{out}}^{\mathrm{g}}}},& {\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}={\displaystyle \frac{Q_{\text{in}}^{\mathrm{g}}}{Q_{\text{out}}^{\mathrm{g}}}}{\displaystyle \frac{S_{\text{in}}^{\mathrm{o}}-S_{\text{in}}^{\mathrm{g}}}{S_{\text{in}}^{\mathrm{o}}-S_{\text{out}}^{\mathrm{g}}}}.\end{array}\end{array}$$

Mass and volume conservation also require

$$\begin{array}{}\text{(A3)}& \begin{array}{r}{\mathit{\alpha }}_{\text{in}}^{\mathrm{r}}+{\mathit{\alpha }}_{\text{in}}^{\mathrm{e}}=\mathrm{1},\\ {\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}+{\mathit{\alpha }}_{\text{out}}^{\mathrm{e}}=\mathrm{1}\end{array}\end{array}$$

and

$$\begin{array}{}\text{(A4)}& \begin{array}{r}{S}_{\text{out}}^{\mathrm{g}}\le S_{\text{out}}^{\mathrm{o}}<S_{\text{in}}^{\mathrm{o}},\\ S_{\text{out}}^{\mathrm{g}}<S_{\text{in}}^{\mathrm{g}}\le S_{\text{in}}^{\mathrm{o}}.\end{array}\end{array}$$

TEF budgets are not exact, and even at a steady state some drift still occurs within the mixing zone; therefore, we make minor adjustments to the TEF transports, ensuring Eqs. (A1), (A3), and (A4) are satisfied before solving Eq. (A2) (e.g., MacCready et al., 2021; Hager et al., 2022), but the resultant error in ${\mathit{\alpha }}_{\text{out}}^{\mathrm{r}}$ was ≤ 0.04 % of the reported value.

Figure A1 Schematic illustrating each variable in efflux–reflux theory (Eq. A1).

Appendix B: Collection and processing of hydrographic data

Hydrographic properties were observed via helicopter-based XCTD casts within Kangiata Sullua (four profiles in August 2014 and five in August 2019) and by shipboard CTD casts in Disko Bay (four profiles in August 2014 and five in August 2019). In both years, an additional XCTD profile was obtained in Disko Bay for comparison with CTD casts. Shipboard CTD casts were conducted using an RBR XR-620 in 2014 and an RBRconcerto in 2019. All casts were averaged into 1 m bins and smoothed using a running mean filter with a 10 m window (comparable to our model vertical resolution) before undertaking the analysis described in Sect. 4.2. Further details about the 2014 field campaign are described in Beaird et al. (2017).

Appendix C: Parameterization flowchart

Figure C1 Flowchart illustrating the step-by-step process for computing AMberg and AMmelt.