|Mackay, J D, Jackson, C R, and Wang, L. 2014. AquiMod User Manual (v1.0). Nottingham, UK, British geological Survey. (OR/14/007).|
AquiMod has three modules that represent the soil, unsaturated and saturated zone hydrology. Each module has been designed so that it can incorporate a number of possible structures (components), each based on a different conceptual representation of the process being considered, but all of which adhere to the same generalised structure outlined in section 1.2. Each component uses one or more parameters which can be modified to change the behaviour of the model. These parameters and the mathematical algorithms that AquiMod employs are described below. Please note also, that all component parameters and recommended calibration ranges are listed in Appendix 1.
Soil Zone Module
Currently there is only one soil zone component available in the AquiMod software (Table 2).
|1||FAO||Drainage from the base of the soil zone is calculated using a soil water balance method based on a simplification of the algorithm developed by the UN Food and Agricultural Organisation (Allen et al., 1998).|
This component simulates soil moisture as a function of vegetation and soil properties. The soil column is conceptualised as a bucket with a maximum volume of water available to plants after the soil has drained to its field capacity (Figure 3.1a). This is termed the total available water (TAW) which is calculated as:
where Zr is the estimated maximum root depth [L] of the vegetation, and FC and WP are the soil field capacity and wilting point, respectively.
As the soil moisture content decreases, it becomes more difficult for vegetation to extract moisture from the soil matrix. The proportion of TAW that can easily be extracted before this point is reached is conceptualised as readily avilable water (RAW), which is calculated as:
where p is the depletion factor of the vegetation.
The water balance of the soil zone is a function of the rainfall input and evaporative flux from the soil and can be written as:
where SMD is the soil moisture deficit [L], PPTN is the total precipitation input [L] and AET is the actual evapotranspiration [L], which is calculated as a function of potential evapotranspiration (PET) [L] and the SMD (Figure 3.1b):
When the soil zone becomes saturated any rainfall in excess of the evaporative demand of the vegetation becomes excess water (EXW). This excess water is then split between soil drainage and surface runoff. The proportion of the excess water (EXW) that drains to the unsaturated zone is:
where SD is the soil drainage and BFI is the baseflow index, which defines the average proportion of stream flow that a river receives from groundwater discharge. Surface runoff is then calculated as the remainder of the excess water.
Unsaturated Zone Module
Currently there is only one unsaturated zone component available in the AquiMod software (Table 3.2).
|1||Weibull||Drainage from the soil is attenuated through the unsaturated zone using a Weibull distribution transfer function.|
Recharge is distributed over a number of time-steps, n, and the proportion of soil drainage for each time-step is calculated using a two-parameter Weibull probability density function:
where k > 0 is the shape parameter and λ > 0 is the scale parameter of the distribution. The λ parameter primarily controls the location of the peak in the probability density function while k controls the density of the function around the peak (Figure 3.2). The resulting distribution is scaled such that the discrete integral of f is equal to unity and consequently the recharge for each time-step (Rt) is spread over the selected number of time-steps, n:
where α is the scaling parameter. The Weibull function can represent exponentially increasing, exponentially decreasing, and positively and negatively skewed distributions. It is used because it allows the exploration of different distributions, whilst being smooth, which is considered to be more physically justifiable than randomly selected monthly weights.
Saturated Zone Module
There are currently five different saturated zone components available in the AquiMod software as summarised in Table 3.3.
|1||Q3K3S1||A three-layer aquifer representation, each of variable thickness and permeability.|
|2||Q2K2S1||A two-layer aquifer representation, each of variable thickness and permeability.|
|3||Q1K1S1||A single layer aquifer representation.|
|4||Q1T1S1||A confined aquifer representation with a fixed transmissivity.|
|5||VKD||An aquifer with a ‘cocktail glass’ representation of the change in hydraulic conductivity with depth.|
This component represents flow in the saturated zone by a rectangular block of aquifer with dimensions ∆x and ∆y denoting its length and width [L] respectively. A mass balance calculation is performed at each time-step to calculate the new groundwater head:
where R is recharge input [L T-1], Q is the total groundwater discharge [L3 T-1], A [L3 T-1] is any additional abstraction, S is the storage coefficient (dimensionless), and dh is the change in groundwater head [L] over time, dt [T].
This component accommodates three layers of variable thickness and permeability (Figure 3.3). Each layer is independent and has its own discharge outlet. The deepest layer represents groundwater which flows out of the catchment via subsurface flow paths. The two upper outlets are lumped representations of surface discharge points in the catchment including rivers and springs. Both may flow intermittently if the groundwater head falls below the outlet elevation.
The total groundwater discharge is the sum of discharge from all layers in the saturated zone which is calculated using an equation of the form:
where i is the layer number for m layers and ∆hi [L] is the difference between the groundwater head and the elevation of the aquifer outlet point. Importantly due to the explicit form of Equation 9 used, the groundwater head at the previous time-step, h* [L] is used:
where zi is the outlet elevation. Transmissivity, Ti [L2 T-1] is a function of the hydraulic conductivity [L T-1], ki, and is calculated using the following piece-wise function:
Using this component, the hydraulic conductivity can be configured to increase or decrease with depth according to the user’s preference. However, for some aquifers it may be preferable to have an increasing hydraulic conductivity with elevation. As such, an additional α parameter also exists for this component which can be switched to force this condition when using randomly sampled Monte Carlo parameter sets (See Appendix 1).
This component employs the same structure according to equations 8-11, however, here only two layers are defined. Each layer is independent and has its own discharge outlet (Figure 3.4). The deepest layer represents groundwater which flows out of the catchment via subsurface flow paths. The upper outlet is a lumped representation of surface discharge points in the catchment including rivers and springs. It may flow intermittently if the groundwater head falls below the outlet elevation. As with the Q3K3S1 component, the hydraulic conductivity can be forced to increase with elevation when using randomly sampled parameter sets by specifying the α parameter.
This component employs the same structure according to equations 8-11, however, here only one layer is defined with a single discharge outlet (Figure 3.5).
Like the Q1K1S1 saturated zone component, this component also employs a single outlet (Figure 3.6). However, the transmissivity is assumed to be fixed and independent of the hydraulic head, h. This is equivalent to a confined aquifer representation.
As an alternative to using specified layers and outlets, the user can also implement a variable conductivity with depth (VKD) profile. Here a ‘cocktail glass’ representation is can be employed which has two distinct sections (Figure 3.7). In the lower section, between the base of the aquifer, z1 and the elevation of the point of inflection, zp, the hydraulic conductivity is constant. In the upper section above zp, the hydraulic conductivity increases linearly with elevation. Subsequently, the transmissivity is calculated as:
where K1 is the hydraulic conductivity at the base of the aquifer and m is the gradient of the hydraulic conductivity profile above the point of inflection. The total discharge is then calculated as before:
- Allen, R G, Pereira, L S, Raes, D & Smith, M. 1998. Crop evapotranspiration — Guidelines for computing crop water requirements — FAO Irrigation and drainage paper 56. Food and Agriculture Organization of the United Nations.