Groundwater interactions (GSWIT)

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.

Note The Groundwater interaction module was previously known as the Groundwater surface water interaction tool (GSWIT).

Figure 59 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.

Source v2.10

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.

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
• 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.

The types of river-aquifer connections are shown in Figure 60. Rivers are commonly grouped as gaining or losing, whereby the former gains water from the aquifer (Figure 60.A) and the latter loses water to the aquifer (Figure 60.B-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 60.A) whereas a losing river might have either a saturated or an unsaturated connection. Figure 60.B shows a losing river with saturated connection and Figure 60.D 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 Dun (Figure 60.C) below the riverbed.

Dun is a function of K1/K2 and B/D, where K1 and K2 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 Dun below the riverbed (Figure 61, depth < Dun).

When the depth to the water table is larger than Dun 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 61, Dun < depth < Dml). When the depth to the water table exceeds Dml, the flux attains a maximum value and becomes independent of depth to water table (Figure 61, depth > Dml 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.

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:

Step 8

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 60); 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

Where:

Ql is the GW-SW exchange flow rate at low flow conditions [L3/T]

hwt is the water table level [L],

hr 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 [L2/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 62 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):

Step 6

Where:

Qs 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 [L2/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 296 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 117

Where:

fev is the evaporative flow rate [L3/T],

Wa is the inundated floodplain area [L2],

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 49

Where:

finf(t) is the infiltrative flow rate [L3/T],

Wa is the inundated floodplain area [L2],

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

Δhw is the head gradient in the wetland

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

Equation 160

Figure 12

Where:

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

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

hw(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 finf(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 Wv becomes zero (the inundated floodplain area empties). The discharge resulting from an instantaneous recharge source finf(t) is given by (Knight et al, 2005):

Figure 46

Where:

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

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

D is aquifer diffusivity [L2/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:

Theory

Where:

CuQob(t) is the cumulative discharge due to overbank return flow [L3/T],

Dob is the discharge due to overbank return flow given by Equation 301,

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

m is the total number of elements in the finf(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 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:

Stability criteria

Where:

Dp 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 [L3/T],

D is aquifer diffusivity [L2/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

Table 51

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:

Table 36

Where:

CuQp(t) is the cumulative depletion flux resulting from a number of pumps [L3/T],

Dp(t) is the depletion flux for an individual pump [L3/T] given by Equation 303, 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:

Table 29

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

Where:

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

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

D is aquifer diffusivity [L2/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):

Figure 29

Where:

α is an integration variable,

λ is the streambed resistance coefficient, and

F(α,t) & 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

,...,...,...

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

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

with

Equation

At large times,

Equation 318

Equation

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

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 303) 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 finf (of Equation 301) representing deep drainage from irrigation. A number of recharge sources can be modeled with the overall flux estimated using:

Equation

Where:

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

DIr(t) is the flux resulting from return flow from an irrigation development [L3/T] (using Equation 303 for constant step recharge, or, Equation 301 for time-variable recharge), and

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

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. An example of an ET-decline function is shown in Figure 63 and is given by (Shah et al, 2007):

Equation

Equation

Where:

ETa 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

Where:

QET is the ET flow rate [L3/T],

L is the length of the link [L],

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

ETa is actual evapotranspiration flux [L/T] defined by Equation 315 & Equation 316.

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.

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.

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

Where:

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

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

Ql is the flow rate at low-flow conditions [L3/T],

Qs is the flow rate due to bank storage [L3/T],

CuQob is the flow rate due to return from overbank flow [L3/T],

CuQp is the flow rate due to groundwater pumping [L3/T],

CuQIr is the flow rate due to return from irrigation [L3/T], and

QET is the flow rate due to ET losses [L3/T].

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

Where:

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

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

Sconc is the average concentration of salt in the groundwater system [M/L3).

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:

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 295. 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 295 describes a linear relationship between head gradient and flux and thus applies to the straight-line portion of the relationship shown in Figure 61 for a depth to watertable of less than Dun.

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.

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.

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.

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 64 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 65 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.

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.

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.

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>Dml) 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 Dml is sufficient for calculating flux. When the depth to groundwater is between Dun and Dml the relationship between flux and head gradient becomes non-linear (Figure 61). 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.

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 hwt(t), Equation 295):

• Use the known groundwater head hwt(t) and the stage height hr(t) at that time to estimate the flux using Equation 295.
• 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 hwt and the stage height hr(t) at that time to estimate the flux using Equation 295.
• 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 hwt and an average stage height hr representative of low flow conditions to estimate the flow rate using Equation 295.

2. within bank (bank storage):

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

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 302.

ii. Groundwater pumping:

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

iii. Irrigation recharge:

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

iv. ET:

• the ET flow rate is estimated using Equation 317.

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

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

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.

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

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 (m2/day)

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

The exchange flux for saturated connection comprises the following components:

Natural interaction:

Equation

Where:

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

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

GWwb 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

GWob 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

Where:

GWsd 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

Where:

GWIr 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

Where:

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

GWET-pr is a constant predevelopment ET flux

GWET-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.

Note that for post-development conditions, PET, Ex, and Lc can be time-variable (reflecting climate change and/or land use change).

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

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

Where:

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

GWpre-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 requirements are summarised in Table 52.

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.

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

** As defined in Equation 295 (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.