Modeling saline-fluid flow through subglacial channels

Subglacial hydrological systems impact ice dynamics, biological environments, and sediment transport. Previous numerical models of channelized subglacial flow have focused on freshwater in temperate ice without considering variable fluid chemistry and properties. Saline fluids can exist in cold glacier systems where freshwater cannot, making the routing of these fluids critical for understanding their influence on geochemical and physical processes in relevant glacial environments. This study advances previous efforts by modeling saline fluid in cold glacier systems, where variable fluid chemistry significantly influences melt rates and drainage processes. We model the drainage of a hypersaline subglacial lake through an ice-walled channel, highlighting the impact of salinity on channel evolution. The model results show that, in subglacial systems at salinity-dependent melting points, channel walls grow more slowly when fluids have higher salt concentrations, leading to significantly lower discharge rates. At higher salinities, more energy is required to warm the fluid to the new melting point as the brine is diluted, which reduces the energy available for melting the channel walls. We also highlight the impact of increased fluid density on subglacial drainage and the importance of accounting for accurate suspended sediment concentrations when modeling outburst floods. This model provides a framework to assess the impact of fluid chemistry and properties on the spatial and temporal variations of fluid flux.

1 Introduction

Subglacial hydrology is of fundamental importance to the dynamics and evolution of ice masses (Flowers, 2018; Morlighem et al., 2014). The presence, distribution, and geometry of the subglacial water system have direct effects on rates of ice sliding, ice mass flux, erosion, and deposition (e.g., Russell et al., 2006; Bell et al., 2007; Stearns et al., 2008; Siegfried et al., 2016; Larsen and Lamb, 2016; Seroussi et al., 2017; Carrivick and Tweed, 2019; Keisling et al., 2020). The subglacial hydrological system affects the distribution and character of subglacial biological communities and influences water and nutrient flux into surrounding waterbodies (e.g., Neal, 2007; Kjeldsen et al., 2014; Meerhoff et al., 2019; Mikucki et al., 2004; Vick-Majors et al., 2020). Whether the goal is to project future sea level rise, understand glacier bed forms, model ocean circulation, or investigate potential extra-planetary habitats through Earth analogs, we require an understanding of the distribution and dynamics of subglacial hydrological systems (e.g., Nienow et al., 2017; Forte et al., 2016).

Subglacial lakes have been observed to drain episodically through outburst floods and less catastrophically through longer-lived drainage events. There is a significant body of work on modeling the drainage of glacial lakes (e.g., Röthlisberger, 1972; Nye, 1976; Spring and Hutter, 1981; Fowler, 1999; Clarke, 2003; Evatt et al., 2006; Kingslake, 2015; Schoof, 2020; Jenson et al., 2022). Many subglacial hydrology models have assumed that drainage from a subglacial lake occurs through ice-walled channels at the ice–bed interface (e.g., Nye, 1976; Fowler, 1999; Clarke, 2003; Evatt et al., 2006). Collectively, this work has focused on freshwater at the pressure melting point and neglected the consideration of water chemistry.

The chemistry of the subglacial water influences the character of the hydrological system. Depression of the pressure melting point through an increased solute concentration is one potential mechanism to explain the presence of subglacial water in locations with a subglacial temperature below, sometimes significantly below, the pressure melting point (Mikucki et al., 2015). Locations with observable saline discharge occur in both Antarctic and Arctic settings, such as at Blood Falls, Taylor Glacier, East Antarctica, and at Borup Fiord Pass Glacier, Ellesmere Island, Canada (Trivedi et al., 2018; Lyons et al., 2019). The salinity of the englacial brine feeding Blood Falls is approximately 125 psu (compared to ≈ 35 psu for seawater), but the precise geometry of the subglacial brine system beneath Taylor Glacier is not fully understood (Badgeley et al., 2017; Lyons et al., 2019). Hubbard et al. (2004) inferred that a zone 3–6 km upglacier from the terminus contained saturated sediments or ponded water based on radar data, and widespread hypersaline groundwater has been detected as far as 5.7 km upglacier from the terminus using transient electromagnetic techniques (Mikucki et al., 2015). Hypersaline lakes have been inferred to exist beneath the Devon Ice Cap in the Canadian High Arctic using airborne radio echo sounding data, with salinity predicted to be in the range of 140 to 160 psu (Rutishauser et al., 2018, 2022), although a more recent study by Killingbeck et al. (2024) using seismic and electromagnetic data argues that the bed from which the subglacial lakes were inferred is dry or frozen. Despite locations of known saline-fluid flow in subglacial environments, the effects of increased salinity on the geometry and flow in subglacial hydrological systems remain unknown.

Subglacial water composition is expected to impact the geometry and dynamics of the subglacial hydrological system. For instance, the hydraulic potential field is modified through the density of the fluid; saline fluid can have a significantly different flow path compared to freshwater for the same glacier geometry (Badgeley et al., 2017; Rutishauser et al., 2022). Channel size and evolution are also expected to differ as the result of fluid chemistry.

An understanding of the impact fluid chemistry has on subglacial systems is important for mapping and classifying subglacial hydrological features using radar. The size, continuity, and electrical conductivity of subglacial channels determine the detectability of subglacial features by radar remote sensing. Constraints provided by modeling inform radar system design decisions such as power requirements, center frequency, and antenna geometry (Scanlan et al., 2022). The ongoing development of multi-polarization radar systems and radar processing algorithms increases the detectability of variations in subglacial hydrological organization (Scanlan et al., 2022, 2020). The expected geometry of subglacial features is an important specification for the design of radar systems (Scanlan et al., 2022). The response in the geometry of subglacial features to changes in the discharge, position along a flowline, and aqueous chemistry provides constraints for the technological and scientific development of new radar systems.

Basal thermal regimes have been shown to impact the solutes, nutrients, and microbes found in the subglacial systems (Dubnick et al., 2020), and subglacial fluid flow can transport these materials, leading to a change in the geomicrobiology of local and nearby environments (Mikucki et al., 2004). We hypothesize that both the basal thermal regime and solute concentrations influence the subglacial hydrological system by altering the effective pressure and fluid flux (which, in turn, will influence geomicrobiology). A better understanding of the flow dynamics in cold ice is important for characterizing the distinct biogeochemistry in saline subglacial systems.

By modeling saline fluid flow through cold ice, we seek to address the following questions: how significant is the effect of salinity on channel wall melt rates? How does the salt concentration change along the channel in response to the melting of the channel walls? And in what systems is the consideration of fluid chemistry and fluid properties important for understanding subglacial hydrology? To answer these questions, we mathematically investigate channel evolution and discharge over time in response to variable fluid chemistry. The results show that the radius of an evolved channel is smaller for saline fluid than the freshwater equivalent when both glacier and lake systems are at their respective salinity- and pressure-dependent melting points. The smaller channel cross-sections affect the temporal and spatial evolution of fluid flux for saline fluid. The differences related to fluid chemistry are greatest for high discharge rates, which are generated by high-volume lakes and channels with steep bed slopes and circular geometry. Additionally, we find that fluid density has a substantial influence on discharge rates, which has relevance to suspended sediments and outburst floods.

Table 1 List of model parameters and variables. Values of constants are specified in brackets.

2 Model description

We construct a lake drainage model in which the water flows from a subglacial lake through an R channel (Röthlisberger, 1972; Nye, 1976). In contrast to previous approaches, we allow for varying salt concentrations in fluid flowing from a subglacial lake. We follow the implementation of Fowler (1999) and Kingslake and Ng (2013). In particular, the equations describing channel evolution and the conservation of mass, momentum, and energy are identical to those in Fowler (1999), with the differences being (i) the fluid density and the melting point of the fluid are functions of salinity, (ii) the assumptions around fluid temperature are related to salinity, and (iii) an additional equation is needed to solve for salinity along the channel and in time. The model equations from Fowler (1999) are described with our assumptions and in our notation in Sect. 2.1, and changes to the model equations are described in Sect. 2.2. For a list of model variables and parameters, along with the consistently used parameter values, see Table 1.

Figure 1 Schematic of simple glacier geometry and subglacial hydrological system with a R channel draining a subglacial lake.

2.1 Model equations

We assume a subglacial conduit on an inclined bed slope beneath ice of constant thickness (Fig. 1). The negative basic hydraulic gradient is the sum of the terms related to glacier geometry:

$$\begin{array}{}\text{(1)}& \mathit{\psi }={\mathit{\rho }}_{\mathrm{b}}g\sin B-{\displaystyle \frac{\partial P_{\mathrm{i}}}{\partial s}},\end{array}$$

where B is the conduit slope (assumed to be constant along the channel and the same as the bed slope), s is the along-flow coordinate parallel to the bed, ρ_b is the density of the fluid, and g is gravitational acceleration (Kingslake and Ng, 2013, Eq. 5). The ice-overburden pressure P_i (in Pa) is given by P_i=ρ_igH, where H is the glacier thickness. Since we assume that the ice thickness is constant, the change in ice-overburden pressure along the channel is zero, and ψ=ρ_bgsin B.

Assuming an existing channel, the channel walls open due to melt and close due to creep closure. Together, these govern the rate of change of the conduit cross-sectional area S with respect to time t:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial S}{\partial t}}={\displaystyle \frac{m}{{\mathit{\rho }}_{\mathrm{i}}}}-KS{N}^n,\end{array}$$

where m is the melt rate (in $\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) (Fowler, 1999, Eq. 2.1). $K=\mathrm{2}A/\left(n^n\right)$ is a function of the Glen's flow law parameter A and exponent n (Evatt et al., 2006). The effective pressure N is the difference between the ice-overburden pressure and the water pressure in the channel. We calculate A as a function of ice temperature using the Arrhenius relation and relevant calibrated values (Cuffey and Paterson, 2010, Eqs. 3.35 and 3.36).

Mass conservation relates the rate of change of the conduit area to the spatial gradient in discharge Q and the production of water due to melt (Fowler, 1999, Eq. 2.2), such that

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial S}{\partial t}}+{\displaystyle \frac{\partial Q}{\partial s}}={\displaystyle \frac{m}{{\mathit{\rho }}_{\mathrm{w}}}}.\end{array}$$

We assume turbulent flow and use Manning's formula to empirically relate the total potential hydraulic gradient (negative basic hydraulic gradient and the effective pressure gradient) to friction along the channel. Note that Manning's formula was derived for freshwater, but, due to a lack of empirical data with saline fluid in ice, we use the Manning's friction factor for freshwater in ice, as used in Fowler (1999). Regardless of salinity, large uncertainties exist in Manning's formula, even in freshwater systems, due to the large variability in the value of Manning's friction factor (Pohle et al., 2022). The equation for the conservation of momentum is then

$$\begin{array}{}\text{(4)}& \mathit{\psi }+{\displaystyle \frac{\partial N}{\partial s}}=f{\mathit{\rho }}_{\mathrm{b}}g{\displaystyle \frac{Q^{\mathrm{2}}}{S^{\mathrm{8}/\mathrm{3}}}},\end{array}$$

where f is the friction factor (Fowler, 1999, Eq. 2.3). For a circular ice-walled channel, $f=(\mathrm{4}\mathit{\pi }{)}^{\mathrm{2}/\mathrm{3}}{\mathcal{R}}^{\mathrm{2}}$ and ℛ=n_i = 0.06 ${\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}\mathrm{s}$, where ℛ is the hydraulic roughness, and n_i is the roughness of ice (Clarke, 2003). For a semi-circular channel at the bed, $f=\left(\mathrm{2}\right(\mathit{\pi }+\mathrm{2}{)}^{\mathrm{2}}\mathit{\pi }{)}^{\mathrm{2}/\mathrm{3}}{\mathcal{R}}^{\mathrm{2}}$, where $\mathcal{R}=\mathit{\pi }/(\mathrm{2}+\mathit{\pi })n_{\mathrm{i}}+(\mathrm{1}-(\mathit{\pi }/(\mathrm{2}+\mathit{\pi }))n_{\mathrm{b}}$, and n_b = 0.16 ${\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}\mathrm{s}$ is the roughness of the bed material (Fowler, 1999, Eq. 2.24).

The equation for the conservation of energy is

$$\begin{array}{}\text{(5)}& \begin{array}{rl}Q\left(\mathit{\psi }+{\displaystyle \frac{\partial N}{\partial s}}\right)=& {\mathit{\rho }}_{\mathrm{b}}{\mathit{\sigma }}_{\mathrm{b}}\left(S{\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{b}}}{\partial t}}+Q{\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{b}}}{\partial s}}\right)+m\mathcal{L}\\ & +m({\mathit{\theta }}_{\mathrm{b}}-\widehat{\mathit{\theta }}),\end{array}\end{array}$$

where θ_b is the temperature of the brine, $\widehat{\mathit{\theta }}$ is the melting point of the ice, σ_b is the specific heat capacity of the brine, and ℒ is the latent heat of the fusion for ice (Fowler, 1999, Eq. 2.4). Following Röthlisberger (1972) and Nye (1976), we neglect the heat transfer equation, which is equivalent to assuming that any heat generated from flow is instantaneously transferred to the channel walls. Consequently, we also need to make an assumption regarding fluid temperature and the melting point, which are related to the salinity of the fluid.

2.2 Consideration of brine

As the brine flows from the subglacial lake, any melt occurring at the channel walls will add freshwater and dilute the brine along the channel. We use a partial differential equation to describe the concentration of salt $\widehat{\mathit{\beta }}$ (kg m⁻³) at position s and time t in response to changes in the channel cross-sectional area and discharge. The fluid moves along the channel at velocity v, which gives the flux of salts per square meter:

$$\mathit{\varphi }=v\widehat{\mathit{\beta }}-D{\displaystyle \frac{\partial \widehat{\mathit{\beta }}}{\partial s}},$$

where D is the diffusion coefficient. The mean velocity of the fluid is given by $v=Q/S$. We calculate a Péclet number of (Pe)>10⁸, which suggests that advection dominates diffusion in fluid flow, and we assume that diffusion is negligible. With this assumption, the flux of salt moving through a channel cross-section with area S is

$$\mathit{\varphi }={\displaystyle \frac{Q\widehat{\mathit{\beta }}}{S}}.$$

Assuming that there is no brine added along the channel and that there is no accretion on the channel walls, the salt concentration equation is

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial }{\partial t}}\left(\widehat{\mathit{\beta }}S\right)={\displaystyle \frac{\partial }{\partial s}}\left(-Q\widehat{\mathit{\beta }}\right).\end{array}$$

In the case of a semi-circular channel, contact with the ground could be a source of salts in the fluid flow. Although we do not include it here, such a mechanism could be accounted for in our model as a source term in Eq. (6).

The salt concentration $\widehat{\mathit{\beta }}$ discussed above is in kg m⁻³ in order to be compatible with the model. These values for salinity are converted to a standard unit for measuring salinity in practical salinity units (psu) before calculating the density and melting point in Eqs. (8) and (9), respectively, using the conversion $\widehat{\mathit{\beta }}$ (kg  m⁻³) = 1000β(ρ_b)⁻¹ (psu). The salinity in the lake is constant in time since no fluid is being added to the lake, which gives the boundary condition $\widehat{\mathit{\beta }}(\mathrm{0},t)=\mathrm{1000}\mathit{\beta }(\mathrm{0},t)({\mathit{\rho }}_{\mathrm{b}}{)}^{-\mathrm{1}}$, where β(0,t) (in psu) is prescribed at the beginning of the simulation. The dilution of the brine along the channel is minimal, and so, at the beginning of the simulation, we assume that the salt concentration in the channel is equal to the concentration in the lake; that is, $\widehat{\mathit{\beta }}(s,\mathrm{0})=\widehat{\mathit{\beta }}(\mathrm{0},t)$.

Fluid properties such as the density, the specific heat capacity, and the melting point of ice are functions of salinity. The specific heat capacity of brine ($\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}$) is calculated with the salinity (psu) and temperature of the brine (°C) in the lake as follows:

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{\mathrm{b}}=\mathrm{4217.4}-\mathrm{3.72}{\mathit{\theta }}_{\mathrm{b}}-\mathrm{7.64}\mathit{\beta },\end{array}$$

where the first-order terms from Eq. (A3.11) in Gill (1982) are used.

As salt concentration changes, the density of brine and the melting point of ice also vary spatially and temporally. The density of the brine (in kg m⁻³) as a function of salt concentration β (psu) under 1 bar using the FREeZing CHEMistry (FREZCHEM) model from Wolfenbarger et al. (2022), disregarding higher-order terms, results in

$$\begin{array}{}\text{(8)}& {\mathit{\rho }}_{\mathrm{b}}=\mathrm{1000}+\mathrm{0.763}\mathit{\beta }.\end{array}$$

Note that we do not account for changes in density due to pressure or temperature. Using the same FREZCHEM model, we calculate the melting point of ice due to salinity and adjust for water pressure. The melting point (in °C) of ice in contact with saline fluid at pressure P_w is

$$\begin{array}{}\text{(9)}& \widehat{\mathit{\theta }}=-c_{\mathrm{b}}\mathit{\beta }-c_{\mathrm{n}}{P}_{\mathrm{w}}=-c_{\mathrm{b}}\mathit{\beta }-c_{\mathrm{n}}(P_{\mathrm{i}}-N),\end{array}$$

where c_b = 6.05 × 10⁻² °C psu⁻¹ and c_n = 7.45 × 10⁻⁸ °C Pa⁻¹. We assume that the lake and surrounding ice system are in thermal equilibrium, which requires it to be that, at the lake, the ice and brine temperatures, θ_i and θ_b, respectively, are equal and at the salinity- and pressure-dependent melting point $\widehat{\mathit{\theta }}$. For a given salt concentration in the lake, we calculate the melting point at the lake and set the ice and brine temperatures to be equal to that temperature at the lake. We assume that the ice temperature remains constant in time and along the channel.

We assume that the brine temperature and the melting point are equal (Röthlisberger, 1972; Werder et al., 2013) and evolve in response to the changes in salinity along the channel and in time. Given that expected changes in the brine temperature are small and that we are modeling slow drainage events, we make the assumption that the brine temperature is in a quasi-steady state ($\frac{\partial {\mathit{\theta }}_{\mathrm{b}}}{\partial t}=\mathrm{0}$). Using the assumption that fluid temperature is equal to the melting point (${\mathit{\theta }}_{\mathrm{b}}=\widehat{\mathit{\theta }}$) along with Eq. (9) gives the simplified equation for the conservation of energy:

$$\begin{array}{}\text{(10)}& Q\left(\mathit{\psi }+{\displaystyle \frac{\partial N}{\partial s}}\right)={\mathit{\rho }}_{\mathrm{b}}{\mathit{\sigma }}_{\mathrm{b}}Q\left(-c_{\mathrm{b}}{\displaystyle \frac{\partial \mathit{\beta }}{\partial s}}+c_{\mathrm{n}}{\displaystyle \frac{\partial N}{\partial s}}\right)+mL.\end{array}$$

The term on the left-hand side of Eq. (5) is the total work done, which must be balanced by the sum of the energy (i) needed to raise the brine temperature to the new salinity- and pressure-dependent melting point and (ii) lost to melting due to latent heat.

We assume a circular channel for most simulations, but we do compare the effect of circular vs. semi-circular channels in Sect. 3. The main differences between these assumptions are that, in the semi-circular case, (1) the fluid is flowing along the bed, and, therefore, the roughness of the bed must be accounted for instead of the roughness of ice, and (2) the substrate may contain some salts. We do not account for (2) in our model. We do account for (1) through the friction factor, which appears in Eq. (4) and changes depending on the channel geometry and roughness of the channel walls or the bed.

2.3 Channel boundary conditions

The only fluid flux from the subglacial lake is the brine flowing out of the channel. Thus, the rate of change of lake volume V is given by

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{d}V}{\mathrm{d}t}}=-Q(\mathrm{0},t).\end{array}$$

We assume a box-shaped lake, which gives the following lake hypsometry:

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{h}{h_{\mathrm{i}}}}={\displaystyle \frac{V}{V_{\mathrm{i}}}},\end{array}$$

where h is the depth of the lake, h_i is the initial lake depth, and V_i is the initial lake volume. For the treatment of more complicated lake geometries, see Kingslake (2013).

Implicitly differentiating gives

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\mathrm{d}V}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}V}{\mathrm{d}h}}{\displaystyle \frac{\mathrm{d}h}{\mathrm{d}t}}={\displaystyle \frac{V_{\mathrm{i}}}{h_{\mathrm{i}}}}{\displaystyle \frac{\mathrm{d}h}{\mathrm{d}t}},\end{array}$$

and, by substitution, the lake depth evolves with time following

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathrm{d}h}{\mathrm{d}t}}={\displaystyle \frac{h_{\mathrm{i}}}{V_{\mathrm{i}}}}(-Q(\mathrm{0},t\left)\right).\end{array}$$

We assume that the lake drains slowly enough that the ice roof drops with the lake depth following Evatt et al. (2006); thus, as the lake drops, the effective pressure at the lake is still the difference between the ice-overburden pressure and the fluid pressure in the lake. The boundary condition where the conduit meets the lake is $N(\mathrm{0},t)=\mathrm{0}$ (Evatt et al., 2006). We impose a Neumann boundary condition at the end of the channel, where

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\partial }{\partial s}}N(s,t){|}_{s=L}=\mathrm{0}.\end{array}$$

We choose this boundary condition (opposed to N=0) in order to solve the system numerically in a more efficient way (see Appendix A for details). Neumann boundary conditions on effective pressure at the end of the channel have been used to solve similar systems of equations without having an influence on the qualitative results (Kingslake, 2015; Evatt et al., 2006).

2.4 Summary of model equations

The full model contains five unknowns ($N,S,m,Q$, and β) and five model equations (Eqs. 2–6) which are solved sequentially. The model equations contain the derived variables $\widehat{\mathit{\theta }}$, ρ_b, and ψ, which depend on salinity. The model equations written in terms of the salinity-dependent derived variables are listed below.

$$\begin{array}{c}{\displaystyle }\text{Channel evolution:}{\displaystyle \frac{\partial S}{\partial t}}={\displaystyle \frac{m}{{\mathit{\rho }}_{\mathrm{i}}}}-KS{N}^{\mathrm{3}}\\ {\displaystyle }\begin{array}{rl}\text{Conservation of energy: }& Q\left(\mathit{\psi }+{\displaystyle \frac{\partial N}{\partial s}}\right)={\mathit{\rho }}_{\mathrm{b}}{\mathit{\sigma }}_{\mathrm{b}}Q\\ & \times \left(-c_{\mathrm{b}}{\displaystyle \frac{\partial \mathit{\beta }}{\partial s}}+c_{\mathrm{n}}{\displaystyle \frac{\partial N}{\partial s}}\right)+mL\end{array}\\ {\displaystyle }\text{Conservation of mass: }{\displaystyle \frac{\partial S}{\partial t}}+{\displaystyle \frac{\partial Q}{\partial s}}={\displaystyle \frac{m}{{\mathit{\rho }}_{\mathrm{b}}}}\\ {\displaystyle }\text{Salt concentration: }{\displaystyle \frac{\partial }{\partial t}}\left(\widehat{\mathit{\beta }}S\right)={\displaystyle \frac{\partial }{\partial s}}\left(-Q\widehat{\mathit{\beta }}\right)\\ {\displaystyle }\text{Conservation of momentum: }\mathit{\psi }+{\displaystyle \frac{\partial N}{\partial s}}=f{\mathit{\rho }}_{\mathrm{b}}g{\displaystyle \frac{Q^{\mathrm{2}}}{S^{\mathrm{8}/\mathrm{3}}}}\end{array}$$

Figure 2 The impact of brine on discharge, velocity, channel radius, and changes in salt concentration and melting points for various initial salt concentrations, shown by the colors in (a). (a) Discharge at the lake outlet as the lake drains. (b) Velocity at the lake outlet over time in days. (c) Channel radius along the length of the channel at the time the lake has emptied. The dashed line is the initial channel radius of 0.25 m. (d) The solid lines (associated with the left axis) indicate the difference between the final salt concentration along the channel after the lake has emptied and the initial salt concentrations shown in the legend of (a). The dashed lines (right axis) refer to the difference along the channel between the melting points at the end of the simulation and initial melting points.

These five equations are non-dimensionalized and are solved numerically, as described in Appendix A. The system of equations is solved using a constant grid spacing of 50 m and a constant time step of less than 1 s for each simulation (smaller time steps are required for higher salinities). After the solution to the salt concentration equation is obtained at each time and space step, the melting point $\widehat{\mathit{\theta }}$ and the density ρ_b are updated along with the basic hydraulic gradient ψ, which is a function of density, using Eqs. (9), (8), and (1), respectively.

3 Results

Idealized model simulations were run to investigate the impact of brine on discharge rates, channel radius, effective pressure, and duration of lake drainage. The model runs until open channel flow occurs, at which point the model run ends. Unless otherwise specified, the parameters used for the baseline simulations are as follows: ice thickness above channel H = 100 m, initial lake volume V_i = 5 × 10⁵ m³, initial lake depth h_i = 10 m, channel length s = 1000 m, initial channel radius r = 0.25 m, bed and conduit slope B = 3°, and circular channel geometry. A range of different values were explored for each parameter listed in Table 2 while holding all other parameters to be equal to the baseline simulation values.

To investigate the effect of saline fluid, we ran six scenarios with $\mathit{\beta }=\mathit{\{}\mathrm{0},\mathrm{5},\mathrm{10},\mathrm{25},\mathrm{50},\mathrm{125}\mathrm{psu}\mathit{\}}$ to explore the range of possible outcomes. Based on the salinity, we set the ice temperature and initial melting point of the ice and initial brine temperature to $\widehat{\mathit{\theta }}=\mathit{\{}-\mathrm{0.06}$, −0.37, −0.67, −1.58, −3.09, −7.63 °C}, respectively. The discharge rates are lower for fluids with higher salt concentrations (Fig. 2a). Higher salt concentrations decrease the peak velocity reached and increase the amount of time required to reach peak velocity (Fig. 2b). The peak velocity and drainage duration change non-linearly with increased salt concentrations. That is, a difference in salinity of 5 psu compared to freshwater has a substantial influence on fluid velocities, while a difference in salinity of the same amount at higher salinities results in much smaller velocity changes (Fig. 2b). After the lake has drained, the channel radius is larger than the initial channel radius for all salinities tested (Fig. 2c). At the end of the simulations, the channel radius is greatest for freshwater, which is more than twice the channel radius of a lake with a salinity of 125 psu.

The salt concentration decreases linearly along the channel in all scenarios, with the most significant changes occurring in cases with higher salinities (Fig. 2d). As the salinity decreases, the melting point increases along the channel (Fig. 2d).

Figure 3 Discharge over time for different channel geometries, lake volumes, initial channel radius, and channel (bed) slope for lakes with freshwater (dashed lines) and 10 psu (solid lines). (a) Channels with semi-circular geometries are denoted in blue, while circular geometries are in yellow. (b) Lake volumes for V_i = 1 × 10⁶ m³ (red line), V_i = 5 × 10⁵ m³ (blue line), and V_i = 1 × 10⁵ m³ (yellow line) are plotted on top of one another. (c) The initial channel radius is varied such that the blue line is r = 0.2 m, the red line is r = 0.3 m, the yellow line is r = 0.4 m, and the purple line is r = 0.5 m. (d) Channel slope is varied by 2° (in yellow), 3° (in blue), and 4° (in red).

We systematically vary parameters to explore the sensitivity of the model and the impact of channel geometry, lake volume, initial channel radius, and bed slope on discharge and the duration of drainage (Fig. 3). The channel geometry (circular vs. semi-circular) changes the time of lake drainage, as well as the peak discharge for freshwater and for a lake with a salinity of 10 psu (Fig. 3a). For a semi-circular channel, the difference in time until lake drainage is more significant than the difference in peak discharge between the freshwater and brine scenarios. For the circular channel, the difference in peak discharge is greater than the difference in the duration of time until lake drainage. In both scenarios, the circular channel drains the lake in less than half the time required by a semi-circular channel, and the peak discharge is over 3 times as high when the channels are circular, which is partially due to the double initial cross-sectional area for those simulations.

The volume of the lake impacts the discharge and the timing of drainage by extending the amount of time the model is run (Fig. 3b). The discharge curves for all lake volumes follow the same curve until the smaller lakes drain. The lake continues to drain for greater lake volumes. For a freshwater lake, the peak discharge roughly triples when comparing V_i = 1 × 10⁵ m³ with V_i = 5 × 10⁵ m³ and V_i = 5 × 10⁵ m³ with V_i = 1 × 10⁶ m³, where, as for a lake with a salinity of 10 psu, the peak discharge approximately doubles. Larger lake volumes extend the time until the lake drains, which results in non-linear increases in the differences between freshwater and saline fluid. The differences in timing and peak discharge for saline and freshwater lakes are largest for greater lake volumes (Fig. 3b).

The initial channel radius impacts the amount of time until the lake drains and the peak discharge (Fig. 3c). With a salinity of 10 psu and an initial channel radius of 0.2 m, the lake drains substantially slower (≈ 10 d) compared to an initial channel radius of 0.3 m. For each 0.1 m increase in the initial channel radius, the peak discharge increases non-linearly for both saline fluid and freshwater, but the rate of increase is less pronounced for saline fluid. Increasing the bed and channel slope increases the peak discharge and decreases the time taken to reach that peak (Fig. 3d). We varied freshwater (β = 0 psu) and brine (β = 10 psu) along with the channel slope and found that, for lower slopes, there is a larger difference in timing between brine and freshwater, and, for greater slopes, there is a larger difference in peak discharge between brine and freshwater. Higher bed slopes lead to higher discharge rates more quickly (Fig. 3d). For all choices of parameters, we find that brine decreases the discharge rates and increases the time until lake drainage compared to freshwater.

Table 2 Prescribed model parameter variables and descriptions with baseline simulation values and ranges explored.

We explored lake depths of 5–15 m, ice thicknesses of 100–1000 m, and channel lengths of 100–5000 m and found that these parameters do not substantially impact the results for the parameter combinations described (data not shown). For all initial conditions on effective pressure, the effective pressure along the channel tends towards zero over time. For a list of parameters and the range of values explored, see Table 2.

As the walls of the channel melt and the brine is diluted, the density of the brine is not constant along the channel. We have accounted for this in our simulations, but neglecting these changes does not have a substantial influence on the results because the changes are very small (data not shown).

4 Discussion

The results of this model suggest that the consideration of brine in relevant glacial systems is important for capturing the dynamics of drainage through ice-walled channels: a failure to consider salt even for low salinities leads to substantially different estimates of channel formation and drainage rates and timescales (Fig. 2). The consideration of salinity is more important when considering systems with high discharge rates (high lake volume and steep bed slopes; see Fig. 3b and d). Large circular channels also tend to lead to more substantial differences in the fluid dynamics between high and low concentrations of salt (Fig. 3a and c).

4.1 Effects of salinity

The presence of salt in fluid decreases the melting point and increases the density of the fluid, both of which have implications for how fluid drains through an R channel from a subglacial lake. An obvious impact is that saline fluid can remain liquid at sub-zero temperatures, and a subglacial lake system can exist below the pressure melting point. The differences related to the salinity-dependent melting point also influence how the channel grows in time. The presence of salt in the system tends to decrease the amount of energy available for melting the channel walls because energy is needed to change the brine temperature to the increasing melting point.

For freshwater systems, a positive feedback allows for dynamic channel growth, where higher discharge rates generate more energy for melting and opening the channels, allowing for even higher discharge rates. This feedback is the mechanism responsible for developing an efficient drainage system with large channels and for glacier lake outburst floods. A system with saline fluid at a sub-zero temperature reduces this positive feedback. The melting of the channel walls results in an increase in the melting point of ice due to the changes in salinity after the addition of fresh meltwater. While melting of the channel wall increases the cross-sectional area and, therefore, the discharge, more energy is now required to increase the brine temperature to the new melting point, and less energy is available for melting.

The ratio of fresh meltwater from the channel walls to saline discharge from the lake determines brine concentration. A channel with a smaller cross-sectional area has a higher surface-area-to-volume ratio compared to a larger channel, leading to a higher ratio of meltwater to saline discharge and larger changes in brine concentration along the channel. Channels are smaller for higher salinities (Fig. 2c). Brine concentration and the resulting changes in the melting point vary linearly along the channel, with larger gradients observed at higher salinities (Fig. 2d). Larger spatial gradients in brine concentration require more energy to raise the brine temperature to the new salinity- and pressure-dependent melting point along the channel. This reduces the energy available for melting, which strongly inhibits rapid channel growth (Eq. 10). The effect of inhibited channel growth is largest at higher salinities, where the initial brine temperature is lowest (−7.63 °C). However, the system is more sensitive at lower salinities, where even small increases in salinity lead to significant changes, as the initial brine temperature is close to 0 °C (−0.37 to −0.67 °C) (Fig. 2c). Inhibited channel growth due to salinity results in more gradual increases in velocity over time (Fig. 2b).

Figure 4 (a) Discharge curves for drainage from a temperate freshwater system with densities equivalent to those of freshwater and brine with a salinity of 125 psu. (b) Discharge curves for brine with β = 5, 10, 25, 50, 125 psu but all with a constant density equal to that of freshwater (1000 kg m⁻³). Simulations from the full model where both density and the effects of saline fluid at sub-zero temperatures are considered are plotted with the dashed lines, which are hardly visible.

Density impacts channel growth in the opposite way. A fluid with a higher salt concentration has a higher density. As a denser material moves through a gravitational potential, more energy is generated and available for melting. In Fig. 4a, a temperate freshwater system ($\widehat{\mathit{\theta }}$ = −0.06 °C) is modeled with the density of freshwater and a fluid density of ρ = 1098 kg m⁻³, equivalent to the density of brine with a salinity of 125 psu. The fluid with the higher density results in a higher peak discharge (Fig. 4a).

Modeling all other changes related to salinity (including the treatment of the melting point and fluid temperature) while holding the density of the brine constant and equal to that of freshwater (ρ = 1000 kg m⁻³) results in the discharge rates shown in Fig. 4b. Without the consideration of accurate brine densities, there is almost no difference in the peak discharge, even for the highest salinities with the highest densities (Fig. 4b). Higher fluid densities amplify the dynamic feedback, leading to channel growth and higher discharge rates (Fig. 4a). However, when the channel growth is reduced by the effects of salinity; sub-zero initial brine temperatures; and, thus, changes in the melting point, this dynamic feedback does not occur, and the influence of density is limited. In part, this is due to the fact that density changes almost linearly with salinity (Eq. 8), while the impact of salinity on discharge is non-linear.

It is important to note that the increased density generates more energy for melting only when the flow is gravity driven (i.e., a downward-sloped inclined channel). For example, if the flow of a saline fluid is driven uphill by ice-overburden pressure, the fluid velocity will be slower compared to the fluid velocity of freshwater.

Due to the substantial influence of changes in the melting point on energy available for melting, the considerations of thermodynamics and assumptions around the brine temperature have large impacts on the results when modeling saline fluid in channelized systems. Assuming constant brine temperature (which is opposed to assuming that the brine temperature remains equal to the melting point, as we do here) results in substantially different results, where the influences of density become the first-order effect. However, this assumption is likely to be unrealistic and differs from the convention in previous studies. In other models of subglacial channel flow that do not explicitly include temperature, the temperature is set to follow the pressure melting point of the ice walls within the channel (e.g., Röthlisberger, 1972; Werder et al., 2013).

4.2 Drainage types of saline systems

There are limited observations of saline outflows and subglacial lake drainage events, which makes it difficult to understand the hydrological state of such systems. Each drainage system has a unique combination of (i) sources and concentrations of salt, (ii) brine and ice temperatures, (iii) bed topography and glacier geometry determining the hydraulic potential, and (iv) the mechanism initiating drainage, all of which play a role in the state of the subglacial system.

Our model results show that, when the brine temperature remains equal to the melting point, the sub-zero saline fluid reduces the positive feedback that leads to the rapid channel growth and increases in fluid velocity. Additionally, our results show that effective pressures along the entire channel tend towards zero over time; that is, the modeled channel is a high-water-pressure system. Consistent high water pressure and reduced channel growth may suggest that saline systems do not tend toward channelization and may exist as distributed systems. However, further modeling of a distributed system with sub-zero saline fluids would be required to characterize such behavior.

Alternatively, Badgeley et al. (2017) hypothesize that the mechanism of drainage at Blood Falls is not by means of the opening and closing of R channels but rather by means of brine injection into basal crevasses. More observational and modeling work is needed to understand more about how saline fluid behaves in a spatially connected subglacial hydrologic system.

4.3 Implications for outburst floods

As shown in Sect. 4.1, density can have a significant effect on fluid flow regardless of salinity. Outburst floods often result in disproportionate amounts of suspended sediment, which increases the density of the water (Snorrason et al., 2002; Church, 1972). Discharge from outburst floods is typically on the order of 100–1000 m³ s⁻¹ (Walder and Costa, 1996, Table 1) and can contain suspended sediment concentrations (SSCs) of 70.7 g L⁻¹ (Beecroft, 1983; Old et al., 2005) and, in some extreme cases, over 400 g L⁻¹ (Maizels, 1997). The density of sediments ρ_s depends on the rock type and clast size but typically ranges from 2350 to 2760 kg m⁻³ (Frederick et al., 2016; Guan et al., 2015; Chikita, 2004). The combined fluid density ρ_c of water with suspended sediments is related to suspended sediment concentrations as follows:

$$\begin{array}{}\text{(16)}& {\mathit{\rho }}_{\mathrm{c}}={\mathit{\rho }}_{\mathrm{w}}+({\mathit{\rho }}_{\mathrm{s}}-{\mathit{\rho }}_{\mathrm{w}}){\displaystyle \frac{\mathrm{SSC}}{{\mathit{\rho }}_{\mathrm{s}}}}.\end{array}$$

We model an outburst flood from a subglacial lake with freshwater at 0 °C and at a volume of V_i=1x10⁷ m³, with all other parameters being equal to those in the baseline simulation (Table 2). We vary the suspended sediment concentrations to be $\mathrm{SSC}=\mathit{\{}\mathrm{0},\mathrm{63},\mathrm{125},\mathrm{188},\mathrm{251}\mathit{\}}$ g L⁻¹ to arrive at the combined fluid densities of ${\mathit{\rho }}_{\mathrm{c}}=\mathit{\{}\mathrm{1000},\mathrm{1040},\mathrm{1080},\mathrm{1120},\mathrm{1160}\mathit{\}}$ kg m⁻³ to simulate different suspended sediment loading (for ρ_s = 2760 kg m⁻³). There is a significant difference between the peak discharge of floods with a lower-density fluid (Q≈ 115 m³ s⁻¹) and that of floods with a higher-density fluid (Q≈ 140 m³ s⁻¹), as well as in the timing of the peak discharge (Fig. 5). Neglecting to account for sediment loading and accurate water densities could lead to inaccurate results when modeling outburst floods where gravity assists in driving fluid flow. The influence of density is particularly important for systems with high gravitational potential energy and high fluid density, which, together, can significantly increase fluid flow rate.

Figure 5 Outburst flood hydrographs for water with different suspended sediment concentrations and, therefore, varying fluid densities.

5 Conclusions

We have presented a subglacial hydrology model which includes the consideration of saline fluid. Salt allows fluid to exist below the freezing point of freshwater and increases the density of the fluid. We show that, if a channel exists, hypersaline fluid can flow through an ice-walled channel when the brine and ice are at the salinity- and pressure-dependent melting point. Our results suggest that a higher salt concentration and commensurate lower brine temperature decrease the peak discharge and increase the time for a fixed volume to drain from a subglacial lake.

The main driver of the decreased discharge rates for higher salinities is the influence of changes in the salinity-dependent melting point. As the salinity along the channel changes due to melting of the channel walls, the melting point changes and, consequently, so does the energy needed to melt the channel walls. While more energy is generated when a higher-density fluid moves through a gravitational potential, this increase in energy for higher salinities is minimal compared to the energy needed to raise the brine temperature to the melting point.

This study shows that accounting for fluid properties is crucial for accurately modeling subglacial hydrology in relevant systems (i.e., saline fluid, more or less dense fluid). For a lake with a salinity of 125 psu, which is approximately the measured value at Blood Falls, and an initial brine temperature of −7.63 °C, the peak velocity reached is 40 % lower and the lake drains 9 % slower than for a freshwater lake. The duration of drainage is most sensitive to the initial channel radius and channel geometry, while peak discharge is most sensitive to lake volume, channel slope, and channel geometry (circular vs. semi-circular). We explored the influence of varying the fluid density in relation to suspended sediment loads on outburst floods and found that peak discharge is significantly higher for a high-density fluid (22 % higher than pure water when the fluid density is 1160 m³ s⁻¹). In our model, sub-zero saline fluid reduced channel growth, which may have implications for the type of subglacial drainage system (channelized vs. distributed) that forms with such fluids.

We make a number of simplifying assumptions in the model and use arbitrary parameters for ice thickness, lake volume, channel length, bed slope, and initial channel radius due to a lack of available data on subglacial hypersaline systems. We also impose brine temperatures that remain equal to the ice melting point. In this light, the model presented here should be viewed as an initial exploration of the impact of brine on the dynamics of a subglacial hydrological system. Additional modeling efforts are needed to provide a thorough sensitivity and stability analysis. Further research is required to understand the initiation of drainage in cold saline environments and the influence of fluid density on drainage networks and outburst floods.

Appendix A: Numerics

After non-dimensionalization, the model equations to be solved numerically are the following.

Dimensionless model equations

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}\text{Channel evolution: }{\displaystyle \frac{\partial S}{\partial t}}=& Q\left(\mathit{\psi }+{\mathit{\gamma }}_{\mathrm{1}}{\displaystyle \frac{\partial N}{\partial s}}+{\mathit{\gamma }}_{\mathrm{2}}{\displaystyle \frac{\partial \mathit{\beta }}{\partial s}}\right)\\ & -S{N}^{\mathrm{3}}.\end{array}\text{(A2)}& {\displaystyle }\text{Conservation of mass: }{\displaystyle \frac{\partial Q}{\partial s}}=\mathit{ϵ}\left(rm-{\displaystyle \frac{\partial S}{\partial t}}\right).\text{(A3)}& {\displaystyle }\text{Conservation of momentum: }{\displaystyle \frac{\partial N}{\partial s}}={\displaystyle \frac{\mathrm{1}}{\mathit{\delta }}}\left({\displaystyle \frac{Q^{\mathrm{2}}}{S^{\mathrm{8}/\mathrm{3}}}}-\mathit{\psi }\right).\text{(A4)}& {\displaystyle }\text{Salt concentration: }\mathit{\lambda }{\displaystyle \frac{\partial }{\partial t}}\left(\widehat{\mathit{\beta }}S\right)={\displaystyle \frac{\partial }{\partial s}}\left(-\widehat{\mathit{\beta }}Q\right).\text{(A5)}& {\displaystyle }\text{Boundary conditions: }N(\mathrm{0},t)=\mathrm{0},\text{(A6)}& {\displaystyle }{\displaystyle \frac{\partial }{\partial s}}N(s,t){|}_{s=\mathrm{1}}=\mathrm{0},\text{(A7)}& {\displaystyle }\widehat{\mathit{\beta }}(\mathrm{0},t)=\widehat{\mathit{\beta }}(\mathrm{0},\mathrm{0}),\text{(A8)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}h}{\mathrm{d}t}}=\mathit{\zeta }Q(\mathrm{0},t).\end{array}$$

Model and scaling parameters:

$$\begin{array}{c}{\displaystyle }{N}_{\mathrm{0}}=(K{t}_{\mathrm{0}}{)}^{-\mathrm{1}/\mathrm{3}},m_{\mathrm{0}}={\displaystyle \frac{Q_{\mathrm{0}}{\mathit{\psi }}_{\mathrm{0}}}{\mathcal{L}}},S_{\mathrm{0}}={\left({\displaystyle \frac{f{\mathit{\rho }}_{\mathrm{b}}g{Q}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{\psi }}_{\mathrm{0}}}}\right)}^{\mathrm{3}/\mathrm{8}},\\ {\displaystyle }{t}_{\mathrm{0}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}{S}_{\mathrm{0}}}{m_{\mathrm{0}}}},r={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{b}}}},\mathit{ϵ}={\displaystyle \frac{s_{\mathrm{0}}{m}_{\mathrm{0}}}{Q_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{i}}}},\\ {\displaystyle }\mathit{\delta }={\displaystyle \frac{N_{\mathrm{0}}}{s_{\mathrm{0}}{\mathit{\psi }}_{\mathrm{0}}}},{\mathit{\gamma }}_{\mathrm{1}}={\displaystyle \frac{N_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{b}}{\mathit{\sigma }}_{\mathrm{b}}{c}_{\mathrm{n}}}{s_{\mathrm{0}}{\mathit{\psi }}_{\mathrm{0}}}},{\mathit{\gamma }}_{\mathrm{2}}={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{b}}{\mathit{\sigma }}_{\mathrm{b}}{c}_{\mathrm{b}}}{s_{\mathrm{0}}}},\mathit{\lambda }={\displaystyle \frac{S_{\mathrm{0}}{s}_{\mathrm{0}}}{t_{\mathrm{0}}{Q}_{\mathrm{0}}}},\\ {\displaystyle }\mathit{\zeta }={\displaystyle \frac{t_{\mathrm{0}}{h}_{\mathrm{i}}{Q}_{\mathrm{0}}}{p{V}_{\mathrm{i}}{h}_{\mathrm{0}}}}.\end{array}$$

We use the subscripts $j=\mathrm{1},\mathrm{2},\mathrm{\dots }n$ to denote the grid points along the channel which are separated by Δs and the superscripts $i=\mathrm{1},\mathrm{2},\mathrm{\dots }m$ to denote time steps separated by Δt. Note that the equation for the conservation of energy has been substituted into the channel evolution equation to give Eq. (A1).

To solve Eq. (A8), we follow Kingslake (2013) in using the forward Euler method to evolve the lake depth forward in time.

$$h^{i+\mathrm{1}}=h^i+\mathrm{\Delta }t\mathit{\zeta }{Q}_{\mathrm{1}}^i$$

Similarly, we solve Eq. (A1) using the same method. The channel cross-sectional area S is moved forward in time at all grid points using

$$\begin{array}{rl}{S}_j^{i+\mathrm{1}}=& S_j^i+\mathrm{\Delta }t(Q_j^i\left(\mathit{\psi }+{\mathit{\gamma }}_{\mathrm{1}}{\displaystyle \frac{N_j^i-N_{j-\mathrm{1}}^i}{\mathrm{\Delta }s}}+{\mathit{\gamma }}_{\mathrm{2}}{\displaystyle \frac{{\mathit{\beta }}_j^i-{\mathit{\beta }}_{j-\mathrm{1}}^i}{\mathrm{\Delta }s}}\right)\\ & -S_j^i\left(N_j^i{)}^{\mathrm{3}}\right)\end{array}$$

for $j=\mathrm{1},\mathrm{2},\mathrm{\dots }n$.

To evolve these two equations forward in time, the discharge and effective pressure at time step i are needed. These variables can be found simultaneously by solving the mass and momentum equations.

We follow Fowler (1999) and Kingslake (2013) in assuming that ϵ is small enough to neglect the terms containing ϵ in Eq. (A2). With parameter values m₀, s₀ = 1000 m, and ρ_i = 917 kg m⁻³, ϵ is on the order of 10⁻³, and, thus, we neglect these terms, which simplifies Eq. (A2) to

$$\begin{array}{}\text{(A9)}& {\displaystyle \frac{\partial Q}{\partial s}}=\mathrm{0}.\end{array}$$

This is equivalent to assuming that any melt generated along the channel is small in comparison to the volume of fluid moving through the channel from the lake. Solving Eq. (A3) for the discharge Q and evaluating at the end of the channel gives

$$\begin{array}{}\text{(A10)}& Q^{\mathrm{2}}(\mathrm{1},t)=\left(\mathit{\delta }{\displaystyle \frac{\partial N}{\partial s}}(\mathrm{1},t)+\mathit{\psi }(\mathrm{1},t)\right)S(\mathrm{1},t{)}^{\mathrm{8}/\mathrm{3}}.\end{array}$$

From Eq. (A6) and considering the assumption that the discharge is always positive (flowing out of the lake),

$$\begin{array}{}\text{(A11)}& Q(\mathrm{1},t)=\sqrt{S(\mathrm{1},t{)}^{\mathrm{8}/\mathrm{3}}\mathit{\psi }(\mathrm{1},t)}.\end{array}$$

From Eqs. (A11) and (A9), $Q(s,t)=Q(\mathrm{1},t)$, and we arrive at the following equation for discharge as a function of time:

$$Q\left(t\right)=\sqrt{S(\mathrm{1},t{)}^{\mathrm{8}/\mathrm{3}}\mathit{\psi }(\mathrm{1},t)}).$$

We solve this equation by calculating the discharge profile at each grid point using

$$\begin{array}{}\text{(A12)}& Q_j^i=\sqrt{(S_n^i{)}^{\mathrm{8}/\mathrm{3}}{\mathit{\psi }}_n^i}.\end{array}$$

To calculate the effective pressure along the channel, we start with the boundary condition at the lake given in Eq. (A5) and use Eq. (A3) to iterate

$$\begin{array}{}\text{(A13)}& N_{j+\mathrm{1}}^i=N_j^i+{\displaystyle \frac{\mathrm{\Delta }s}{\mathit{\delta }}}\left({\displaystyle \frac{(Q_j^i{)}^{\mathrm{2}}}{(S_j^i{)}^{\mathrm{8}/\mathrm{3}}}}-{\mathit{\psi }}_j^i\right)\end{array}$$

from $j=\mathrm{1},\mathrm{2},\mathrm{\dots }n$.

To solve the dimensionless brine equation of Eq. (A4), we use an upwind difference scheme such that

$${\widehat{\mathit{\beta }}}_j^{i+\mathrm{1}}={\widehat{\mathit{\beta }}}_j^i-\mathrm{\Delta }t\left({\displaystyle \frac{{\mathit{\beta }}_j^i}{S_j^i}}\left({\displaystyle \frac{S_j^i-S_j^{i-\mathrm{1}}}{\mathrm{\Delta }t}}\right)+{\displaystyle \frac{Q_j^i}{\mathit{\lambda }{S}_j^i}}\left({\displaystyle \frac{{\widehat{\mathit{\beta }}}_j^i-{\widehat{\mathit{\beta }}}_{j-\mathrm{1}}^i}{\mathrm{\Delta }s}}\right)\right).$$

After calculating the non-dimensional salt concentration, we re-dimensionalize the salt concentration to be in units of kg m⁻³ and convert this to psu using

$$\begin{array}{}\text{(A14)}& {\mathit{\beta }}_j^{i+\mathrm{1}}={\displaystyle \frac{{\widehat{\mathit{\beta }}}_j^{i+\mathrm{1}}{\widehat{\mathit{\beta }}}_{\mathrm{0}}\mathrm{1000}}{{{\mathit{\rho }}_{\mathrm{b}}}_j^i}}.\end{array}$$

The density of brine along the channel can be updated with the new salt concentration using Eq. (8). Similarly, the basic hydraulic gradient is updated using the new brine density. In order to solve the system, initial conditions are needed for cross-sectional area, effective pressure, discharge, and salinity along the channel. At the beginning of the simulation, we assume a constant cross-sectional area and salt concentration along the channel. We use the initial cross-sectional area at the end of the channel to calculate the initial discharge curve along the channel using Eq. (A12). We assume that the initial effective pressure profile linearly increases from N=0 at the lake to N=P_i at the end of the channel.