FjordRPM v1.0: a reduced-physics model for efficient simulation of glacial fjords

Interactions between ice masses and the ocean are key couplings in the global climate system. In many cases, these interactions occur through glacial fjords, which are long, deep, and narrow troughs connecting the open ocean to marine-terminating glaciers. By controlling the fluxes of ocean heat towards the ice sheet and ice sheet freshwater towards the ocean, glacial fjords play an important role in modulating ice sheet mass loss and the impacts of freshwater on ocean circulation. Yet, these dynamics occur at small scales that are challenging to resolve in earth system models and hence are often ignored, represented in an ad-hoc manner, or studied using expensive high-resolution models that are limited in scope. Here, we propose a means of capturing glacial fjord dynamics at negligible computational expense in the form of a “reduced-physics” model (FjordRPM) that resembles a “1.5-dimensional” or box model. We describe the design and physical parameterisations in the model and demonstrate its ability to capture important modes of glacial fjord circulation by comparing it against a general circulation model in idealised and realistic simulations. We suggest that the model is a useful tool for understanding fjord dynamics and a promising approach for representing glacial fjord processes within large-scale models or climate and sea level projection efforts.

1 Introduction

There is increasing recognition that ice–ocean interaction is a key exchange in our earth system. Ocean heat flux towards ice sheet margins is linked to ice mass loss and sea level contribution (Straneo and Heimbach, 2013), ice sheet freshwater has the potential to influence large-scale ocean dynamics (Böning et al., 2016), and ice–ocean processes fertilise ecosystems (Oliver et al., 2023). One prominent system in which these interactions occur is glacial fjords, which link marine-terminating glaciers to the open ocean. Glacial fjords are found in glaciated regions around the world such as the West Antarctic Peninsula, Svalbard, and Alaska and are particularly numerous in Greenland. Around two-thirds of the mass loss from the Greenland Ice Sheet since the 1970s has occurred at marine-terminating glaciers that flow into glacial fjords (Mouginot et al., 2019), and there is concern that the associated increased freshwater flux to the ocean may be capable of weakening the Atlantic Meridional Overturning Circulation (Frajka-Williams et al., 2016; Jackson et al., 2023; van Westen et al., 2024). As such, the ability to capture the effects of glacial fjord dynamics in models used for climate and ice sheet projections is crucial.

Glacial fjords, with a typical width of 2–10 km, are, however, much too small to resolve in earth system or global ocean models and generally too small to include in even regional ocean models (e.g. Zuo et al., 2019; Nguyen et al., 2021). When considering the impact of the ocean on the Greenland Ice Sheet in large-scale simulations (e.g. Goelzer et al., 2020), glacial fjord dynamics have, therefore, either been ignored or represented in a very basic manner (Slater et al., 2020), missing important details of how subglacial discharge, icebergs, shelf winds, and sills all modify the properties of ocean waters reaching marine-terminating glaciers (Mortensen et al., 2014; Jackson et al., 2014; Straneo and Cenedese, 2015; Hager et al., 2022; Davison et al., 2020). Similarly, when considering the impact of Greenland Ice Sheet freshwater on ocean dynamics, models often impose the freshwater as an unmixed and distributed surface flux (e.g. Golledge et al., 2019) or make an a priori choice of depth distribution (e.g. Gillard et al., 2016), neither of which is faithful to our understanding of how ice sheet liquid freshwater enters the ocean as a highly diluted, subsurface flux that can be delayed in time (Beaird et al., 2018; Sanchez et al., 2023) and of how a significant portion of the solid ice flux melts inside fjords (Moon et al., 2018; Moyer et al., 2019). In sum, better representation of glacial fjord processes within large-scale studies is needed to fully understand ice–ocean interactions and their impact on ice masses and climate.

Here, we introduce a means of approximating glacial fjord dynamics at negligible computational cost, which we name FjordRPM (“Fjord Reduced-Physics Model”). The model describes the volume-mean fjord temperature and salinity in a number of discrete, fixed, vertically stacked layers. Exchange processes between the fjord, glacier, shelf, and icebergs are parameterised using known or adapted expressions. The resulting model lies somewhere between a box model and a “1.5-dimensional” model (i.e. a one-dimensional model with lateral fluxes), and hence we adopt the terminology “reduced-physics model”. Models of this type have a history of providing insight into ocean systems (e.g. Babson et al., 2006; Gillibrand et al., 2013) and of representing coastal processes within larger-scale ocean models (e.g. Sun et al., 2017). It is our hope that FjordRPM can be of similar use within the field of ice sheet–ocean interaction.

The paper proceeds as follows. We first describe the formulation of FjordRPM and the parameterisation of exchanges (Sect. 2), followed by the method of solution, the structure of the code, and the required inputs (Sect. 3). We demonstrate that the model can capture fundamental modes of glacial fjord circulation by comparing to simulations conducted using the established general circulation model MITgcm (Sect. 4). We conclude by reflecting on the strengths and weaknesses of the model and proposing avenues for future work (Sect. 5).

2 Model description

2.1 Preliminaries

Glacial fjords are host to a range of processes that serve to influence a fjord's circulation and properties (e.g. Straneo and Cenedese, 2015). In this first iteration of FjordRPM, we focus on a small number of key processes that we deem to be of primary importance (Fig. 1). Two of these are specific to glacial environments and are recognised to be of particular significance in glacial fjords. Firstly, the input of subglacial discharge at the glacier grounding line drives buoyant plumes that mix with and upwell ambient waters before reaching the fjord surface or finding neutral buoyancy at an intermediate depth (Jenkins, 2011; Jackson et al., 2017; Hewitt, 2020). Secondly, the melting of icebergs calved into the fjord from tidewater glaciers acts as a heat sink and freshwater source while driving convective upwelling adjacent to icebergs (Moon et al., 2018; Davison et al., 2020; Cenedese and Straneo, 2023). In addition to these processes, we include the exchange of water between the fjord and continental shelf, as driven by pressure gradients in response to spatial and temporal variation in water properties (e.g. Jackson et al., 2014), and vertical mixing and advection between model layers.

Figure 1 Schematic of FjordRPM. The model layers (j) have the static properties of volume, V_j, and iceberg surface area, I_j, and the dynamic properties of temperature, T_j(t), and salinity, S_j(t). The layers have horizontal exchanges of volume (X=V), heat (X=T), and salt (X=S) with the plume ($Q_j^{Xp}$) and the shelf ($Q_j^{Xs}$). They have vertical exchanges due to iceberg melt and upwelling ($Q_j^{Xi}$), vertical advection ($Q_j^{Xv}$), and vertical mixing ($Q_j^{Xk}$). The forcing boundary conditions are the subglacial discharge, Q_sg(t), and the shelf temperature, T^s(z,t), and salinity, S^s(z,t).

We split the fjord (Fig. 1) into N vertically stacked layers, indexed by j from the surface (j=1) to the bottom (j=N). Each layer has a volume V_j, which is fixed in time. At present, the model assumes a cuboid-shaped fjord of length L and width W, such that V_j=WLH_j, where H_j is the thickness of layer j. Each layer also has a static iceberg area in contact with the ocean I_j and a time-varying temperature T_j(t) and salinity S_j(t) that can be considered as an average over the volume of the layer. The number of layers and their volumes are user-defined; the layers do not have to have equal volume or thickness (e.g. see Fig. 1), though the sum of the layer thicknesses has to equal the fjord depth. Fewer layers means less vertical resolution, whereas more layers limits the time step of the model. If the fjord has a sill, then the model will slightly adjust the thicknesses to ensure that the sill depth coincides with a layer boundary (e.g. Fig. 1), simplifying the formulation of the fjord–shelf exchange.

The layers experience time-dependent exchanges (Fig. 1) with the subglacial discharge-driven plumes (entrainment and intrusion), with the shelf (fjord–shelf exchange), with icebergs (melting and upwelling), and with adjacent layers (vertical mixing and vertical advection). Note that the model allows for multiple plumes because multiple marine-terminating glaciers can terminate in a single fjord. Other glacial fjord exchange processes and drivers of circulation, such as ambient melting of the glacier, surface fluxes, sea ice, winds, and tides, were considered lower priority and are not yet implemented, but this is not intended to imply that these processes are not important in some situations. Exchange fluxes are denoted by $Q_j^{Xy}$, where j denotes the layer where the exchange is happening, X denotes the flux quantity being exchanged (V for volume, T for heat, or S for salt) and y denotes the process driving exchange (p for plume, s for shelf, i for icebergs, k for vertical mixing, or v for vertical advection).

For fluxes that represent exchange between adjacent layers (iceberg-driven upwelling, vertical mixing, and vertical advection), a further line of explanation is required. These fluxes naturally take a general form $Q_{j+\mathrm{1},j}^{Xy}$, where j+1 and j are the layers exchanging. These fluxes can, however, be written in the single-subscript format previously introduced by noting that the net effect on layer j is the difference between the exchange with layer j+1 and the exchange with layer j−1; that is, $Q_j^{Xy}=Q_{j+\mathrm{1},j}^{Xy}-Q_{j,j-\mathrm{1}}^{Xy}$. In this expression, there is no exchange between the bottom layer and the sea floor or between the surface layer and atmosphere ($Q_{N+\mathrm{1},N}^{Xy}=Q_{\mathrm{1},\mathrm{0}}^{Xy}=\mathrm{0}$).

For sign convention, a volume flux $Q_j^{Vy}$ is defined as positive if it is adding volume to a layer. A volume flux between layers $Q_{j+\mathrm{1},j}^{Vy}$ is defined as positive if it removes volume from the deeper layer j+1 and adds it to the shallower layer j. To give some examples, $Q_{\mathrm{2}}^{Vp}$ is the volume exchange between the plume and layer 2, which is positive if the plume is adding volume to layer 2, and $Q_{\mathrm{5}}^{Vs}$ is the volume exchange between fjord layer 5 and the shelf, which is positive if it is directed from the shelf into the fjord.

Estimating many of the exchanges requires knowledge of the density differences between two water masses. To do this, we employ a linear equation of state $\mathit{\rho }(T,S)={\mathit{\rho }}_{\text{ref}}\left[\mathrm{1}+{\mathit{\beta }}_{\mathrm{S}}(S-S_{\text{ref}})-{\mathit{\beta }}_{\mathrm{T}}(T-T_{\text{ref}})\right]$. Here, ρ_ref is the reference density of water with salinity S_ref and temperature T_ref, and β_S and β_T are the haline contraction and thermal expansion coefficients, respectively. The relative buoyancy of two water masses a and b is then defined by

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{g}_{\mathrm{ab}}^{\prime }=g{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}-{\mathit{\rho }}_{\mathrm{b}}}{{\mathit{\rho }}_{\text{ref}}}}=g\left[{\mathit{\beta }}_{\mathrm{S}}\left(S_{\mathrm{a}}-S_{\mathrm{b}}\right)-{\mathit{\beta }}_{\mathrm{T}}\left(T_{\mathrm{a}}-T_{\mathrm{b}}\right)\right],\end{array}$$

in which g is the gravitational acceleration and the reference properties have cancelled out.

We now describe in turn how we parameterise exchange with the plume (Sect. 2.2), with the shelf (Sect. 2.3), with icebergs (Sect. 2.4), due to vertical mixing (Sect. 2.5), and due to vertical advection (Sect. 2.6). For reference, a full list of model variables and parameters, including their units, is given in Tables A1 and A2.

2.2 Exchange with the plume

If there is subglacial discharge emerging from beneath a marine-terminating glacier, this will drive an upwelling plume adjacent to the glacier that entrains water at depth and intrudes into the fjord at a level of neutral buoyancy (Cowton et al., 2015; Carroll et al., 2015; Stevens et al., 2016; Mankoff et al., 2016). Within FjordRPM, this is represented by a flux leaving the deeper layers and a flux entering a shallower layer (Fig. 1). To quantify these fluxes, we employ a now-standard approach that has been termed “plume-melt theory” by Jackson et al. (2022), consisting of the classic theory for buoyant plumes (Morton et al., 1956) coupled to the standard three-equation parameterisation for submarine melt (Holland and Jenkins, 1999). If there are multiple discharge plumes emerging from a single glacier or multiple glaciers with discharge plumes draining into a single fjord, the following approach is applied separately for each discharge plume. Note that submarine melting and melt-driven convection outside of subglacial discharge plumes are not currently accounted for in FjordRPM.

We assume a line plume geometry (Jackson et al., 2017) with a fixed plume width in the across-fjord direction of W_p. As the plume rises, plume volume, momentum, heat, and salt flux per unit width evolve according to

$$\begin{array}{ll}\text{(2a)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left(b_{\mathrm{p}}{u}_{\mathrm{p}}\right)={\mathit{\alpha }}_{\mathrm{p}}{u}_{\mathrm{p}}+m_{\mathrm{p}},\text{(2b)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left(b_{\mathrm{p}}{u}_{\mathrm{p}}^{\mathrm{2}}\right)=b_{\mathrm{p}}{g}_{\mathrm{p}}^{\prime }-C_{\mathrm{d}}{u}_{\mathrm{p}}^{\mathrm{2}},{\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left(b_{\mathrm{p}}{u}_{\mathrm{p}}{T}_{\mathrm{p}}\right)={\mathit{\alpha }}_{\mathrm{p}}{u}_{\mathrm{p}}T+m_{\mathrm{p}}{T}_{\mathrm{b}}\\ \text{(2c)}& {\displaystyle }& {\displaystyle }-C_{\mathrm{d}}^{\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_T{u}_{\mathrm{p}}\left(T_{\mathrm{p}}-T_{\mathrm{b}}\right),{\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left(b_{\mathrm{p}}{u}_{\mathrm{p}}{S}_{\mathrm{p}}\right)={\mathit{\alpha }}_{\mathrm{p}}{u}_{\mathrm{p}}S+m_{\mathrm{p}}{S}_{\mathrm{b}}\\ \text{(2d)}& {\displaystyle }& {\displaystyle }-C_{\mathrm{d}}^{\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_S{u}_{\mathrm{p}}\left(S_{\mathrm{p}}-S_{\mathrm{b}}\right),\end{array}$$

in which b_p is the plume width in the along-fjord direction, u_p is the plume (vertical) velocity, T_p is the plume temperature, S_p is the plume salinity, α_p is the plume entrainment coefficient, $g_{\mathrm{p}}^{\prime }=g\left[{\mathit{\beta }}_{\mathrm{S}}\left(S-S_{\mathrm{p}}\right)-{\mathit{\beta }}_{\mathrm{T}}\left(T-T_{\mathrm{p}}\right)\right]$ is the relative buoyancy of the fjord and plume, and T and S are the fjord temperature and salinity (see also Table A1). All other variables and parameters relate to the submarine melt rate induced by the plume, m_p, which is given by three equations that balance heat and salt flux through the ice–ocean boundary layer and maintain the boundary layer at the in situ freezing point (e.g. Jenkins, 2011):

$$\begin{array}{}\text{(3a)}& {\displaystyle }& {\displaystyle }{m}_{\mathrm{p}}l+m_{\mathrm{p}}{c}_{\mathrm{i}}\left(T_{\mathrm{b}}-T_{\mathrm{i}}\right)=c_{\mathrm{w}}{C}_{\mathrm{d}}^{\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_T{u}_{\mathrm{p}}\left(T_{\mathrm{p}}-T_{\mathrm{b}}\right),\text{(3b)}& {\displaystyle }& {\displaystyle }{m}_{\mathrm{p}}{S}_{\mathrm{b}}=C_{\mathrm{d}}^{\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_S{u}_{\mathrm{p}}\left(S_{\mathrm{p}}-S_{\mathrm{b}}\right),\text{(3c)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{b}}={\mathit{\lambda }}_{\mathrm{1}}{S}_{\mathrm{b}}+{\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}\left|z\right|,\end{array}$$

where T_b and S_b are the temperature and salinity of the ice–ocean boundary layer, l is the latent heat of ice melt, c_i and c_w are the heat capacities of ice and seawater, C_d is the plume–ice drag coefficient, Γ_T and Γ_S are heat and salt transfer coefficients, and λ_i denotes the constants that control the linearised freezing point (see also Table A2).

The boundary conditions for Eq. (2) are the fjord temperature T and salinity S, which are taken from the FjordRPM layers. The plume is initiated by a flux of subglacial discharge Q_sg emerging from beneath the glacier at the grounding line depth H_gl, with initial salinity S_sg=0 and temperature $T_{\text{sg}}={\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}{H}_{\text{gl}}$ (i.e. the in situ freezing point). The initial plume width and velocity are set by ensuring that the plume has a balance of buoyancy and momentum (e.g. Slater et al., 2016) and an initial flux of Q_sg, giving $b_{\mathrm{p}}={\left({\mathit{\alpha }}_{\mathrm{p}}{Q}_{\text{sg}}^{\mathrm{2}}/g_{\mathrm{p}}^{\prime }{W}_{\mathrm{p}}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{3}}$ and $u_{\mathrm{p}}=Q_{\text{sg}}/b_{\mathrm{p}}{W}_{\mathrm{p}}$.

At each FjordRPM time step, Eqs. (2) and (3) are numerically integrated using an Euler scheme on the vertical grid defined by the interfaces between the FjordRPM model layers (Fig. 1). To provide an example based on Eq. (2a), if the plume width, velocity, and melt rate at the interface between layers j+1 and j are $b_{\mathrm{p}}^{j+\mathrm{1},j}$, $u_{\mathrm{p}}^{j+\mathrm{1},j}$, and $m_{\mathrm{p}}^{j+\mathrm{1},j}$, then the plume volume flux at the next interface up, between layers j and j−1, is estimated by

$$\begin{array}{ll}{\displaystyle }{b}_{\mathrm{p}}^{j,j-\mathrm{1}}{u}_{\mathrm{p}}^{j,j-\mathrm{1}}=& {\displaystyle }{b}_{\mathrm{p}}^{j+\mathrm{1},j}{u}_{\mathrm{p}}^{j+\mathrm{1},j}\\ \text{(4)}& {\displaystyle }& {\displaystyle }+H_j\left({\mathit{\alpha }}_{\mathrm{p}}{u}_{\mathrm{p}}^{j+\mathrm{1},j}+m_{\mathrm{p}}^{j+\mathrm{1},j}\right).\end{array}$$

Beginning from the grounding line, this numerical integration continues to shallower layers until the plume finds neutral buoyancy and intrudes into the fjord. The layer of plume neutral buoyancy is defined as the layer where the plume–fjord relative buoyancy, $g_{\mathrm{p}}^{\prime }$, switches from positive to negative, indicating that the plume is negatively buoyant in that layer and will stop rising and intrude horizontally. If the plume reaches the fjord surface without becoming negatively buoyant, then it intrudes into the surface layer. In this manner, we obtain the plume velocity u_p and submarine melt rate m_p at each of the layer interfaces.

The volume flux per unit plume width entrained from the fjord into the plume is the first term on the right-hand side of Eq. (2a), while the associated heat and salt fluxes are the first terms on the right-hand side of Eqs. (2c) and (2d). For layers where the plume is rising, the exchange fluxes between the model layers and the plume (Fig. 1) are therefore given by

$$\begin{array}{}\text{(5a)}& {\displaystyle }& {\displaystyle }{Q}_j^{V_{\mathrm{p}}}=-{\mathit{\alpha }}_{\mathrm{p}}{W}_{\mathrm{p}}{u}_{\mathrm{p}}^{j+\mathrm{1},j}{H}_j,\text{(5b)}& {\displaystyle }& {\displaystyle }{Q}_j^{T_{\mathrm{p}}}=-{\mathit{\alpha }}_{\mathrm{p}}{W}_{\mathrm{p}}{u}_{\mathrm{p}}^{j+\mathrm{1},j}{H}_j{T}_j,\text{(5c)}& {\displaystyle }& {\displaystyle }{Q}_j^{S_{\mathrm{p}}}=-{\mathit{\alpha }}_{\mathrm{p}}{W}_{\mathrm{p}}{u}_{\mathrm{p}}^{j+\mathrm{1},j}{H}_j{S}_j,\end{array}$$

where we know $u_{\mathrm{p}}^{j+\mathrm{1},j}$ from the numerical integration. For any layers that are deeper than the grounding line or shallower than the layer where the plume intrudes into the fjord, all plume fluxes are 0.

In the layer where the plume intrudes into the fjord, say j₀, the flux from the plume into the fjord consists of the subglacial discharge, the submarine meltwater, and all of the water that has been entrained into the plume. Thus, the volume flux from the plume into the layer of neutral buoyancy is

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }{Q}_{j_{\mathrm{0}}}^{V_{\mathrm{p}}}=Q_{\text{sg}}+W_{\mathrm{p}}\sum _j{H}_j\left[{\mathit{\alpha }}_{\mathrm{p}}{u}_{\mathrm{p}}^{j+\mathrm{1},j}+m_{\mathrm{p}}^{j+\mathrm{1},j}\right],\end{array}$$

where the sum runs over all layers j between the grounding line and the layer of neutral buoyancy. To obtain the heat flux exchange, we apply the same principle, but we first use Eq. (3a) to rewrite Eq. (2c) as

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}}\left(b_{\mathrm{p}}{u}_{\mathrm{p}}{T}_{\mathrm{p}}\right)={\mathit{\alpha }}_{\mathrm{p}}{u}_{\mathrm{p}}T+m_{\mathrm{p}}{T}_{\text{eff}},\end{array}$$

where

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{T}_{\text{eff}}=T_{\mathrm{b}}-{\displaystyle \frac{l}{c_{\mathrm{w}}}}-{\displaystyle \frac{c_{\mathrm{i}}}{c_{\mathrm{w}}}}\left(T_{\mathrm{b}}-T_{\mathrm{i}}\right)\end{array}$$

is the effective meltwater temperature (e.g. Jenkins, 2011). Noting the relative size of the terms, we furthermore approximate $T_{\text{eff}}\approx -l/c_{\mathrm{w}}$. The heat flux from the plume into the layer of neutral buoyancy can then be calculated as

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{Q}_{j_{\mathrm{0}}}^{T_{\mathrm{p}}}=Q_{\text{sg}}{T}_{\text{sg}}+W_{\mathrm{p}}\sum _j{H}_j\left[{\mathit{\alpha }}_{\mathrm{p}}{u}_{\mathrm{p}}^{j+\mathrm{1},j}{T}_j+m_{\mathrm{p}}^{j+\mathrm{1},j}{T}_{\text{eff}}\right].\end{array}$$

Lastly, the salt flux from the plume into the layer of neutral buoyancy can be calculated as

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{Q}_{j_{\mathrm{0}}}^{S_{\mathrm{p}}}=W_{\mathrm{p}}\sum _j{H}_j\left[{\mathit{\alpha }}_{\mathrm{p}}{u}_{\mathrm{p}}^{j+\mathrm{1},j}{S}_j\right],\end{array}$$

where we have used the fact that the salinity of subglacial discharge, S_sg, is 0 and the fact that the last two terms on the right-hand side of Eq. (2d) cancel according to Eq. (3b). The sums in Eqs. (9) and (10) again run over all layers j between the grounding line and the layer of neutral buoyancy.

By these definitions, the plume removes water from the deeper layers, mixes it with subglacial discharge and submarine meltwater, and puts the resulting mixture back into a layer closer to the surface, representing the upwelling and transformation of fjord waters by the plume and driving an overturning circulation. While the plume exchange fluxes in a given layer may be large, note that because all of the volume flux entrained from deeper layers is put back into a shallower layer, the net volume flux over all layers associated with the plume is simply the subglacial discharge and the submarine melt flux; that is,

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }\sum _{j=\mathrm{1}}^N{Q}_j^{V_{\mathrm{p}}}=Q_{\text{sg}}+\sum _{j=\mathrm{1}}^N{W}_{\mathrm{p}}{H}_j{m}_{\mathrm{p}}^{j+\mathrm{1},j}\equiv Q_{\text{sg}}+Q_{\text{sm}}.\end{array}$$

Similarly, all of the heat and salt entrained from deeper layers is put back into a shallower layer, and the net heat and salt flux associated with the plume come only from the subglacial discharge and submarine melting.

2.3 Exchange with the shelf

We assume that exchange between the shelf and fjord is driven by pressure gradients. The pressure difference, Δp_j, between the shelf and fjord in layer j may be quantified as

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle \frac{\mathrm{\Delta }{p}_j}{{\mathit{\rho }}_{\text{ref}}}}=g\mathrm{\Delta }\mathit{\eta }+\int {\displaystyle \frac{g\mathrm{\Delta }\mathit{\rho }\left(z\right)}{{\mathit{\rho }}_{\text{ref}}}}\mathrm{d}z,\end{array}$$

where Δη is the sea surface height difference between the shelf and fjord, Δρ(z) is the density difference between the shelf and fjord, and the integral runs over all depths from the surface down to the midpoint of layer j. The integrand in Eq. (12) is the relative buoyancy between the shelf and fjord, given by

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }{g}_{\mathrm{s}j}^{\prime }=g\left[{\mathit{\beta }}_{\mathrm{S}}\left(S_j^{\mathrm{s}}-S_j\right)-{\mathit{\beta }}_{\mathrm{T}}\left(T_j^{\mathrm{s}}-T_j\right)\right],\end{array}$$

where $S_j^{\mathrm{s}}$ and $T_j^{\mathrm{s}}$ are the shelf temperature and salinity profiles mapped onto the FjordRPM grid (Fig. 1). That is,

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }{S}_j^{\mathrm{s}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{H_j}}\underset{j}{\int }{S}^{\mathrm{s}}(z,t)\mathrm{d}z,\end{array}$$

and the integral runs over the depth range covered by layer j. An equivalent expression applies for temperature. If we now define a quantity ϕ, given for each layer by

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }{\mathit{\varphi }}_j=g_{\mathrm{s}j}^{\prime }{H}_j/\mathrm{2}+\sum _{k=\mathrm{1}}^{j-\mathrm{1}}{g}_{\mathrm{s}k}^{\prime }{H}_k,\end{array}$$

then Eq. (12) may be written as $\mathrm{\Delta }{p}_j/{\mathit{\rho }}_{\text{ref}}=g\mathrm{\Delta }\mathit{\eta }+{\mathit{\varphi }}_j$. That is, ϕ_j is the integral in Eq. (12) evaluated using the model layers. The shelf-to-fjord pressure gradient is then finally estimated by dividing by the fjord length, L. We assume that the sea surface height difference adjusts quickly to ensure that the net flux into the fjord is 0 and thus treat the barotropic pressure term separately below.

To estimate the shelf–fjord exchange, previous studies have considered cases where this pressure gradient is (i) balanced by friction and mixing (Geyer and MacCready, 2014; Sanchez et al., 2023); (ii) in some form of geostrophic balance (Zhao et al., 2021); or (iii) under hydraulic control (Nilsson et al., 2023). The relevant balance for a given fjord will depend on the fjord and sill geometry, the stratification, and the circulation. We have derived our exchange parameterisation with regimes (i) and (ii) in mind.

To proceed, we consider the along-fjord momentum balance, neglecting the acceleration and advection terms. Under regime (i), the Coriolis term is additionally neglected, and the pressure gradient can be balanced by vertical stress divergence, obtaining an exchange velocity that scales as $\frac{H^{\mathrm{2}}}{A}\frac{{\mathit{\varphi }}_j}{L}$, in which H is the fjord depth and A is a vertical eddy viscosity (e.g. Geyer and MacCready, 2014). Alternatively, the pressure gradient can be said to be balanced by more general sources of friction, obtaining an exchange velocity that scales as $\frac{\mathrm{1}}{r}\frac{{\mathit{\varphi }}_j}{L}$, in which $\mathrm{1}/r$ is a frictional timescale (Sanchez et al., 2023). Under regime (ii), Zhao et al. (2021) argued that boundary currents should set up an along-fjord pressure head of similar magnitude to the across-fjord pressure head, which leads to a velocity scaling like $\frac{\mathrm{1}}{f}\frac{{\mathit{\varphi }}_j}{L}$, where f is the Coriolis parameter. Generalising these cases, we take an exchange velocity of the form $C_{\mathrm{0}}\frac{{\mathit{\varphi }}_j}{L}$ and treat C₀ as an empirical tunable parameter.

We then define the fjord–shelf volume exchange fluxes for above-sill layers as

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }{Q}_j^{Vs}=W{H}_j\left(u_{\mathrm{b}}+C_{\mathrm{0}}{\displaystyle \frac{{\mathit{\varphi }}_j}{L}}\right),\end{array}$$

where WH_j is the cross-sectional area of contact between layer j and the shelf and u_b is a barotropic velocity that ensures overall conservation of the fjord volume (see below). C₀ controls the strength of the exchange and accounts for the fjord–shelf balance that is being assumed. For layers that are below the sill (Fig. 1), there is no direct exchange between that layer and the shelf.

Regarding the barotropic velocity u_b, note that if the sum of all the volume fluxes into the fjord is not 0, there will be a change in the total fjord volume. This could be accounted for by dynamically tracking the sea surface height of the fjord relative to the shelf and its impact on fjord-to-shelf pressure gradients. However, the adjustment timescale for the sea surface height is very fast (on the order of minutes), and therefore tracking sea surface height requires a very small model time step that limits the rest of the model. In addition, the sea surface height rapidly approaches a steady state in which the fjord–shelf volume fluxes have adjusted to ensure there is no net volume flux into the fjord. On the timescales of the model time step (hours to days), it is therefore a very good approximation to choose u_b in Eq. (16) to ensure there is no net volume flux into the fjord. This is equivalent to the “rigid lid” approach sometimes applied in general circulation models. Other than the shelf fluxes, the only contribution to net volume change is from the plume fluxes, where the net volume flux comprises the subglacial discharge and submarine meltwater and is given by Eq. (11). Thus, if there are N_above above-sill layers, conservation of overall fjord volume means

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }{Q}_{\text{sg}}+Q_{\text{sm}}+\sum _{j=\mathrm{1}}^{N_{\text{above}}}{Q}_j^{Vs}=\mathrm{0},\end{array}$$

from which, combining with Eq. (16), we can obtain the barotropic velocity, u_b, as

$$\begin{array}{ll}{\displaystyle }{u}_{\mathrm{b}}=& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{W{\sum }_{j=\mathrm{1}}^{N_{\text{above}}}{H}_j}}\\ \text{(18)}& {\displaystyle }& {\displaystyle }\times \left(Q_{\text{sg}}+Q_{\text{sm}}+{\displaystyle \frac{W{C}_{\mathrm{0}}}{L}}\sum _{j=\mathrm{1}}^{N_{\text{above}}}{H}_j{\mathit{\varphi }}_j\right).\end{array}$$

The associated heat and salt fluxes are obtained by noting that if, for layer j, the volume flux is directed out of the fjord, then the salinity associated with the volume flux should be that of the layer, while if the volume flux is directed into the fjord, then the salinity should be that of the shelf. Thus, the salt fluxes associated with exchange with the shelf are given by

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }{Q}_j^{Ss}=\left\{\begin{array}{ll}{Q}_j^{Vs}{S}_j& \text{if }{Q}_j^{Vs}<\mathrm{0}\\ Q_j^{Vs}{S}_j^{\mathrm{s}}& \text{if }{Q}_j^{Vs}\ge \mathrm{0}\end{array}\right.,\end{array}$$

with similar expressions for heat, completing the definition of the shelf fluxes.

2.4 Exchange with icebergs

Iceberg melting freshens and cools fjords (Hager et al., 2024; Abib et al., 2024), and the input of the buoyant meltwater gives melt-driven convective plumes rising up the sides of icebergs (Cenedese and Straneo, 2023). These plumes could be represented by plume-melt theory (Magorrian and Wells, 2016), as for the subglacial discharge plumes in Sect. 2.2, but we would require many such plumes and would have to decide, on the basis of iceberg characteristics, at what depth to initiate these plumes. We are not convinced that such complexity is merited within a reduced-physics model, and so we have aimed to take a simpler approach without losing the key physics.

Starting with iceberg melting, numerous parameterisations exist for iceberg melt rate (e.g. Neshyba and Josberger, 1980; FitzMaurice et al., 2017; Cenedese and Straneo, 2023). Ideally, we would use the three-equation melt rate parameterisation (Eq. 3); however, it is not straightforward to estimate the relative velocity of the icebergs and ocean that would be required. We therefore instead take the simple approach of assuming that the melt rate of icebergs in layer j is proportional to the thermal forcing of that layer, with a constant of proportionality M₀. Once scaled for the layer iceberg surface area, I_j, the iceberg melt flux is given by

$$\begin{array}{}\text{(20)}& {\displaystyle }{\displaystyle }{Q}_j^{\text{melt}}=M_{\mathrm{0}}\left(T_j-T_j^f\right)I_j,\end{array}$$

where M₀ is a constant and $T_j^f={\mathit{\lambda }}_{\mathrm{1}}{S}_j+{\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}{z}_j$ is the linearised in situ freezing point, which depends on the layer salinity S_j and the mean depth of the layer, $z_j=H_j/\mathrm{2}+{\sum }_{n=\mathrm{1}}^{j-\mathrm{1}}{H}_n$.

To estimate the resulting upwelling flux, we follow scalings from Magorrian and Wells (2016) for melt-driven convection. Leaving the details to Appendix A, the scale of the (vertical) convection velocity resulting from iceberg melting in layer j can be taken to be

$$\begin{array}{}\text{(21)}& {\displaystyle }{\displaystyle }{v}_j={\left[{\displaystyle \frac{Q_j^{\text{melt}}{g}_{j,\text{melt}}^{\prime }{H}_j}{{\mathit{\alpha }}_i{I}_j}}\right]}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where α_i is the entrainment coefficient and $g_{j,\text{melt}}^{\prime }$ is the relative buoyancy of the layer and meltwater, given by

$$\begin{array}{}\text{(22)}& {\displaystyle }{\displaystyle }{g}_{j,\text{melt}}^{\prime }=g\left[{\mathit{\beta }}_S{S}_j-{\mathit{\beta }}_{\mathrm{T}}\left(T_j-T_{\text{eff}}\right)\right],\end{array}$$

in which $T_{\text{eff}}=-l/c_{\mathrm{w}}$ is the effective meltwater temperature as in Sect. 2.2.

By making an analogy with a line plume of width $W_j=I_j/H_j$, the entrainment flux into the rising melt-driven convection plume is then

$$\begin{array}{}\text{(23)}& {\displaystyle }{\displaystyle }{Q}_j^{\text{ent}}={\mathit{\alpha }}_i{v}_j{H}_j{W}_j={\mathit{\alpha }}_i{v}_j{I}_j.\end{array}$$

In a stratified fjord, however, the melt-driven convection cell will reach neutral buoyancy some height after initiating. The length scale that controls this height can be estimated, again following scalings in Magorrian and Wells (2016) and detailed in Appendix A, by

$$\begin{array}{}\text{(24)}& {\displaystyle }{\displaystyle }{l}_{j+\mathrm{1},j}^{\text{ice}}={\displaystyle \frac{v_{j+\mathrm{1}}^{\mathrm{2}}}{H_{j+\mathrm{1}}}}{\displaystyle \frac{H_j+H_{j+\mathrm{1}}}{\mathrm{2}{g}_{j+\mathrm{1},j}^{\prime }}}.\end{array}$$

If $l_{j+\mathrm{1},j}^{\text{ice}}>H_{j+\mathrm{1}}$, then any melt-driven convection cell that initiates in layer j+1 will upwell to layer j. If, however, $l_{j+\mathrm{1},j}^{\text{ice}}=H_{j+\mathrm{1}}/\mathrm{2}$ (say), then only melt-driven convection cells initiating in the upper half of layer j+1 would upwell over the layer boundary to layer j. Assuming there is an equal probability of melt-driven convection cells initiating at any point in the layer, the resulting upwelling flux would be half of that in Eq. (23). The final iceberg upwelling volume flux from layer j+1 to j is thus given by Eq. (23) scaled by a fraction f^ice that accounts for how far iceberg meltwater is able to upwell:

$$\begin{array}{}\text{(25)}& {\displaystyle }{\displaystyle }{Q}_{j+\mathrm{1},j}^{Vi}=\text{min}\left[\mathrm{1},{\displaystyle \frac{l_{j+\mathrm{1},j}^{\text{ice}}}{H_{j+\mathrm{1}}}}\right]Q_{j+\mathrm{1}}^{\text{ent}}\equiv f_{j+\mathrm{1},j}^{\text{ice}}{Q}_{j+\mathrm{1}}^{\text{ent}}.\end{array}$$

Note that $f_{N+\mathrm{1},N}^{\text{ice}}=f_{\mathrm{1},\mathrm{0}}^{\text{ice}}=\mathrm{0}$ because there is no upwelling into the bottom layer or from the surface layer. The associated temperature and salt fluxes are given by multiplying by the temperature and salinity of the layer from which the volume is upwelling.

Having accounted for the upwelling driven by iceberg melting, we return to the melting itself. The impact of melting ice on ambient ocean waters is equivalent to adding freshwater of effective temperature T_eff to the ambient (Jenkins, 2011). In addition, the volumes of iceberg meltwater are much smaller than the layer volumes, such that we may treat the meltwater as a virtual flux, meaning that we account for its impact on layer temperature and salinity but not on layer volume. Under these assumptions, and if all of the iceberg meltwater in layer j stays in layer j, then the layer heat flux associated with the iceberg melt flux $Q_j^{\text{melt}}$ is $-Q_j^{\text{melt}}l/c_{\mathrm{w}}$, and the salt flux is $-Q_j^{\text{melt}}{S}_j$.

Putting it all together, the net volume flux into layer j due to iceberg melting is

$$\begin{array}{}\text{(26)}& {\displaystyle }{\displaystyle }{Q}_j^{Vi}=Q_{j+\mathrm{1},j}^{Vi}-Q_{j,j-\mathrm{1}}^{Vi},\end{array}$$

where the terms on the right-hand side are given by Eq. (25). The net heat flux into layer j is

$$\begin{array}{ll}{\displaystyle }{Q}_j^{Ti}=& {\displaystyle }{Q}_{j+\mathrm{1},j}^{Vi}{T}_{j+\mathrm{1}}-Q_{j,j-\mathrm{1}}^{Vi}{T}_j-f_{j+\mathrm{1},j}^{\text{ice}}{Q}_{j+\mathrm{1}}^{\text{melt}}{\displaystyle \frac{l}{c_{\mathrm{w}}}}\\ \text{(27)}& {\displaystyle }& {\displaystyle }-\left(\mathrm{1}-f_{j,j-\mathrm{1}}^{\text{ice}}\right)Q_j^{\text{melt}}{\displaystyle \frac{l}{c_{\mathrm{w}}}},\end{array}$$

where the first term is the upwelling of ambient water from the layer below, the second term is the upwelling of ambient water to the layer above, the third term is the upwelling of meltwater from the layer below, and the fourth term is the meltwater generated in the present layer that stays in the present layer. Finally, the net salt flux into layer j is

$$\begin{array}{ll}{\displaystyle }{Q}_j^{Si}=& {\displaystyle }{Q}_{j+\mathrm{1},j}^{Vi}{S}_{j+\mathrm{1}}-Q_{j,j-\mathrm{1}}^{Vi}{S}_j-f_{j+\mathrm{1},j}^{\text{ice}}{Q}_{j+\mathrm{1}}^{\text{melt}}{S}_j\\ \text{(28)}& {\displaystyle }& {\displaystyle }-\left(\mathrm{1}-f_{j,j-\mathrm{1}}^{\text{ice}}\right)Q_j^{\text{melt}}{S}_j.\end{array}$$

Note that, by definition, ${\sum }_{j=\mathrm{1}}^N{Q}_j^{Vi}=\mathrm{0}$; hence, there is no fjord-scale volume flux associated with iceberg melting. In doing the same sum for heat and salt flux, the first two terms on the right-hand sides of Eqs. (27) and (28) cancel because upwelling vertically redistributes heat and salt but does not change the fjord-scale content. The second two terms on the right-hand sides of Eqs. (27) and (28) do not, however, cancel, giving a net cooling and freshening due to iceberg melting.

2.5 Vertical mixing

Exchange between layers occurs due to vertical mixing at the layer interfaces. We assume that this causes an exchange of heat and salt but no net exchange of volume and calculate the vertical turbulent tracer fluxes as the product of an eddy diffusivity, K_z, and the vertical gradient of the tracer. Once scaled for the area of contact between layers j+1 and j, this gives a salt exchange of

$$\begin{array}{}\text{(29)}& {\displaystyle }{\displaystyle }{Q}_{j+\mathrm{1},j}^{Sk}={\mathrm{WLK}}_z{\displaystyle \frac{\partial S}{\partial z}}=\mathrm{2}{\mathrm{WLK}}_z{\displaystyle \frac{S_{j+\mathrm{1}}-S_j}{H_{j+\mathrm{1}}+H_j}},\end{array}$$

and an equivalent expression applies for temperature (we make no distinction between the eddy diffusivity of heat and salt). This eddy diffusivity is estimated as a function of the local Richardson number, Ri, using the KPP scheme (Large et al., 1994):

$$\begin{array}{}\text{(30)}& {\displaystyle }{\displaystyle }{K}_z=K_{\mathrm{b}}+\left\{\begin{array}{ll}{K}_{\mathrm{0}}& \mathit{\text{Ri}}\le \mathrm{0}\\ K_{\mathrm{0}}{\left[\mathrm{1}-{\left(\mathit{\text{Ri}}/{\mathit{\text{Ri}}}_{\mathrm{0}}\right)}^{\mathrm{2}}\right]}^{\mathrm{3}}& {\mathit{\text{Ri}}}_{\mathrm{0}}>\mathit{\text{Ri}}>\mathrm{0}\\ \mathrm{0}& \mathit{\text{Ri}}\ge {\mathit{\text{Ri}}}_{\mathrm{0}}\end{array}\right.,\end{array}$$

where Ri₀=0.7 is the Richardson number above which shear-driven mixing is suppressed, $K_{\mathrm{0}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ m² s⁻¹ scales the shear-driven mixing, and $K_{\mathrm{b}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹ is the background vertical mixing associated with internal waves (Table A2). We adopt these default values from Large et al. (1994), but we note evidence from Greenland fjords that smaller values may be more appropriate (Bendtsen et al., 2021). The Richardson number takes the usual definition $\mathit{\text{Ri}}=N^{\mathrm{2}}/(\partial u/\partial z{)}^{\mathrm{2}}$ and is estimated at the interface between layers j+1 and j as

$$\begin{array}{ll}{\displaystyle }\mathit{Ri}& {\displaystyle }={\displaystyle \frac{\mathrm{2}{g}_{j+\mathrm{1},j}^{\prime }/(H_{j+\mathrm{1}}+H_j)}{{\left[\mathrm{2}\left(u_{j+\mathrm{1}}-u_j\right)/\left(H_{j+\mathrm{1}}+H_j\right)\right]}^{\mathrm{2}}}},\\ \text{(31)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{g_{j+\mathrm{1},j}^{\prime }\left(H_{j+\mathrm{1}}+H_j\right)}{\mathrm{2}{\left(u_{j+\mathrm{1}}-u_j\right)}^{\mathrm{2}}}},\end{array}$$

where $g_{j+\mathrm{1},j}^{\prime }$ is the relative buoyancy of the two layers according to Eq. (1). Lastly, the horizontal velocity scale for each layer, u_j, is estimated using the plume and shelf exchange fluxes as

$$\begin{array}{}\text{(32)}& {\displaystyle }{\displaystyle }{u}_j={\displaystyle \frac{Q_j^{Vp}-Q_j^{Vs}}{\mathrm{2}W{H}_j}}.\end{array}$$

By these definitions, tracer vertical mixing is inhibited by stratification between layers but enhanced by velocity shear driven by the exchange between layers and the plume or shelf. Like for the iceberg fluxes, these layer-to-layer tracer fluxes can be converted into net fluxes for each layer using

$$\begin{array}{}\text{(33)}& {\displaystyle }{\displaystyle }{Q}_j^{Sk}=Q_{j+\mathrm{1},j}^{Sk}-Q_{j,j-\mathrm{1}}^{Sk},\end{array}$$

and an equivalent expression applies for net heat fluxes.

2.6 Conservation of layer volume

By the described parameterisations, there are volume fluxes into or out of the layers due to the plume, shelf, and icebergs. If, for a given layer, these fluxes do not sum to 0, then in a model with fixed layer volumes, there must be an additional volume flux that ensures conservation of layer volume. By analogy with models that calculate vertical velocities by imposing incompressibility, we take this additional volume flux to be vertical advection between layers. As an illustrative example, consider a fjord with a plume and a shallow sill that significantly restricts the inflow and outflow. The water that is upwelled by the plume is not able to be replaced directly (i.e. horizontally) by inflow from the shelf. There must therefore be downwelling (or reflux) within the fjord to replace the water entrained into the plume. This downwelling would be represented by the vertical flux now described.

Suppose that the net volume flux imbalance in layer j by the parameterisations described so far is $Q_j^{\text{imb}}$, given by

$$\begin{array}{}\text{(34)}& {\displaystyle }{\displaystyle }{Q}_j^{\text{imb}}=Q_j^{Vp}+Q_j^{Vs}+Q_j^{Vi}.\end{array}$$

For the surface layer (j=1), and because there can be no volume flux to the atmosphere, we ensure conservation of volume by imposing a flux from the second layer into the surface layer given by $Q_{\mathrm{2},\mathrm{1}}^{Vv}=-Q_{\mathrm{1}}^{\text{imb}}$. To ensure conservation of volume for the second layer, we must then impose a flux from the third layer into the second layer given by $Q_{\mathrm{3},\mathrm{2}}^{Vv}=-Q_{\mathrm{2}}^{\text{imb}}+Q_{\mathrm{2},\mathrm{1}}^{Vv}$. We can continue iterating to get the required vertical fluxes for all boxes with the general expression

$$\begin{array}{}\text{(35)}& {\displaystyle }{\displaystyle }{Q}_{j+\mathrm{1},j}^{Vv}=-\sum _{n=\mathrm{1}}^j{Q}_n^{\text{imb}}.\end{array}$$

The associated heat and salt fluxes depend on the direction of the flux. For salt flux, this is written as

$$\begin{array}{}\text{(36)}& {\displaystyle }{\displaystyle }{Q}_{j+\mathrm{1},j}^{Sv}=\left\{\begin{array}{ll}{Q}_{j+\mathrm{1},j}^{Vv}{S}_{j+\mathrm{1}}& \text{if }{Q}_{j+\mathrm{1},j}^{Vv}\ge \mathrm{0}\\ Q_{j+\mathrm{1},j}^{Vv}{S}_j& \text{if }{Q}_{j+\mathrm{1},j}^{Vv}<\mathrm{0}\end{array}\right..\end{array}$$

That is, if the flux is directed from layer j+1 to j, then the relevant salinity is S_j+1, whereas if the flux is directed from layer j to j+1, the relevant salinity is S_j. As for icebergs and vertical mixing, these layer-to-layer fluxes can finally be converted to net fluxes for each layer using

$$\begin{array}{}\text{(37)}& {\displaystyle }{\displaystyle }{Q}_j^{Sv}=Q_{j+\mathrm{1},j}^{Sv}-Q_{j,j-\mathrm{1}}^{Sv},\end{array}$$

with equivalent expressions applying for volume and heat fluxes.

2.7 Evolution equations

Having now defined all of the required fluxes (Fig. 1), we turn to the equations giving the time evolution of the model. Because layer volumes are fixed in time, the only evolution equations are for the salt and heat content of each layer. For layer salt content, we have

$$\begin{array}{}\text{(38)}& {\displaystyle }{\displaystyle }{V}_j{\displaystyle \frac{\mathrm{d}{S}_j}{\mathrm{d}t}}=Q_j^{Sp}+Q_j^{Ss}+Q_j^{Si}+Q_j^{Sk}+Q_j^{Sv},\end{array}$$

where the terms are, respectively, the salt fluxes associated with the plume(s) (Eqs. 5c and 10), shelf (Eq. 19), icebergs (Eq. 28), vertical mixing (Eq. 33), and vertical advection (Eq. 37). For layer heat content, we have the equivalent

$$\begin{array}{}\text{(39)}& {\displaystyle }{\displaystyle }{V}_j{\displaystyle \frac{\mathrm{d}{T}_j}{\mathrm{d}t}}=Q_j^{Tp}+Q_j^{Ts}+Q_j^{Ti}+Q_j^{Tk}+Q_j^{Tv}.\end{array}$$

These evolution equations, together with all of the flux definitions throughout Sect. 2.2–2.6, form a closed set of equations that can be numerically integrated in time.

2.8 Density inversions

Before describing the numerics of FjordRPM, we note one last way in which layer temperature and salinity can evolve, which is that we enforce total mixing of two adjacent layers if the upper layer becomes denser than the lower layer. Physically, this represents convection. Specifically, at each time step and for each pair of adjacent layers j+1 and j, we check the buoyancy jump $g_{j+\mathrm{1},j}^{\prime }$ between the layers using Eq. (1). For any pair of layers having $g_{j+\mathrm{1},j}^{\prime }<\mathrm{0}$, we reset the salinity of the layers to

$$\begin{array}{}\text{(40)}& {\displaystyle }{\displaystyle }{S}_{j+\mathrm{1}}=S_j={\displaystyle \frac{S_{j+\mathrm{1}}{V}_{j+\mathrm{1}}+S_j{V}_j}{V_{j+\mathrm{1}}+V_j}},\end{array}$$

with an equivalent equation for temperature.

3 Numerical implementation

3.1 Code structure and method of solution

FjordRPM is coded in MATLAB, though future plans include a translation to an open-source language such as Python and to Fortran for ease of coupling with earth system models. We make use of a number of structures to efficiently carry around the required variables: “p” contains all the model parameters, “f” contains the forcing boundary conditions, “a” contains the initial conditions, and “s” contains the solution. The time vector for the model is denoted “t”. Full lists of the variables in these structures are included for reference in Tables A1 and A2. The code is arranged in a modular fashion (Fig. 2) and proceeds as follows. The top-level routine (run_model) takes as input the structures p, t, f, and a and outputs the solution s. Within this routine, the code first makes a number of checks on the inputs (check_inputs), such as ensuring they have the correct dimensions. The solution variables are initialised and the initial conditions read in (initialise_variables). The forcing boundary conditions may be provided on a depth and time grid that differs from the model layers and time stepping, so some spatial and temporal interpolation of the forcings may be required (bin_forcings).

Figure 2 Schematic of the FjordRPM code.

The model then enters the main time-stepping loop (Fig. 2). At each time step t_i, with dynamic variables T_j(t_i) and S_j(t_i), the first action is to check for and resolve any density inversions (homogenise_unstable_layers; Sect. 2.8). The code then uses the temperature and salinity of the layers at the current time step [T_j(t_i),S_j(t_i)], together with all of the parameters and variables that do not vary in time, to calculate the heat and salt fluxes described in Sect. 2.2–2.6 (compute_fluxes and functions therein). With these fluxes, the model makes a forward Euler time step of Eqs. (38) and (39) to obtain the temperature and salinity $\left[T_j\left(t_{i+\mathrm{1}}\right),S_j\left(t_{i+\mathrm{1}}\right)\right]$ of the layers at the next time step (step_solution_forwards). The final stage of the time-stepping loop is to check that the sum of the volume flux exchange for each layer, which should be 0, is below a certain tolerance (check_solution). If it is not, which only happens if the model becomes unstable (see Sect. 3.3), the simulation is aborted. Once the time-stepping loop has finished, a final function (get_final_output) processes the solution into a final form for saving as output.

A special note is required here for the plume fluxes. Variability in the plume dynamics is determined by the subglacial discharge and fjord water properties. Because subglacial discharge and fjord water properties will generally vary only on timescales of days, while the model time step may be several hours, it is not essential that we “renew” the plume dynamics every single time step. Indeed, integrating the plume-melt theory (Eqs. 2 and 3) is one of the slowest parts of the model, and so updating the plume dynamics only every so often can substantially speed up the model. To implement this, note that the exchange fluxes between the plume and fjord (Eqs. 5, 6, 9, and 10) depend on the plume velocity, the induced submarine melt rate, and the fjord properties. Thus, we can use the same plume velocity and induced submarine melt rate for a number of time steps while using the fjord properties from the current time step. Essentially, this decouples the dynamics of the plume from the fjord properties for a number of time steps but maintains conservation of heat and salt because the temperature and salinity that enter the plume fluxes are always from the current time step. The parameter run_plume_every defines the number of time steps between updates of the plume dynamics.

3.2 Required inputs

As inputs, the model first requires specification of all of the physical parameters contained in the structure “p”. Default values for all of these parameters are provided in a function default_parameters; some, like the heat capacity of seawater c_w, are unlikely to need to be changed, while others, such as the shelf exchange parameter C₀, are more likely to be varied by the user. The geometry of the problem must also be specified, including the length, width, and depth of the fjord, the presence and depth of a possible sill, and the grounding line depth of any glaciers with plumes.

The user must also specify the number and thickness of the model layers. The layers do not all have to have the same thickness, so the user could choose to have higher vertical resolution at depths of interest. A large number of layers gives the solution better vertical resolution but will require a shorter time step (see next section) and increase the computation time, though this is not normally an issue unless doing an ensemble of hundreds of simulations. We recommend roughly 20–80 layers, giving a vertical resolution of O(10 m) in Greenland's fjords. When there is a sill, the model may make a small adjustment to the layer thicknesses (in the initialise_variables routine) to ensure that the sill depth lies at a layer boundary, simplifying the calculation of the shelf fluxes.

The forcing boundary conditions must be provided in the structure “f”. These include the subglacial discharge for each plume as a function of time and the shelf temperature and salinity as a function of depth and time. Finally, the initial conditions must be specified in the structure “a”. These include the initial temperature and salinity of the layers and also the surface area of icebergs in each layer. The latter is, for now, classified as an initial condition rather than a forcing because it does not vary in time.

3.3 Stability of time-stepping

The last required input is the time vector, “t”, encompassing the duration and time step for the simulation. Like most models, a time step that is too large can lead to numerical instability. The complexity of the right-hand side of Eqs. (38) and (39), however, makes it difficult to derive a general limit on the stable time step. If we take into account only the mixing fluxes, Eq. (39) becomes a 1D heat equation with a spatially and temporally variable diffusivity. With a forward Euler scheme, as employed here, stable time steps of the 1D heat equation must satisfy $\mathrm{\Delta }t\le (\mathrm{\Delta }z{)}^{\mathrm{2}}/\mathrm{2}{K}_z$, where Δz would be our layer thickness and K_z is the eddy diffusivity of heat. With K_z at a reasonable maximum possible value ($\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ m² s⁻¹, Sect. 2.5) and Δz=10 m, the maximum stable time step is around 0.1 d. However, in FjordRPM, this condition guarantees neither stability nor instability, as K_z is a function of the Richardson number and only reaches such a high value in a few layers and/or at a few time steps and because vertical mixing is just one of the total five fluxes.

Another handle on stability is the fraction of a layer volume that is exchanged in a given time step. For example, if a vertical volume flux is Q (Sect. 2.6), then the fraction of layer volume V that is exchanged in a time step Δt due to this flux is $Q\mathrm{\Delta }t/V$. Writing the volume flux as Q=WLv, with v the vertical velocity, and the volume as V=WLΔz, with Δz the layer thickness, the fraction of layer volume exchanged in a time step can be expressed as $v\mathrm{\Delta }t/\mathrm{\Delta }z$, which is an advective Courant–Friedrichs–Lewy condition. Investigations have shown that instability is often associated with this fraction exceeding ∼0.5; and with a typical vertical velocity scale of 10⁻⁴ m s⁻¹ and a layer thickness Δz=10 m, this gives a maximum stable time step of around 0.5 d. But again, this condition is not always relevant, as instability can also arise from horizontal fluxes or the vertical mixing or from interactions between these fluxes. Furthermore, we do not have these fluxes in advance of running the simulation, so this does not provide an easy way of choosing the time step.

In practice, the model is sufficiently cheap to run such that the pragmatic way to choose the time step is to try a number of values and check the solution for instability. As a guide, we have found that a stable time step is on the order of 0.1–4 d for realistic forcings and layer thicknesses of 10 m. Such simulations take no more than a few seconds per model year on a laptop.

3.4 Character of the model

Before moving onto example results, we provide a few insights into the character of the model. Although there are many parameters and possible choices, the key dynamics boil down to just a handful of parameters – essentially one per exchange flux process.

• The plume fluxes depend principally on the subglacial discharge Q_sg and plume width W_p. For a given subglacial discharge, the main “adjustable” parameter is the plume width. A smaller value will decrease the flux of plume waters that are input to the layer at neutral buoyancy, and this neutral buoyancy layer will tend to be closer to the surface. Conversely, a larger plume width will increase the flux of plume waters, and the neutral buoyancy level will tend to be deeper.

• The only adjustable parameter in the shelf fluxes is the constant C₀, which determines the strength of the fjord–shelf exchange for a given fjord–shelf pressure gradient. Large values of C₀ give quick communication between the fjord and shelf. Small values of C₀ inhibit fjord–shelf exchange, meaning that, for example, plume waters struggle to exit the fjord, and there is recirculation and retention within the fjord. This tends to vertically redistribute waters in the fjord, leading to a fjord–shelf exchange that is weaker and more distributed in the vertical compared to that obtained with a large value of C₀.

• Within the iceberg fluxes, the main adjustable parameter is M₀, which controls the iceberg melt rate for a given temperature. For a given layer temperature, an increase in M₀ will give a higher iceberg melt rate and therefore a higher iceberg melt flux. This has a greater cooling and freshening effect on fjord waters, will drive greater upwelling of fjord waters by icebergs, and will, in turn, drive greater fjord–shelf exchange.

• Lastly, in the vertical mixing, the main adjustable parameter is K₀, which controls the magnitude of vertical mixing of heat and salt. A higher value smoothes out vertical gradients in temperature and salinity, which can then contribute to the vertical structure of the fjord–shelf exchange. There are no parameters associated with the vertical advective fluxes because these are calculated as a residual of all the other fluxes.

In summary, for a given geometry and forcings, the main adjustable parameters controlling the dynamics are the plume width W_p, the shelf exchange parameter C₀, the iceberg melt parameter M₀, and the vertical mixing parameter K₀.

4 Validation against MITgcm

We provide an initial validation of FjordRPM by comparing it with the established general circulation model MITgcm in both idealised (Sect. 4.1) and realistic (Sect. 4.2) cases. In the idealised cases, we consider a simple fjord–shelf geometry and simulate three fundamental modes of glacial fjord circulation: (i) buoyancy-driven circulation resulting from the input of subglacial discharge; (ii) intermediary circulation driven by shelf variability; and (iii) iceberg melt-driven circulation. In the realistic case, we consider a nearly 3-year simulation of Sermilik Fjord, SE Greenland (Sanchez et al., 2024). While MITgcm is not “the truth”, it is a full three-dimensional hydrodynamic model, its dynamics are much more comprehensive than FjordRPM, and it has been widely used in fjord studies (e.g. Xu et al., 2012; Sciascia et al., 2013; Carroll et al., 2015; Cowton et al., 2016; Slater et al., 2018; Zhao et al., 2022; Hager et al., 2022; Sanchez et al., 2024). It therefore provides a benchmark against which to test FjordRPM. We have allowed ourselves to vary a single FjordRPM parameter to improve the comparison to the MITgcm simulations – this is the shelf exchange parameter, which takes the same value $C_{\mathrm{0}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}}$ s across all the FjordRPM simulations shown. All other parameter values used are given in Table A2.

FjordRPM has no awareness of across-fjord variability, and the only awareness of along-fjord variability is that we have an estimate of horizontal velocity at the glacier (Sect. 2.2) and fjord mouth (Sect. 2.3). Thus, the main quantities to compare with FjordRPM are the along- and across-fjord-averaged temperature and salinity from MITgcm and the fjord–shelf fluxes at the fjord mouth from MITgcm.

4.1 Idealised simulations

The idealised MITgcm simulations (Fig. 3) consider a fjord that is 60 km long and uniformly 6 km wide and 800 m deep. At the fjord head, there is a single glacier with grounding line depth 800 m, and at the fjord mouth, there is a sill of depth 400 m. The fjord opens out onto a shelf that is uniformly 600 m deep, 182 km wide (i.e. in the along-fjord direction), and 225 km long (i.e. in the across-fjord direction); note that only part of the shelf is shown in Fig. 3a. The horizontal resolution is 250 m within the fjord, telescoping to around 3 km at the far reaches of the shelf. The vertical resolution is 10 m near the surface, increasing to 50 m at the bottom. Temperature and salinity at the boundaries of the shelf (Fig. 3b and c) are, depending on the simulation, restored to either “summer” or “winter” conditions, with summer conditions coming from a conductivity–temperature–depth (CTD) profile on the continental shelf close to Sermilik Fjord, SE Greenland, and winter conditions coming from an expendable CTD in the mouth of the same fjord (Straneo et al., 2011). The initial conditions throughout the domain are taken to be the same as these far-field boundary conditions, but the simulations are run for long enough that there should be no memory of the initial conditions. The intermediary circulation simulations involve a time-varying wind stress applied to the shelf (Fig. 3d; further details in Sect. 4.1.2), and the iceberg melt-driven simulations require an iceberg concentration (Fig. 3e; further details in Sect. 4.1.3). The simulations take place on an f-plane with $f=\mathrm{1.35}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹. In sum, these simulations are intended to be representative of large glacial fjords in Greenland.

Figure 3 Setup of the idealised MITgcm simulations. (a) Bathymetry of part of the model domain. The red dashed square shows the region from which we take the shelf boundary conditions for FjordRPM. (b) Temperature and (c) salinity to which the MITgcm simulations are restored at the north, south, and east boundaries for summer and winter cases. (d) Wind stress applied to the shelf region in the intermediary circulation simulations. (e) Iceberg ice–ocean surface area concentration profile used in the iceberg simulations.

The equivalent FjordRPM simulations use the same fjord geometry, with 60 layers of 13.3 m thickness and a time step of 0.1 d. The “shelf” profiles T^s(z,t) and S^s(z,t) used to force FjordRPM are extracted from the MITgcm simulations in a 6 km×6 km region on the shelf adjacent to the fjord (Fig. 3a). This is because shelf dynamics mean that conditions in this region can differ significantly from the MITgcm boundary conditions (Fig. 3b and c), particularly for the intermediary circulation simulations, and we wish FjordRPM to experience the same shelf forcing as the fjord in MITgcm does.

4.1.1 Buoyancy-driven circulation

Our first test case is the buoyancy-driven circulation driven by the input of subglacial discharge from beneath tidewater glaciers, thought to be dominant in Greenland's glacial fjords during summer. In MITgcm, the input of subglacial discharge and the resulting plume are represented using the “Iceplume” package (Cowton et al., 2015), and the boundary conditions at the edge of the MITgcm domain are the summer conditions (Fig. 3b and c). In both MITgcm and FjordRPM, the plume width is W_p=500 m, and we consider three values of subglacial discharge: Q_sg=100, 300, or 900 m³ s⁻¹. The subglacial discharge is held constant in time for the full 400 d length of the simulations. This is clearly longer than a melt season, but it enables the simulations to reach a steady state, which provides a good state for comparison. There is no wind stress over the shelf and no icebergs in the fjord. The plots use properties averaged over days 390–400 of the simulations.

The result obtained using MITgcm in the Q_sg=300 m³ s⁻¹ case is shown in Fig. 4. The dynamics are as expected from previous studies, with up-fjord flow at depth due to entrainment into the plume and down-fjord flow of plume waters closer to the surface, with an intensification of these currents over the sill. In temperature and salinity, the sill blocks the deep warm and salty water on the shelf from getting into the fjord. Properties above the sill are relatively similar between the fjord and shelf, though upwelling by the plume results in a slight warm and salty anomaly in the outflowing layer.

Figure 4 MITgcm simulation of steady-state fjord dynamics obtained under sustained subglacial discharge of Q_sg=300 m³ s⁻¹. (a) Along-fjord velocity averaged over the width of the fjord. Positive velocities are directed from fjord to shelf. (b) Temperature and (c) salinity averaged over the width of the fjord. For x-axis values greater than 0, we are out on the shelf and there is no fjord width to average over, but we average over the same range of y values as if we were in the fjord.

A comparison of MITgcm and FjordRPM is shown in Fig. 5. Simulations for all three values of subglacial discharge are shown but are qualitatively similar. FjordRPM captures the profile of the fjord–shelf exchange velocity over the sill very well, including the inflow at depth, the shallower outflow, and the negligible exchange close to the surface (Fig. 5a). Even details in MITgcm, such as the roughly Gaussian velocity profile of the outflow and the intensification of inflow just above the sill, are present in FjordRPM. The match in fjord-average temperature and salinity is also very good (Fig. 5c and d); in particular, FjordRPM captures the divergence in fjord-to-shelf properties below the sill depth, the warm and salty anomaly in the outflowing layer, and the match of fjord-to-shelf properties close to the surface. FjordRPM also skilfully captures the impact of varying subglacial discharge, mirroring the strengthening of the circulation and shallowing of the outflowing layer under increasing subglacial discharge.

Figure 5 Comparison of MITgcm and FjordRPM when simulating subglacial discharge buoyancy-driven circulation. (a) Fjord–shelf exchange velocity (i.e. along-fjord velocity at the sill, averaged across the width of the fjord), calculated as $-Q_j^{Vs}/W{H}_j$. Positive velocities are directed from fjord to shelf. (b) Overturning stream function at the sill. The small nonzero residual at the surface is the net flux due to the input of subglacial discharge into the fjord. (c) Temperature and (d) salinity, averaged over the length and width of the fjord. The grey dashed horizontal line in all panels denotes the sill depth.

Turning to the differences between MITgcm and FjordRPM, which are minor, we note that the outflows in FjordRPM are consistently slightly deeper than in MITgcm and that the intensification of the inflow just above the sill is stronger in MITgcm. This results in a slightly different shape of the overturning stream function (Fig. 5b). Overall, however, FjordRPM captures the dynamics present in MITgcm extremely well.

4.1.2 Intermediary circulation

The second test case is intermediary circulation driven by variability in conditions on the shelf. The MITgcm boundary conditions are the winter profiles from Fig. 3b and c because the classic instance of intermediary circulation is driven by strong winds over the continental shelf and this is more common in winter (e.g. Jackson et al., 2014). To produce these dynamics in MITgcm, we impose a temporally varying southwards wind stress on the “shelf” portion of the domain in which, with a return time of 10 d, there is a wind event that lasts 5 d with a stress that peaks at 1 Pa after 2.5 d (Fig. 3d). After around five such cycles (total 50 d), the fjord in MITgcm reaches a state where the dynamics repeat themselves every 10 d and the initial conditions have been forgotten. There are no icebergs and no subglacial discharge in this simulation. We analyse the period between days 40 and 50 of the simulation.

In MITgcm, for the first 4 d of the 10 d cycle (Fig. 3d), the southward wind deepens shelf isopycnals close to the coast, driving flow into the fjord in the upper layer and out of the fjord at depth (Fig. 6a), associated with negative values of the overturning stream function approaching 150 mSv (Fig. 6c). As the fjord adjusts to the shelf, the isopycnals in the fjord deepen, and the layer of cold water at the surface thickens (Fig. 6c and e). Once the wind relaxes, after 4 d of the cycle, the circulation reverses, and the fjord rebounds to the pre-wind state (Fig. 6a, c, and e).

Figure 6 Comparison of MITgcm and FjordRPM when simulating the fjord response to wind-driven shelf variability. Panels (a), (c), and (e) show the MITgcm results, and panels (b), (d), and (f) show the FjordRPM results. Note that only the depth range above the sill is shown. (a, b) Hovmöller plots of exchange velocity over the sill, averaged across the width of the fjord. Positive velocities are directed out of the fjord, and the contours show potential density. (c, d) Equivalent plots of the overturning stream function over the sill. (e, f) Hovmöller plots of temperature averaged across the length and width of the fjord.

When forced by the time-varying shelf properties taken from the fjord mouth in the MITgcm simulations (red box in Fig. 3a), FjordRPM captures the same key features of the intermediary circulation, including the direction, timing, and vertical structure of the circulation (Fig. 6b, d, and f). There are differences here, though, with the FjordRPM fjord–shelf exchange generally more sluggish and less surface-intensified than in MITgcm (Fig. 6a vs. b), giving a weaker overturning stream function (Fig. 6c vs. d) and less pronounced deepening then shallowing of isopycnals and isothermals (Fig. 6). The inflowing and outflowing layers are consistently separated by the 26.8 kg m⁻³ potential density contour (Fig. 6a and b); using this to separate the flow into upper and lower layers and then plotting the upper layer volume flux show that FjordRPM captures the timing of the circulation well but underpredicts the strength of the circulation by around 30 % compared to MITgcm (Fig. 7a). As a consequence, the mean temperature of the fjord through the wind event varies less in FjordRPM than in MITgcm (Fig. 7b).

Figure 7 Further comparison of MITgcm and FjordRPM in simulating intermediary circulation. (a) Volume flux between the 26.8 kg m⁻³ potential density contour (Fig. 6a, b) and the surface. Positive values are directed into the fjord. (b) Fjord temperature, averaged over the fjord length, width, and depth, through the 10 d wind cycle.

Improving the match between MITgcm and FjordRPM for these dynamics is not easy. The obvious parameter with which to control the fjord–shelf exchange in FjordRPM is C₀. Neglecting the feedback of C₀ on the pressure gradient, increasing C₀ will strengthen the fjord–shelf exchange (Eq. 16). In this case, however, we found that increasing C₀ does not significantly strengthen the fjord–shelf exchange. We hypothesise that this is because the exchange is, or becomes, large enough to dampen the fjord–shelf pressure gradients that drive the flow. We can therefore say that FjordRPM captures the intermediary circulation relatively well but appears slightly sluggish compared to MITgcm; we return to this point in the discussion.

4.1.3 Iceberg melt-driven circulation

Lastly, we compare the ability of MITgcm and FjordRPM to capture circulation driven by iceberg melting. In MITgcm, the icebergs are represented by the “Iceberg” package (Davison et al., 2020), in which the icebergs are thermodynamically active but not mechanically active. That is, the simulation accounts for the cooling and freshening effect of iceberg meltwater but not for the drag exerted on the circulation by the presence of icebergs. MITgcm is configured to have the same iceberg melt parameterisation as FjordRPM (Eq. 20). The distribution of icebergs in the fjord is set by the iceberg surface area concentration profile $\mathit{\eta }\left(z\right)=\frac{\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}}{WL}\exp (-|z|/\mathrm{100})$, where L=60 km is the fjord length and W=6 km is the fjord width (Fig. 3e). The surface area of icebergs in an MITgcm grid cell with dimensions $(\mathit{\delta }x,\mathit{\delta }y,\mathit{\delta }z)$ is then η(z)δx δy δz, with units m². The iceberg concentration in MITgcm is set to be uniform in the horizontal. In FjordRPM, the surface area of icebergs in layer j is I_j=η(z_j)WLH_j, where z_j is the depth of layer j. The total submerged surface area of icebergs in this setup is 200 km², consistent with observational estimates for a large Greenland fjord such as Sermilik (Enderlin et al., 2016). The simulations are run with both the summer and winter boundary conditions from Fig. 3b and c. There is no subglacial discharge and no wind stress on the shelf.

The obtained iceberg melt rates and fluxes in MITgcm and FjordRPM are shown in Fig. 8. The melt rates and fluxes are comparable with observations from Sermilik Fjord (Enderlin et al., 2016); this is no coincidence, as it is how we set our suggested value of the iceberg melt parameter, $M_{\mathrm{0}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}$ $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}(\mathrm{°}\mathrm{C}{)}^{-\mathrm{1}}$ (Table A2). MITgcm and FjordRPM show mostly close agreement in both melt rate and melt flux. Given that the two methods of simulation have been configured with the same iceberg distribution and the same melt parameterisation (which depends only on thermal forcing), this is not surprising and is only really a test of the modelled fjord temperature.

Figure 8 Comparison of MITgcm- and FjordRPM-simulated (a) iceberg melt rate and (b) iceberg melt flux in the simulation of fjord dynamics driven by iceberg melt.

A tougher test of FjordRPM comes from comparing the dynamics induced by the melting icebergs. In MITgcm with the summer shelf boundary conditions, the dominant feature is an overturning cell at intermediate depth with an outflow centred at 200 m (Figs. 9a and 10a and b). We interpret that this overturning cell is driven by deep iceberg melt (Fig. 8) that is able to upwell through the relatively unstratified waters at depth before ceasing to rise and flowing horizontally below the more stratified upper 200 m (Figs. 9c and 10d). The upwelling entrains ambient fjord waters, setting up the circulation cell, and gives a small warm anomaly in the fjord relative to the shelf at 200 m depth (Fig. 10c). The presence of iceberg melt and upwelling below sill depth (Fig. 8) leads to some of the below-sill water being replaced; hence, the fjord becomes slightly cooler and fresher than the shelf just below sill depth (Fig. 10c and d). Finally, the summer MITgcm simulation shows a strong cooling of the fjord relative to the shelf at the surface (Fig. 10c).

Figure 9 MITgcm simulation of steady-state fjord dynamics driven by iceberg melting with summer shelf boundary conditions. (a) Along-fjord velocity, averaged over the width of the fjord. Positive velocities are directed from fjord to shelf. (b) Temperature and (c) salinity, averaged over the width of the fjord.

Figure 10 Comparison of MITgcm and FjordRPM when simulating circulation driven by iceberg melting. (a) Fjord–shelf exchange velocity. (b) Overturning stream function at the sill. (c) Temperature and (d) salinity, averaged over the length and width of the fjord. The grey dashed horizontal line in all panels denotes the sill depth.

In MITgcm with the winter shelf boundary conditions (the winter equivalent of Fig. 9 is not shown), the cooler surface waters mean that iceberg melt rates are lower than in summer (Fig. 8). The overturning stream function has a second maximum located at 80 m depth (Fig. 10b) that we interpret is driven by iceberg upwelling in the less-stratified near-surface waters that are present in winter (Fig. 10d).

Overall, FjordRPM does very well at representing both the summer and winter dynamics (Fig. 10), capturing (i) the main and secondary circulation cells, (ii) the strength of these cells, (iii) the warm anomalies in the fjord at intermediate depths, (iv) the cold anomaly in the fjord at the surface in summer, and (v) the below-sill cooling and freshening. Considering how these dynamics come about, in MITgcm, cells with icebergs are freshened and cooled, and any upwelling results from a solution of the full hydrodynamic equations (note, however, that these plumes of upwelling iceberg meltwater are not properly resolved in a fjord-scale MITgcm simulation). In FjordRPM, upwelling is simply parameterised as described in Sect. 2.4. Given these very different approaches, it is encouraging that FjordRPM matches the MITgcm simulations so well.

4.2 Realistic simulation of Sermilik Fjord

As a final test of FjordRPM, we compare to a nearly 3-year, realistic MITgcm simulation of Sermilik Fjord and the adjacent shelf in SE Greenland (Fig. 11), conducted by Sanchez et al. (2024). In this context, “realistic” means that their simulation used the actual bathymetry of the fjord and shelf and was forced at the glacier by subglacial discharge from a regional climate model (RACMO; Noël et al., 2018) at the surface by an atmospheric reanalysis (ERA5; Hersbach et al., 2020) and at the ocean boundary by an ocean reanalysis (ASTE; Nguyen et al., 2021). Together with observational validation from this well-studied fjord, the simulations in Sanchez et al. (2024) are as close to realistic as we can currently get (aside from the presence of icebergs, which they did not consider). Full details may be found in their paper.

Figure 11 Sermilik Fjord, SE Greenland, as modelled by Sanchez et al. (2024). Note that their model domain included a larger region of the shelf than is indicated by this plot. The white annotated patches show the locations of the three main glaciers that discharge into the fjord. The red dotted line shows the fjord mouth flux gate, and the yellow dashed square shows the region from which we extract properties from MITgcm to force FjordRPM. The inset shows the location in SE Greenland.

We attempt to undertake the same simulation in FjordRPM, making the setup as close as possible given the much-simplified nature of the reduced-physics model. Whereas MITgcm used realistic bathymetry (Fig. 11), FjordRPM has no sill and comprises a cuboid fjord of width 6 km, depth 900 m, and length 101 km (ensuring that the total volume of the fjord is the same in both MITgcm and FjordRPM). Whereas in MITgcm there are three glaciers with plumes in side-arms of the main fjord, FjordRPM has three plumes in the cuboid fjord with respective initiation depths of 650, 160, and 465 m, corresponding to the grounding line depths of the three glaciers. The plume width of 280 m in MITgcm is inherited from the MITgcm model resolution and is set to be the same in FjordRPM, and the subglacial discharge forcing is identical. Unlike MITgcm, FjordRPM has no surface forcing. For the shelf forcing in FjordRPM, we extract properties from the MITgcm simulation close to the fjord mouth (Fig. 11). We configure FjordRPM with 60 layers of 15 m thickness, a time step of 0.1 d, and parameters as in Table A2. The FjordRPM simulation runs in a few seconds on a laptop.

We compare MITgcm and FjordRPM in terms of simulated fjord properties (Fig. 12) and fjord–shelf exchange fluxes (Fig. 13). The fjord properties extracted from MITgcm are averaged over all model points up-fjord of a flux gate at the fjord mouth (Fig. 11), while from FjordRPM, they are simply the layer properties. The exchange fluxes from MITgcm are those passing through the flux gate at the fjord mouth, while from FjordRPM, they are the fjord–shelf exchange fluxes.

Figure 12 Sermilik Fjord properties from March 2015 to December 2017, as modelled by MITgcm (blue, Sanchez et al., 2024) and by FjordRPM when forced by shelf properties from MITgcm (red). (a–c) Temperature averaged over the indicated depths; (d–f) salinity averaged over the same depths. The plots also show the mean absolute error (MAE) of FjordRPM compared to MITgcm.

When compared to MITgcm, FjordRPM captures the evolution of fjord properties relatively well (Fig. 12). It captures the timing and magnitude of the seasonal cycle in both temperature and salinity at all depths. It also recreates the high-frequency variability in properties that arises during winter storm systems (e.g. Fig. 12e). There are differences between MITgcm and FjordRPM – in general, FjordRPM appears warmer and saltier than MITgcm. There are also specific depths and times where the properties diverge, for example, during summer at 50–100 m depth (Fig. 12a), where FjordRPM becomes warmer than MITgcm, likely due to differences in plume neutral buoyancy. We have not dug further into these differences, as to do so properly would require more space than feels appropriate for the present paper, but possible reasons for the differences include the specifics of how the shelf forcing is sampled (Fig. 11) and the lack of realistic hypsometry in FjordRPM.

FjordRPM also does well at capturing fjord-to-shelf exchange at approximately monthly or longer timescales when compared to MITgcm (Fig. 13, thicker lines). FjordRPM captures the seasonal cycle in fjord-to-shelf volume exchange (Fig. 13a), with increased exchange during winter and more quiescent periods during summer. Considering freshwater fluxes, FjordRPM captures the variable exchange of freshwater during winter and the emergence of a net freshwater flux of 1–5 mSv from the fjord to the shelf during summer (Fig. 13b). In fjord-to-shelf heat exchange, FjordRPM has much of the same variability as MITgcm, but in terms of absolute value, FjordRPM exceeds MITgcm during summer, and MITgcm exceeds FjordRPM during winter (Fig. 13c). We suspect that these differences may result from air–sea heat fluxes, which are included in MITgcm but not in FjordRPM – indeed, excluding the surface 30 m from the exchange flux changes the sign of the heat flux in MITgcm during summer, leading to better agreement between the models (Fig. 13d).

Figure 13 Fluxes at the fjord mouth (section indicated in Fig. 11) as modelled by MITgcm (blue, Sanchez et al., 2024) and by FjordRPM when forced by shelf properties from MITgcm (red). (a) Volume exchange obtained by integrating the absolute value of the meridional velocity over the fjord mouth cross-section (i.e. integrating over the fjord width and depth). (b) Freshwater exchange assuming S_ref=34.9 g kg⁻¹, which is a typical value at the bottom of the fjord. Positive values indicate a net flux of freshwater from the shelf into the fjord. (c) Heat exchange. Positive values indicate a net flux of heat from the shelf into the fjord. (d) Same as for panel (c) but excluding the surface 30 m from the integration. In all plots, the thin lines are daily values, while the thick lines are the result of smoothing the daily values with a 30 d centred moving window. Also included is the MAE of FjordRPM compared to MITgcm, with the normal-font weight number calculated on the daily values and the bold number calculated on the smoothed values.

The greatest differences in exchange fluxes between FjordRPM and MITgcm are at weekly or shorter timescales (Fig. 13, thinner lines), where MITgcm shows a greater magnitude of variability, especially during summer. This is reflected in the mean absolute errors of FjordRPM compared to MITgcm, which are significantly larger when incorporating shorter timescales (Fig. 13). We have forced FjordRPM with daily output from the MITgcm simulation – it is possible that FjordRPM would show greater variability at short timescales if we extracted, say, hourly output from MITgcm. However, this “sluggish” nature of FjordRPM appears similar to the idealised intermediary circulation simulations above and hence may represent something more fundamental about the implemented physics. We expand on this in the discussion and conclude at this point that FjordRPM shows good promise at capturing fjord properties and longer-term fjord–shelf exchange fluxes when applied to real fjords.

5 Discussion

Having motivated, designed, coded, and tested FjordRPM, we now discuss its strengths and weaknesses together with potential future work and applications. One aspect that has not been discussed so far is the efficiency of the model. Consisting of around 400 lines of code, the model takes around a second to run a fjord for one model year on a laptop. FjordRPM is clearly very different in character to a general circulation model and is not suited to all purposes, but this compares to hundreds or thousands of CPU hours to run equivalent simulations (such as the validation experiments in this paper) using a general circulation model. The principal factors determining the FjordRPM run time are the number of layers (i.e. the vertical resolution), the time step, and how frequently the subglacial discharge plume dynamics are updated.

When compared to several idealised and a single realistic MITgcm simulation, FjordRPM shows excellent ability to capture fjord properties and circulation. Specifically, FjordRPM shows the ability to capture (i) the influence of a sill on below-sill fjord temperature and salinity (Fig. 5); (ii) the influence of subglacial discharge-driven upwelling on fjord properties in the upper water column (Fig. 5); (iii) the influence of synoptic (Figs. 6, 7, and 12) and seasonal (Fig. 12) shelf variability on fjord properties; and (iv) the role of icebergs in cooling, freshening, and upwelling of fjord properties (Fig. 10). Furthermore, FjordRPM captures very well the exchange between fjord and shelf – in both a depth-resolved and depth-integrated manner – induced by (i) subglacial discharge (Fig. 5), (ii) shelf variability (Figs. 6, 7, and 13), and (iii) iceberg melt (Fig. 10). The listed dynamics comprise many of the processes understood to be important in glacial fjords. There are also aspects in the comparison to MITgcm in which FjordRPM is less successful. FjordRPM is generally sluggish compared to MITgcm in its response to synoptic shelf variability, as seen in both the idealised intermediary circulation simulations (Figs. 6 and 7) and the realistic simulation of Sermilik Fjord (Fig. 13). From the implemented physics, we might expect FjordRPM to be more sluggish in this manner – the fjord–shelf exchange fluxes (Sect. 2.3) are based on a steady-state balance that does not consider an acceleration term. In addition, FjordRPM has no along-fjord variability, so waters that come into the fjord from the shelf are instantly mixed along the full fjord, limiting the fjord–shelf density gradient. In contrast, MITgcm has along-fjord resolution, which may allow it to have high fjord–shelf density gradients purely in the region of the fjord mouth, leading to higher-frequency and higher-magnitude fjord–shelf exchange.

A related issue is how exactly to derive the shelf boundary conditions required by FjordRPM. In the validation against MITgcm, we have sampled these boundary conditions from a region on the shelf very close to the fjord mouth (Figs. 3 and 11). However, there can be spatial variability in water properties, particularly in the presence of variable bathymetry (e.g. Fig. 11). We found that adjusting the definition of the shelf sampling region, say by making the region smaller or larger or moving it around, changes the FjordRPM results quantitatively but not qualitatively. That is, the plume neutral buoyancy level may adjust, the fjord–shelf exchange fluxes may increase or decrease, and the fjord water properties may warm or cool, but the magnitude of these changes is small compared to the main features and variability in the FjordRPM solution.

A broader point, however, is whether it is necessary to have a well-resolved simulation of the continental shelf to provide the shelf boundary conditions for FjordRPM. The comparisons presented in this paper have shown very good agreement between FjordRPM and MITgcm in the fjord when FjordRPM is forced with water properties extracted from MITgcm just outside the fjord mouth, but as illustrated in the MITgcm simulation of Sermilik Fjord (Sanchez et al., 2024), these fjord mouth properties may themselves be substantially modified during cross-shelf transit in ways that are not effectively captured by relatively coarse ocean reanalyses or projections. Future work could seek to quantify to what extent this is a significant issue, that is, what difference it makes to the FjordRPM solution when we force it with a coarse vs. resolved representation of the shelf outside the fjord. If there is a significant difference, a further link in the model chain may then be required, either using higher-resolution regional simulations (Verjans et al., 2023) or simplified models of shelf processes more akin to the philosophy of FjordRPM.

There are similarly many avenues for further development of the model. One obvious aspect is to implement realistic fjord hypsometry. At present, the fjord is assumed to be a cuboid that we can adjust to ensure it has the same overall volume as a real-world fjord (e.g. Sect. 4.2), but we cannot ensure it has the right volume at all depths. Another prominent aspect for improvement is the fjord–shelf exchange parameterisation (Sect. 2.3). Previous work has suggested that, depending on the geometry, stratification, circulation, and shelf variability, fjord–shelf exchange can be dictated by a form of geostrophic balance, hydraulic control, or friction and mixing (Zhao et al., 2021; Nilsson et al., 2023; Sanchez et al., 2023; Bonneau et al., 2024). Jackson et al. (2018) suggest that, because Greenlandic fjord widths are comparable to the deformation radius, three-dimensional dynamics play an important role and that, furthermore, the timescale of shelf variability relative to the adjustment time of the fjord is an important factor. We have found that FjordRPM appears very effective at capturing a range of circulations with a very simple fjord–shelf exchange parameterisation (Sect. 2.3), but its somewhat empirical nature means that it may need to be tuned individually for each fjord. A more sophisticated fjord–shelf exchange parameterisation and a further set of validation experiments could seek to capture various fjord–shelf dynamical balances in one framework and allow a “universal” tuning of this exchange parameterisation.

There are also fjord processes that we have not represented at all but which may be important in particular systems or for particular research questions. These include ambient (outside-of-plume) melting of the glacier front, fluxes of heat and freshwater at the fjord surface, sea ice, wind stress, and tides. Iceberg concentration, currently prescribed by a static ice–ocean surface area, could be made to evolve upon the balance of input from calving, loss to melting, and export to the shelf.

Another avenue for future work concerns validation of the model beyond the initial tests that we have undertaken here. This could take the form of further MITgcm simulations that are targeted to test specific aspects of the exchange parameterisations, but it would be particularly valuable to test FjordRPM against observations, either looking at the model's ability to capture fjord-to-fjord differences (say, by using CTD profiles from many of Greenland's fjords) or its ability to capture variability in depth and time at an individual fjord (say, by using mooring data).

Lastly, we return to the potential applications of FjordRPM. While being mindful of the issues just discussed, the ability to capture fjord temperature under a range of forcings (Figs. 5, 7, 10, and 12) means it could be used to generate ocean boundary conditions for ice sheet models. Similarly, the ability to capture freshwater export to the shelf (Figs. 5, 7, 10, and 13) suggests it could be used to generate freshwater boundary conditions for ocean models. In both directions, it offers a means of coupling ice sheet and ocean models in a practical and efficient manner. Due to the simplicity and transparent nature of the dynamics (e.g. the ability to easily turn on or off, or tweak, individual processes), FjordRPM may also have a role in investigating fjord processes.

6 Conclusions

We have presented the design and initial validation of FjordRPM – a reduced-physics model for glacial fjords. The model consists of a number of vertically stacked layers, each of which extends over the full width and length of the fjord at a given depth. Exchanges of volume, heat, and salt into and out of layers are parameterised to represent important fjord processes, including subglacial discharge-driven upwelling and submarine melting, exchange with the shelf, iceberg-driven upwelling and melting, vertical mixing, and vertical advection. The model thus lies somewhere between a box model and a “1.5-dimensional” model (i.e. a one-dimensional model with lateral exchanges). The model has the capability to represent a sill and multiple glaciers in a single fjord.

We have validated FjordRPM by comparing to idealised and realistic simulations in the established general circulation model MITgcm, finding that FjordRPM is successful in capturing almost all tested aspects of glacial fjord properties and circulation but is sluggish in response to high-frequency shelf variability. FjordRPM is also highly efficient, taking around a second per fjord per model year on a laptop, in contrast to the hundreds or thousands of CPU hours required for MITgcm to run the equivalent simulations. While further development and validation would be valuable, we present FjordRPM as a promising tool for conceptualising fjord dynamics and for linking ice masses and the ocean in ice and climate projection efforts.

Appendix A: Scalings for iceberg-driven upwelling

Within an unstratified body of water (e.g. within each FjordRPM layer), Magorrian and Wells (2016) show in their Eq. (31) that, up to a constant, the upwelling velocity scales as

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\displaystyle }u\sim {\left({\displaystyle \frac{\mathrm{\Delta }{\mathit{\rho }}_u}{{\mathit{\rho }}_I}}g\sin \mathit{\varphi }X\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

in which Δρ_u is the density difference between the upwelling plume and the ambient, ρ_I is a reference density, g is gravity, ϕ is the angle of the interface measured from the horizontal, and X is the along-ice distance. To derive a simple scaling for FjordRPM, we take X to be the layer thickness H_j and assume the iceberg sides are vertical so that sin ϕ=1. Continuing to follow Magorrian and Wells (2016), the density difference is given by

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }{\mathit{\rho }}_u={\displaystyle \frac{M_u}{{\mathit{\alpha }}_i+M_u}}\mathrm{\Delta }{\mathit{\rho }}_{I\mathrm{0}}^{\text{ef}},\end{array}$$

where α_i is the entrainment coefficient as in Sect. 2.4, $M_u=\dot{m}/u$ is the constant that relates velocity and melt rate, and $\mathrm{\Delta }{\mathit{\rho }}_{I\mathrm{0}}^{\text{ef}}$ is the density difference between meltwater and the ambient. Assuming M_u≪α_i and substituting Eq. (A2) into Eq. (A1), we get

$$\begin{array}{}\text{(A3)}& {\displaystyle }{\displaystyle }u\sim {\left({\displaystyle \frac{g\mathrm{\Delta }{\mathit{\rho }}_{I\mathrm{0}}^{\text{ef}}}{{\mathit{\rho }}_I}}{\displaystyle \frac{\dot{m}}{{\mathit{\alpha }}_{i}u}}{H}_j\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

Now, we note that $g\mathrm{\Delta }{\mathit{\rho }}_{I\mathrm{0}}^{\text{ef}}/{\mathit{\rho }}_I$ is $g_{j,\text{melt}}^{\prime }$ as given by Eq. (22) and that $\dot{m}$ can be written as $Q_j^{\text{melt}}/I_j$ (i.e. the melt rate is the melt flux divided by the surface area undergoing melting). We can then solve for u to obtain the estimate of the upwelling velocity that is used in the model (Eq. 21):

$$\begin{array}{}\text{(A4)}& {\displaystyle }{\displaystyle }u\equiv v_j\sim {\left[{\displaystyle \frac{Q_j^{\text{melt}}{g}_{j,\text{melt}}^{\prime }{H}_j}{{\mathit{\alpha }}_i{I}_j}}\right]}^{\mathrm{1}/\mathrm{3}}.\end{array}$$

To obtain the length scale that controls how far the upwelling rises, start from Eq. (37) of Magorrian and Wells (2016):

$$\begin{array}{}\text{(A5)}& {\displaystyle }{\displaystyle }{l}_{\mathit{\rho }}=\mathrm{\Delta }{\mathit{\rho }}_u{\left({\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{a}}}{\partial X}}\right)}^{-\mathrm{1}}.\end{array}$$

Using the same substitutions as for the upwelling velocity, this may be written as (again under M_u≪α_i)

$$\begin{array}{ll}{\displaystyle }{l}_{\mathit{\rho }}& {\displaystyle }={\displaystyle \frac{M_u}{{\mathit{\alpha }}_i+M_u}}\mathrm{\Delta }{\mathit{\rho }}_{I\mathrm{0}}^{\text{ef}}{\left({\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{a}}}{\partial X}}\right)}^{-\mathrm{1}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\dot{m}}{{\mathit{\alpha }}_{i}u}}{\displaystyle \frac{g_{j,\text{melt}}^{\prime }{\mathit{\rho }}_I}{g}}{\left({\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{a}}}{\partial X}}\right)}^{-\mathrm{1}}\\ \text{(A6)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{Q_j^{\text{melt}}{g}_{j,\text{melt}}^{\prime }}{{\mathit{\alpha }}_i{I}_{j}u}}{\left({\displaystyle \frac{\partial g^{\prime }}{\partial X}}\right)}^{-\mathrm{1}},\end{array}$$

where, in the last step, we have used the fact that $g^{\prime }=g\mathit{\delta }\mathit{\rho }/{\mathit{\rho }}_{\text{ref}}$ (Eq. 1). We can then use Eq. (A4) to rewrite this as

$$\begin{array}{}\text{(A7)}& {\displaystyle }{\displaystyle }{l}_{\mathit{\rho }}={\displaystyle \frac{v_j^{\mathrm{2}}}{H_j}}{\left({\displaystyle \frac{\partial g^{\prime }}{\partial X}}\right)}^{-\mathrm{1}}.\end{array}$$

If we finally estimate the derivative for layer j using the first-order finite difference

$$\begin{array}{}\text{(A8)}& {\displaystyle }{\displaystyle \frac{\partial g^{\prime }}{\partial X}}\approx {\displaystyle \frac{g_{j+\mathrm{1},j}^{\prime }}{\frac{\mathrm{1}}{\mathrm{2}}\left(H_j+H_{j+\mathrm{1}}\right)}},\end{array}$$

then we get the final estimate of the upwelling length scale that is used in the model (Eq. 24):

$$\begin{array}{}\text{(A9)}& {\displaystyle }{\displaystyle }{l}_{j+\mathrm{1},j}^{\text{ice}}={\displaystyle \frac{v_{j+\mathrm{1}}^{\mathrm{2}}}{H_{j+\mathrm{1}}}}{\displaystyle \frac{H_j+H_{j+\mathrm{1}}}{\mathrm{2}{g}_{j+\mathrm{1},j}^{\prime }}}.\end{array}$$

Table A1 List of FjordRPM model variables. The first section contains basic variables, the second section contains plume flux variables, the third section is shelf flux variables, the fourth section is iceberg variables, the fifth section is vertical mixing variables, and the last section concerns conservation of layer volume. Starred variables are key inputs or outputs; all others are essentially working variables.

Jenkins (2011)Cowton et al. (2015)Large et al. (1994)Jackson et al. (2017)Large et al. (1994)Large et al. (1994)Enderlin et al. (2016)

Table A2 List of parameters in FjordRPM with suggested values. The first section of the table is intended to be physical constants that should rarely need changing, while the second section contains parameters that should be seen as adjustable, for example, in tuning FjordRPM to MITgcm or observations.