Groundwater interaction module (Link storage routing) - SRG

Background

Interactions between the underlying groundwater and surface water systems affect flow and constituents in these systems. The Groundwater Interaction module enables the analysis of these effects. It determines the exchange flux of water between a river and the underlying aquifer for each link of Source at any time-step. The estimated flux accounts for all the interactions between groundwater (GW) and surface water (SW) along the entire length of the link. The direction of the flux can either be from the river to aquifer or vice versa, that is, the river either loses water to the groundwater or it gains water from the groundwater system. The Groundwater module also calculates salinity (as Total Dissolved Solids, TDS), which is generally a time-invariant constant for the link.

Figure 1 shows the exchange of data between the Groundwater module and the Source Link module.

The model is applied at the link scale. The model assumes a daily time-step.

Version

Source v2.10

Assumptions and constraints

The following assumptions are made:

• The river always interacts with the underlying aquifer. This interaction may be different for every link and is time variable.
• A linear system is assumed whereby the overall GW-SW exchange along each link at any time-step can be estimated by summing the individual interactions resulting from any of the relevant GW-SW processes.
• The type of connection between a river and the underlying aquifer is assessed for each link based on the GW-SW connectivity mapping for the region; it remains constant during the simulation period.
• The aquifer underlying a river is single-layered, unconfined, and has an infinite boundary located at the opposite side of the river.
• Aquifer transmissivity is assumed to be constant along a link (homogeneous aquifer).
• The hydraulic conductance of the river-aquifer interconnection is assumed to be constant along a link but can be different for every link.
• The GW-SW exchange flux is estimated based on best knowledge of initial parameter values; those parameters should be refined during the routing calibration process.

Theory

The hydraulic connection between a river and the underlying aquifer (as informed by connectivity mapping) can either be saturated or unsaturated depending on the location of the water table relative to the river stage. It is useful to classify the connection types using this criterion (rather than gaining/losing) because the flux estimation method and the relevant processes are different for each of those two connection types. Depending on the type of connection, the exchange flux can be attributed to one or more of the following relevant processes:

Natural exchanges between the river and the aquifer as driven by the difference in river and water table levels during and outside storm events; note that this process continually takes place and this is why some level of interaction with the groundwater system always exists:

• Exchange to groundwater recharge (eg. land use change, irrigation)
• Exchange due to groundwater pumping; and
• Exchange due to evapotranspiration.

A linear system is assumed where the cumulative exchange flux is the summation of the individual fluxes resulting from each of the three processes outlined above. The flux can either be directly estimated based on one of the following methods, depending on the data available:

• Head method: Uses knowledge of river conductance and the difference between river stage and the water table within its close vicinity;
• Flux method: Inferred from knowledge of active stresses in the aquifer; or
• A hybrid approach combining both methods.

The level of spatial and numerical complexity of groundwater process represented within the Groundwater - Surface Water Link will be dependent on:

• The type of connection between the river and the underlying aquifer, and
• Data availability, which includes river stage heights, water table levels, and details of the active stresses in the aquifer (eg, pumping).

The appropriate method for flux estimation, whether it be head- or flux-based, is dictated by data availability and connection type; in fact the three (methodology, connection type, and data availability) are inter-dependent.

For example, in situations where the water table level is well below the base of the river such that an unsaturated zone exists between it and the groundwater system, the flux will be predominantly a loss (river recharging the aquifer) and is estimated using a head-based approach. The flux-groundwater level relationship can either be linear or non-linear (up to a maximum loss rate) depending on the length of the unsaturated zone below the river. In this case, only the natural interaction becomes relevant, that is, the type of connection has deemed all the data related to all the other processes (which are now irrelevant) as unnecessary.

At the other end of the spectrum are floodplain situations where the shallow groundwater beneath the floodplain is fully connected to the river. In those cases, predictions of groundwater loss by evapotranspiration, stream depletion due groundwater pumping, bank storage exchange, and the infiltration of overbank flows through wetlands and recharge beneath irrigation areas down to the groundwater and ultimately back to the river all become very relevant. This implicitly indicates higher data requirements including information on surface soil type, vegetation type, flood inundation mapping, bore locations and pumping volumes, in addition to data relating to surface water flow and water table levels. However, if time-series of river stage heights and water table levels for the entire river model simulation exists (thus allowing the head-based method to be used), that precludes the need for detailed data to infer the flux from individual stresses.

Hence, one needs to specify the type of connection for each link and thus identify relevant processes, assess data availability to decide which method is appropriate, and finally collate the required data and proceed with flux estimation.

Types of connections between rivers and acquifers

The types of river-aquifer connections are shown in Figure 2. Rivers are commonly grouped as gaining or losing, whereby the former gains water from the aquifer (Figure 2A) and the latter loses water to the aquifer (Figure 2B-D). However, for the purpose of this model, the river-aquifer connections are grouped into saturated and unsaturated. A gaining river must have a saturated connection with the underlying aquifer (Figure 2A) whereas a losing river might have either a saturated or an unsaturated connection. Figure 2B shows a losing river with saturated connection and Figure 2D shows a losing river with an unsaturated connection. Saturated conditions can still be maintained even after the water table drops below the riverbed level provided that water table level remains higher than some depth D_un (Figure 2C) below the riverbed.

D_un is a function of K₁/K₂ and B/D, where K₁ and K₂ are the hydraulic conductivities of the aquifer and riverbed sediments, respectively, and B/D is the ratio of river width to depth to water table (at some large distance from the river). The exchange flux continues to vary linearly with head gradient as long as the water table is less than D_un below the riverbed (Figure 3, depth < D_un).

When the depth to the water table is larger than D_un below the riverbed, an unsaturated zone starts to develop whereby the hydraulic conductivity starts to decline rapidly thus becoming a function of the soil’s suction. During that stage the flux varies in a non-linear manner with depth to the water table (Figure 3, D_un < depth < D_ml). When the depth to the water table exceeds D_ml, the flux attains a maximum value and becomes independent of depth to water table (Figure 3, depth > D_ml below the riverbed); unit-gradient, gravity flow prevails thereafter.

Classifying the types of river-aquifer connections into saturated and unsaturated is useful because different methods for estimating the exchange flux are needed for each connection type. In the next section, methodologies for evaluating the exchange flux between groundwater and surface water for various types of connections will be presented.

Processes contributing to the GW-SW exchange flux

Natural River-Aquifer Interaction: A river is in continuous interaction with the underlying aquifer. This interaction can be discretised into three components depending on the status of flow in the river:

Interaction during base-flow (non-event, low-flow) conditions;

• Interaction during within-bank flood events; and
• Interaction during overbank flood events.
• The overall natural exchange between the river and the underlying aquifer can be estimated by superimposing the fluxes resulting from the three individual components as follows:

Equation 1

Interaction during low-flow (non-event) conditions: a river may continuously lose water to, or gain water from the underlying aquifer (as shown in Figure 2); neutral cases are also a possibility with no head gradient and zero exchange. If we consider a regional river-aquifer system that is at a steady-state (ie. recharge into the entire aquifer is equal to discharge from the aquifer to the river system), this exchange flux (gain or loss) would remain constant with time. It is given by:

Equation 2

where:

Q_l is the GW-SW exchange flow rate at low flow conditions [L³/T]

h_wt is the water table level [L],

h_r is the river stage level [L],

Δh is the head difference [L],

K is the hydraulic conductivity of the riverbed sediments [L/T],

L is the length of the river link [L],

W is the width of the river link [L],

M is the thickness of the riverbed sediments [L], and

C is the hydraulic conductance of the river-aquifer interconnection [L²/T].

Interaction during within-bank flow events (bank storage): river stage rises as the river flow increases following rainfall events. This triggers a change or reversal in head gradient between the river and the aquifer, which results in water infiltrating into the aquifer. When the river stage returns back to pre-event conditions, this water would subsequently return to the river. The significance of bank storage varies with the size and hydraulic properties of the floodplain that a river interacts with. A highly transmissive floodplain allows a large amount of water to infiltrate into the floodplain, which quickly returns back to the river after the flood-wave subsides; the opposite occurs in the case of low transmissivity floodplains. Depending on the purpose of the modelling and the response of the river-floodplain system in question, one can either account for this phenomenon or not. The bank storage phenomenon is only relevant to river-aquifer systems with a saturated connection. Figure 4 shows a typical flux response for an idealized flood wave where C and Φ represent conductance and storage coefficient of the aquifer, respectively.

Here, we use the formulation derived by Moench et al (1974) and Hall and Moench (1972) to estimate the bank storage flux at any time-step (as implemented by Birkhead and James, 2002):

Equation 3

where:

Q_s is the bank storage flow rate [L/T; discharge from the bank to the river is positive, consistent with sign convention for a gaining river],

T is the aquifer transmissivity [L²/T],

S is the aquifer specific yield,

b is the lateral extent of the floodplain [L]; note that bank storage is effective on both sides of the river bank and as the lateral extent of the floodplain may vary on either side, Equation 3 needs to be implemented independently for each side of the river.

L is the length of the link [L],

h is the river stage height at any step i [L],

dt is the time-step [T],

t is the time level, and

n is the number of times the exponential term needs to be summed; pilot runs have indicated diminishing benefits in accuracy for n>25.

Interaction during over-bank flow events: overbank flow during major flood events results in large quantities of water that travels as sheet flow across the floodplain. After the flood-wave subsides, the overbank water returns back to the river minus the portion that had evaporated and infiltrated through the floodplain surface. Some of the overbank water might be trapped in depressions whereby some of this water will evaporate and the remainder will continue to infiltrate long after the flood-wave subsides. All the infiltrated water travels through the unsaturated zone and eventually recharges the aquifer; this recharge is time variant. This recharge increases river discharge after a time lag, which depends on the orthogonal distance to the river (a) and the aquifer diffusivity (D). The over-bank flow events are mostly relevant to river-aquifer systems with a saturated connection because the lateral flow will not eventuate as river discharge in systems with an unsaturated connection. It is worthwhile noting that very large flood events may significantly raise the water table with the potential to change the connection type between the river and the aquifer. This condition should be assessed by the modeler, and in the event of it happening, then a different scenario should be run with a changed connection type.

The proposed algorithm calculates the evaporative flux from an inundated floodplain area in addition to the infiltrative flux passing through the floodplain bed while maintaining the mass balance of the inundated floodplain area. The evaporative flow rate is given by:

Equation 4

where:

f_ev is the evaporative flow rate [L³/T],

W_a is the inundated floodplain area [L²],

PE is the open water surface evaporation rate [L/T], and

P is the rainfall [L/T].

The infiltrative flow rate is given by:

Equation 5

where:

f_inf(t) is the infiltrative flow rate [L³/T],

W_a is the inundated floodplain area [L²],

K_s is the hydraulic conductivity of the wetland bed [L/T], and

Δh_w is the head gradient in the wetland

The volume of the water stored in the inundated floodplain area (W_v) and water level is calculated at the end of every time-step (t) as follows:

Equation 6

Equation 7

where:

W_v(t) is the volume of the water stored in the inundated floodplain area at the current time-step [L³],

W_v(t-1) is the volume of the water stored in the inundated floodplain area at the previous time-step [L³],

h_w(t) is the water head in the inundated floodplain area at the current time-step [L], and

dt is the time-step [T].

The time series for f_inf(t) becomes a recharge time series (after a time lag in the unsaturated zone, which is neglected in this formulation). The length of the recharge time series depends on the initial volume of water stored in the inundated floodplain area and the rates of evaporation and infiltration, that is, the time series ends when W_v becomes zero (the inundated floodplain area empties). The discharge resulting from an instantaneous recharge source f_inf(t) is given by (Knight et al, 2005):

Equation 8

where:

D_ob(t) is the instantaneous river discharge time series resulting from aquifer recharge sourced from overbank flow [L³/T],

f_inf(t) is the instantaneous recharge sourced from overbank flow during one time-step [L³/T],

D is aquifer diffusivity [L²/T],

a is the orthogonal distance from the centre of the inundated floodplain area to the river [L], and

t is time [T].

The overall discharge is estimated by summing up the individual discharge response for every recharge value at any time, which is given by:

Equation 9

where:

CuQ_ob(t) is the cumulative discharge due to overbank return flow [L³/T],

D_ob is the discharge due to overbank return flow given by Equation 8,

n is a counter for the number of elements in f_inf time series,

m is the total number of elements in the  f_inf(t) time series, and

t is the time [T].

For large flood events that have the potential to change the type of connection between a river and the underlying aquifer, the user should identify the time at which such conditions occur then set a new run with the updated connection type of the relevant river links.

Groundwater pumping

Groundwater pumping is one of the most important processes that impact the exchange flux between groundwater and surface water. Pumping-induced river depletion is defined as the reduction of river flow due to induced infiltration of stream water into the aquifer or the capture of aquifer discharge to the river. This process is only relevant to rivers connected to the aquifer via fully saturated material.

During pumping, while the cone of depression progresses towards a nearby river, groundwater depletion occurs. When this cone reaches the river, the rate of groundwater discharge to the stream reduces, and surface water may even start to infiltrate the aquifer thus marking the start of river depletion. After a long period of pumping, the cone of depression takes its final shape (ie. steady-state is reached), and a portion, or in some cases all, of the pumping will be balanced by a reduction in or reversal of flow from the aquifer to the river. The proportion of pumping met by river depletion in the steady-state case will depend on various factors, including the proximity of the bore to the river compared to the distance between the bore and other stresses to the groundwater system (eg. recharge, ET).

The time required to reach steady-state state varies with the square of the distance between the river edge and the pumping well and varies linearly with aquifer diffusivity. Other important factors that may significantly affect river depletion include riverbed clogging as quantified by the riverbed-aquifer hydraulic conductivity contrast, degree of stream partial penetration, and aquifer heterogeneity. Identifying the time to reach steady-state state is important. There are numerous analytical solutions for river depletion derived for a variety of conceptual systems for pumping; a comprehensive study of those solutions is given by Rassam and Werner (2008). Here, we implement the classical Glover and Balmer (1954) solution for cumulative discharge response in unconfined aquifers, which is given by:

Equation 10

where:

D_p is the depletion discharge rate resulting from groundwater pumping (the fraction of pumped water being sourced from the river) [L3/T],

P is the pumping rate [L³/T],

D is aquifer diffusivity [L²/T],

a is the orthogonal distance between the river and the pump [L],

t is time [T], and

erfc is the complementary error function where erfc(x) = 1-erf(x) and the error function is given by

Equation 11

where m=50 provides adequate accuracy.

Assuming linearity and superposition (Rassam et al, 2004), a number of pumping sources can be modeled with the overall depletion flow rate estimated using:

Equation 12

where:

CuQ_p(t) is the cumulative depletion flux resulting from a number of pumps [L³/T],

D_p(t) is the depletion flux for an individual pump [L³/T] given by Equation 10, and

np is the number of pumping sources active along a link.

The classical Glover and Balmer (1954) solution is valid for a semi-infinite domain where the aquifer is assumed to be infinite in the direction opposite to the river. In areas where a no-flow boundary exists at less than four times the distance to the river, the Glover and Balmer solution would underestimate response. The solution proposed by Knight et al (2005) accounts for the effect of a no-flow boundary, which is given by:

Equation 13

where c represents the distance from the no-flow boundary to the river; sufficient accuracy is obtained when n is equal to 20.

A pump (or a cluster of pumps) may impact one or more links in Source. The flux can be apportioned according to its relative location between two river links. The location of each pump is identified by two coordinates, the orthogonal distance to the river (a) and its position relative to two consecutive nodes defined by a distance (y). The flux can be apportioned as follows (Rassam et al, 2004):

Equation 14

where:

f is the instantaneous flux at any lateral distance from a node per unit length of the river [L²/T]; to obtain the total flux across a link, Equation 14 needs to be integrated along its length,

y is the lateral distance from the upstream node of the link [L],

D is aquifer diffusivity [L²/T],

a is the orthogonal distance between the river and the pump [L], and

t is time [T].

For the case of a semi-confined aquifer, stream depletion may be evaluated as follows (Hunt, 2003):

Equation 15

where:

α is an integration variable,

λ is the streambed resistance coefficient, and

F(α,t) and G(α,t) are functions defined by Equations 47-52 in Hunt (2003).

Stream depletion according to Hunt (2003) is a function of four non-dimensional variables, namely: t*, K*, λ*, and ε, where:

Equation 16

The dimensional variables are defined as follows:

t is time

T is aquifer transmissivity

S is aquifer storativity

L is orthogonal distance from river to well

K’ is aquitard permeability

B’ is aquitard saturated thickness

λ is streambed resistance coefficient

σ is aquitard porosity

Diffuse recharge

Flux to a river from a continuous diffuse recharge source R over land of width L, which is based on solutions in Carslaw and Jaeger (1959; John Knight, personal communications), is given at small times by:

Equation 17

with

Equation 18

At large times,

Equation 19

Equation 20

(non-dimensional time), and the discharge flux D(t) = F(τ) • RL.

Point recharge

Various land use practices such as irrigation and farm dams lead to drainage below the root zone, which eventually recharges the underlying aquifer. This leads to increased discharge to the nearby river, the mechanism is identical to that for overbank flow. The recharge can either be as a step change or variable in time. For the former case, the algorithm for estimating the discharge response for groundwater pumping (Equation 10) is used with the parameter P replaced by the recharge value but with an opposite sign to P (ie. it results in positive discharge flux to the river rather than a negative depletion flux). For the latter case where recharge is time variant, the algorithm for overbank return flow is used with f_inf (of Equation 8) representing deep drainage from irrigation. A number of recharge sources can be modeled with the overall flux estimated using:

Equation 21

where:

CuQ_Ir(t) is the cumulative flux resulting from return flow from a number of irrigation developments [L³/T],

D_Ir(t) is the flux resulting from return flow from an irrigation development [L³/T] (using Equation 10 for constant step recharge, or, Equation 8 for time-variable recharge), and

n_Ir is the number of irrigation developments active along a link.

Evapotranspiration

Evapotranspiration (ET) is a major component of the water budget in vegetated areas that have relatively shallow water tables. In these areas transpiration directly from groundwater by near-shore vegetation can intercept base flow that would otherwise discharge to a stream. Depending on the positioning of the root zone with respect to the water table, the plant can extract water directly from groundwater, from the unsaturated zone, or from both. Actual evapotranspiration remains equal to the potential (PET) value as dictated by the climate down to some depth called the "transition depth" where ET starts to shift from atmospheric-control to soil-moisture control. This implies that below this depth, soil properties and vegetation type would dictate the magnitude of actual ET. At some depth called the "extinction depth", evapotranspiration becomes very low as the vegetation can no longer extract any groundwater.

Figure 5. ET-decline functions for two different land covers

An example of an ET-decline function is shown in Figure 5 and is given by (Shah et al, 2007):

Equation 22

Equation 23

where:

ET_a is the actual evapotranspiration flux [L/T],

PET is the potential evapotranspiration flux [L/T],

d is the depth to groundwater table [L],

b is a decay coefficient [L^-1], and

d′ is the transition depth where ET shifts from atmospheric control to soil-moisture control [L].

The total flow rate due to ET is estimated by:

Equation 24

where:

Q_ET is the ET flow rate [L³/T],

L is the length of the link [L],

b is the lateral extent of the floodplain [L], and

ET_a is actual evapotranspiration flux [L/T] defined by Equations 22 and 23.

Extreme care should be taken to ensure that there is no double accounting for the impacts of any process. For example, when using the flux-based method, the pre-development flux implicitly accounts for ET fluxes. Therefore, one needs only to account for any changes in actual ET that have occurred since that time. That is, if new vegetation has been introduced, the new ET associated with this vegetation should be added. On the other hand, if land has been cleared, or a new groundwater development that captures ET had been introduced since that time, then one needs to deduct that flux from the river loss.

Relevance of processes for saturated/unsaturated connections

The different processes contributing to the exchange flux between a river and the underlying aquifer were discussed in the previous section. However, some of those processes become irrelevant when the connection between the river and underlying aquifer changes from saturated to unsaturated. The flux for unsaturated connection comprises only of the natural within-bank interaction where the flux is driven by the head difference between the river and the aquifer.

Total exchange volume along a link

In the Groundwater module, a linear groundwater system is assumed whereby aquifer transmissivity (T) remains constant in space and time. This means that transient changes in T (eg. as a result of drawdown during pumping) are neglected; this assumption holds as long as the magnitude of the head change is small relative to the saturated aquifer thickness.

The assumption of linearity allows superposition, whereby the impacts of individual stresses can be added to obtain an overall impact on the system. Rassam et al (2004) provided a comprehensive study of the versatility of linearity and superposition. Here, we use the principle of superposition and add the individual flux components resulting from each of the relevant GW-SW interaction processes to estimate the overall exchange flux.

Therefore, the overall flux at any time is estimated by adding the relevant process:

Equation 25

where:

Q_total is the total GW-SW exchange volume along a link during a time-step [L3],

T_s is the length of the time-step [T],

Q_l is the flow rate at low-flow conditions [L³/T],

Q_s is the flow rate due to bank storage [L³/T],

C_uQ_ob is the flow rate due to return from overbank flow [L³/T],

C_uQ_p is the flow rate due to groundwater pumping [L³/T],

C_uQ_Ir is the flow rate due to return from irrigation [L³/T], and

Q_ET is the flow rate due to ET losses [L³/T].

Groundwater salinity

River salinity is assumed to be generated from the groundwater component of flow. The groundwater is assigned a constant salinity (concentration), which may vary along links. This implies that salinity will only be associated with gaining streams (with a saturated connection). The salt load in a link at the end of each time-step is estimated as follows (provided that the net exchange between the groundwater system and river results in a net gain to the river):

Equation 26

where:

S_load is the salt load in the relevant river link [M],

Q_total is the total GW-SW exchange volume along a link during a time-step [L³], and

S_conc is the average concentration of salt in the groundwater system [M/L³).

Methodologies for GW-SW exchange flux estimation

In cases where the exchange flux between a river and the underlying aquifer is known (eg. from a groundwater model), functionality will be made available to import that flux, which precludes the need to estimate the flux. However, in most cases the flux needs to be calculated, using one of the following methods.

Head-based method

From knowledge of the head difference between the river stage and the water table level in the underlying aquifer at any time, the GW-SW exchange flux can be defined on a head basis using Equation 2. The head difference at any time t captures the state of connection between the river and the underlying aquifer and hence represents a realisation of all the stresses in the aquifer up to that particular time t. The adopted sign-convention indicates that a positive Δh refers to a gaining river whereas a negative Δh refers to a losing river. Equation 2 describes a linear relationship between head gradient and flux and thus applies to the straight-line portion of the relationship shown in Figure 3 for a depth to water table of less than D_un.

The advantage of this method is that the flux calculations can start and finish at any time thus coinciding with the simulation (or calibration) period of the river model; this is conceptually correct because of the transient nature of the data. However, the limitation of this method is that the required data is unlikely to be available in many instances.

Flux-based method

Assuming that an initial starting condition can be identified at some time t in the past prior to development (ie. when the direction and magnitude of the exchange flux was known), one can then estimate subsequent fluxes directly from the stresses in the aquifer (eg. a pump) using the appropriate analytical solutions then add them to the initial starting flux. The pre-development flux used should be a long-term average during non-flood event conditions, which is mostly representative of base-flow conditions. Depending on the choice of t and how it relates to when the stresses had started (or will start), the full, or partial impacts of those stresses are superimposed on the initial starting flux to obtain the total exchange flux at any time (including extrapolation into the future). It is a pre-requisite in this method that pre-development initial conditions be identified and used as starting conditions, ie. one must start from a dynamic-equilibrium steady-state state prior to development. This method does not require groundwater head data, which is mostly unavailable. However, there is a considerable effort in identifying all the stresses and collating the relevant data, which are derived from historical records.

Mixed head- and flux-based method

This option is a combination of the previously described methods. Referring back to the flux-based method, one can start at post-development conditions provided that an initial depth to groundwater can be identified and used to calculate an initial starting (transient) flux that accounts for the realisation of all the imposed stresses in the aquifer up to that particular time. Subsequently, the unrealised portion of the impacts from various stresses in the aquifer can then be superimposed on that flux to obtain the cumulative GW-SW exchange flux. Another case where one would need to adopt this method is when pressure heads are being imported from a groundwater model where a key process that leads to significant GW-SW exchanges fluxes might have been neglected such as evapotranspiration (ET) and/or flooding recharge, a very common practice in numerical modelling (eg. using MODFLOW). In those cases, the flux time-series for every (neglected) stress can be calculated and subsequently superimposed on the fluxes imported from the numerical model.

Selecting the appropriate methodology

The process for selecting the appropriate methodology to estimate GW-SW exchange flux is mainly dictated by the type of connection between a river and the underlying aquifer. The final selection of the methodology will be largely impacted by data availability. Figure 6 is a decision-making flowchart that guides the user in making a decision on selecting the appropriate methodology for any situation. The relevant processes for each methodology type and the corresponding required parameters are listed in Table 2. The potential methods for estimating the exchange flux for every connection type are listed in Figure 7 and are discussed in more detail herein.

The most accurate and straightforward method for estimating the flux is based on knowledge of the head difference between the river stage and the watertable in the underlying aquifer at any time. However, in most cases such data may not be available. Alternatively, this data may be derived from a calibrated numerical groundwater model (ie. MODFLOW) that accounts for all the relevant processes active in the field.

Saturated connection, flux-based method

For cases where the head data are not available, a flux-based approach may be adopted. In this approach, one should identify pre-development, steady-state conditions then progressively add the impacts of all stresses that had been acting to date; simulations can also extrapolate into future impacts. This method requires a comprehensive set of data that fully describes all the active stresses in the aquifer along with all the relevant parameters required to model their impact on the aquifer and nearby river.

Saturated connection, mixed head and flux-based method

In some cases, one might opt to adopt a combination of the previous approaches. For example, if groundwater head data is available for a short duration period or has gaps, one can extract a representative head at some time to calculate a transient exchange flux then superimpose subsequent fluxes. The advantage with this method is not having to start from pre-development conditions.

Unsaturated connection

When unsaturated conditions develop beneath a river, some of the processes that contribute to the exchange flux between groundwater and surface water become either less relevant, or irrelevant. Therefore, implementing the flux-based approach based on assessing the impacts of individual stress becomes unreliable.

For deep water tables (depth>D_ml) where the maximum-loss condition has been achieved, the flux becomes a function of the river stage irrespective of variations in groundwater levels. Hence, for such a case knowledge of the magnitude of D_ml is sufficient for calculating flux. When the depth to groundwater is between D_un and D_ml the relationship between flux and head gradient becomes non-linear (Figure 3). For such cases, a time series for groundwater levels would still be ideal, but an average long-term value should be satisfactory as the flux becomes less dependent on depth to water table.

Estimateing GW-SW exchange volume

The GW-SW exchange flux at the current time-step is calculated based on information from the previous time-step (eg. stage height) as follows:

1. Flux time series option:

• The known flux for each time-step is read from a time series.
• The total exchange volume during a time-step is calculated.
• Return to Source with the known exchange volume.

2. Head time series option (known h_wt(t), Equation 2):

• Use the known groundwater head h_wt(t) and the stage height h_r(t) at that time to estimate the flux using Equation 2.
• The total exchange volume during a time-step is calculated.
• Return to Source with the known exchange volume.

3. Link Model option:

a. Assess connectivity, flag (unsaturated/saturated).

b. Unsaturated option:

• • Use an average groundwater head h_wt and the stage height h_r(t) at that time to estimate the flux using Equation 2.
  • The total exchange volume during a time-step is calculated.
  • Return to Source with the known exchange volume.

c. saturated - choose one or more of the following flux components:

i. Natural:

1. at low flow:

• • • Use an average groundwater head h_wt and an average stage height h_r representative of low flow conditions to estimate the flow rate using Equation 2.

2. within bank (bank storage):

• • • use the stage height at two consecutive time-steps to estimate the bank storage flow rate using Equation 3.

3. overbank flow:

• • • flow threshold for overbank flow is examined.
    • algorithm is activated when the flow threshold is exceeded.
    • the discharge flow rate sourced from overbank flow is calculated using Equation 9.

ii. Groundwater pumping:

• • • the depletion flow rate is calculated using Equation 12.
    • the flow may be apportioned across two neighboring links using Equation 14.

iii. Irrigation recharge:

• • • for step change recharge, Equation 10 is used to estimate the flow rate.
    • for time-variable recharge, Equation 11 is used to estimate the flow rate.
    • in both cases the cumulative flow rate from multiple developments is estimated using Equation 21.

iv. ET:

• • • the ET flow rate is estimated using Equation 24.

d. Sum up the relevant flow rate components (a to d):

• • • estimate the total exchange volume during a time-step using Equation 25.
    • return to Source with the known exchange volume.

Calibration process

The algorithms for various types of connections between rivers and aquifers produce a time-series of exchange flux (gain/loss) that is passed to the relevant Source link at every time-step. This flux becomes a new term in the routing formula; some parameters used for estimating the flux must be included in the routing calibration process.

Unsaturated connection

The exchange flux is a function of the head difference and the hydraulic conductance of the river-aquifer interconnection, which is represented mathematically as follows:

Equation 27

where:

Δh(t) is the time-variable head difference between stage and groundwater table, and

C is the hydraulic conductance of the river-aquifer interconnection (m²/day)

The conductance term (C) must be included in routing calibration process.

Saturated connection

The exchange flux for saturated connection comprises the following components:

Natural interaction:

Equation 28

where:

GW_n is the natural exchange flux between the river and the underlying aquifer,

GW_pr is the flux component under pre-development conditions (constant),

GW_wb is within-bank exchange flux (function of head difference between time-variable stage height and a reference (constant) water table height at low-flow conditions Δh(t), and aquifer diffusivity D), and

GW_ob is over-bank exchange flux resulting from return flow from floodplain wetland (function of overbank volume V, potential evaporation PE, aquifer diffusivity D, orthogonal distance between centre of floodplain wetland and river a, and the hydraulic conductivity of the wetland bed k).

Stream depletion

Equation 29

where:

GW_sd is the exchange flux between the river and the underlying aquifer due to groundwater pumping, which is a function of: time-variable pumping rate P(t), aquifer diffusivity D, orthogonal distance between centre of floodplain wetland and river a, and aquifer boundary conditions B.

Irrigation recharge

Equation 30

where:

GW_Ir is the exchange flux between the river and the underlying aquifer due to return flow from irrigation recharge, which is a function of: time-variable deep drainage rate Inf(t), aquifer diffusivity D, and the orthogonal distance between the centre of the irrigation development and river a.

Evapotranspiration

Equation 31

where:

GW_ET is the exchange flux between the river and the underlying aquifer due to evapotranspiration, this is discretised into

GW_ET-pr is a constant predevelopment ET flux

GW_ET-pd is a time-variable (optional) post-development flux,

PET is potential evapotranspiration,

Ex is the extinction depth (below which ET becomes zero),

St is soil type, and

Lc is land cover type.

Total flux

The overall groundwater exchange flux can be grouped into two components for the purpose of calibration, one is time-variant and the other is constant, as follows:

Equation 32

The first constant term comprises the following: (1) the predevelopment natural flux, (2) the predevelopment ET, and (3) the realized impacts of pumping up to the beginning of the calibration period.

The second term comprises the following: (1) the time-variant ET term (which can be dropped as changes in ET are less likely to be modelled), (2) the overbank term representing return flow from wetlands, (3) within-bank exchange flux, and (4) the unrealized impacts the pumping during the calibration period. Therefore, the groundwater flux in the river routing specifications can be written as follows:

Equation 33

where:

GW[t,D] is the time-variable groundwater flux, which is a function of aquifer diffusivity.

GW_{pre-calibration} and D must be included in the routing calibration process.

Salinity

At the end of the routing calibration process, if the refined total flow is a net gain to the river, then a second calibration process for salinity can be carried out whereby the salinity concentration can be varied to achieve the best match with the observed river salinity.

Input data

Input data requirements are summarised in Table 1.

* If head data is derived from groundwater model that has neglected this process, they need to be accounted for.

** As defined in Equation 2 (calculating riverbed conductance requires: K, the hydraulic conductivity of the riverbed sediments; L, the length of the river link; W, the width of the river link; and M, the thickness of the riverbed sediments.

Bibliography

Birkhead, A.L. and James, C.S. (2002). Muskingum river routing with dynamic bank storage. Journal of Hydrology, 264, 113-132.

Glover, R.E. and G.G. Balmer (1954). River depletion resulting from pumping a well near a river, American Geophysical Union Transactions, 35(3), 368-470.

Knight, J., Gilfedder, M., and Walker, G. (2005). Impact of Irrigation and dryland development on groundwater discharge to rivers: A unit response approach to cumulative impacts analysis. Journal of Hydrology, 303(104), 79-91.

Shah, N., Nachabe, M., and Ross, M. (2007). Extinction depth and evapotranspiration from ground water under selected land covers, Ground Water, 45(3), 329-338.

Carslaw, H.S. and Jaeger, J.C. (1959). Conduction of heat in solids, Second edition, Clarendon Press, Oxford.

Hunt, B. (2003). Unsteady stream depletion when pumping from semiconfined aquifer. Journal of Hydrologic Engineering, 8, 12-19.

Rassam, D.W. and Werner (2008). Review of Groundwater - Surface Water Interaction Modelling Approaches and Their Suitability For Australian Conditions, eWater CRC Technical Report .

Rassam, D., Walker, G., and Knight, J. (2004). Applicability of the Unit Response Equation to assess salinity impacts of irrigation development in the Mallee region. CSIRO Technical Report 35/04.