Groundwater models
Authors: Andy Wilkins and Maëlle Grass
The PorousFlow module may be used to simulate groundwater systems. A typical example can be found in Herron et al. (2018). Large and complex models may be built that include the effect of rainfall patterns, spatially and temporally varying evaporation and transpiration, seasonal river flows, realistic topography and hydrogeology, as well as human factors such as mines, water bores, dams, etc.
It is also possible to include the vadose zone (unsaturated flow), to track the transport of tracers and pollutants through the system (multi-component flow), to utilize high-precision equations of state for the groundwater (density and viscosity as functions of porepressure and temperature), to include heat flows, to include the impact of rock deformation (solid mechanics), to include gas flows (multi-phase flows) and geochemical reactions. This page concentrates mostly on the construction of traditional groundwater models and only touches lightly on these more elaborate options, which are explained in the other PorousFlow examples, tutorials and tests.
Groundwater concepts
Because PorousFlow has been built as a general multi-component, multi-phase simulator, its models and physics are most conveniently and commonly expressed in terms of quantities that are alien to many experienced groundwater modellers, such as porepressure, permeability and porosity. This section provides translations between more conventional concepts — hydraulic head, hydraulic conductivity and storativity — and the quantities usually used in PorousFlow simulations. In the following sections, PorousFlow models are presented using the conventional groundwater quantities (hydraulic head, etc) as well as the "natural" quantities (porepressure, etc) but for the reasons mentioned below, users are encouraged to use the "natural" quantities in their models.
Hydraulic head and porepressure
Many experienced groundwater modellers will be used to expressing their models in terms of hydraulic head. For instance, if a borehole is drilled into an aquifer that has a hydraulic head of 5m, then the groundwater in the borehole will rise 5m above the water table. (This is assuming the hydraulic head is measured relative to the local water table, which is the common convention.) This is shown in Figure 1.
Figure 1: Hydraulic head, , of an aquifer is measured by recording the height above the water table (the datum in this case) of groundwater in a borehole.
Hydraulic head determines the direction of groundwater flow. For instance, if the hydraulic head of a sandstone aquifer is greater than a siltstone aquifer, groundwater will attempt to flow from the sandstone to the siltstone. Or, if the hydraulic gradient increases from south to north — that is, the hydraulic head is small in the southern regions and is large in the northern regions — then groundwater will tend to flow from north to south.
PorousFlow models are more naturally expressed in terms of porepressure, which is the pressure of the groundwater. This is related to hydraulic head through: (1) Here:
is porepressure, measured in Pa
is hydraulic head, with respect to a fixed datum such as the elevation above sea level of water table at a monitoring site, measured in m
is depth below the fixed datum, measured in m
is the water density (approximately 1000kg.m) and is the acceleration due to gravity (approximately 10m.s), so Pa.m. (Recall 1Pa = 1N.m = 1kg.m.s.)
Consider a confined aquifer with constant hydraulic head, . This means that the groundwater in a borehole will rise metres above the fixed datum (e.g. the local groundwater level), irrespective of where the borehole is drilled and how deep it is (assuming the well intersects this aquifer and not another one). Eq. (1) shows that porepressure increases with depth! This may be conceptually challenging for modellers used to hydraulic head. Remember, porepressure is the actual pressure of the groundwater. It increases with depth because of the weight of the groundwater sitting above it.
Porepressure is advantageous to use in PorousFlow because users will often like to enhance their models using the advanced inbuilt features of PorousFlow, such as
realistic equations of state for the water (or brine), in which density and viscosity are functions of porepressure
vadose-zone physics (unsaturated flow) which is conventionally expressed in terms of porepressure
coupling with solid mechanics through the effective stress, which is expressed in terms of porepressure
Below, models are expressed in terms of hydraulic head and porepressure, but in order to facilitate easy extension of their models, users are encouraged to utilize porepressure.
Hydraulic conductivity and permeability
Imagine a pipe of cross-sectional area filled with porous material, as shown in Figure 2 Assume there is a constant hydraulic gradient between its ends: Darcy's law states that the flow rate of water, , is proportional to the hydraulic gradient: (2) In this equation:
is the flow rate of water through the pipe, measured in m.s
is the cross-sectional area of the pipe, measured in m
is the hydraulic conductivity of the porous material inside the pipe, measured in m.s
is the hydraulic gradient, and the negative sign indicates water flows from regions of high head to low head.
Figure 2: Darcy flow, , through a pipe of cross-sectional area due to a hydraulic gradient.
To provide extra context for the meaning of hydraulic conductivity, imagine a vertical column of saturated porous material that is open at both ends. This could be the rock above an underground mine, for instance. Gravity will cause the water to move downwards (and be replaced by air at the top end). The hydraulic head in this situation is , where is the elevation (and the porepressure is zero, or atmospheric pressure), so . The flow of water per unit area is therefore That is, under gravity-driven flow, the flow-rate of water controlled by , measured in m(water)/m(area)/s.
The physics of PorousFlow is more naturally presented in terms of permeability rather than hydraulic conductivity. Permeability is an intrinsic quantity of a porous material, in contrast to hydraulic conductivity, which only pertains to water. Permeability is therefore useful in multi-phase systems that may not even include water, and in situations where the fluids do not have constant density. Hence, permeability is usually used in the petroleum literature, geothermal scenarios (where the density of water is varying), unconventional-gas models, etc.
Writing Eq. (2) using the porepressure of Eq. (1) yields In these equations:
it is assumed that depth, , does not depend on . If the axis were oriented vertically then the equations would include instead of just .
is the permeability, measured in m.
is the water dynamic viscosity, measured in Pa.s.
The relationship between the permeability and hydraulic conductivity is therefore: (3) where is the density of water and is its dynamic viscosity.
Permeability is often measured in Darcys instead of its natural unit of m. The conversion is
Both permeability and hydraulic conductivity are tensors, so PorousFlow demands that users enter the full tensorial quantity, even though for most cases this is of the form
permeability = 'horizontal_x 0 0 0 horizontal_y 0 0 0 vertical'
Storativity and porosity
A porous material like rock is conceptualised as consisting of solid rock grains and void space. Groundwater flows through the void space. Both storativity and porosity parameterise the amount of groundwater present in a volume of porous material.
Porosity is defined as: Here
is the porosity, which is dimensionless
is the volume of connected void space, measured in m. (There may also be unconnected void spaces, which are ignored in this presentation since they are unimportant for fluid flow and only impact the solid-mechanical deformations.)
is the volume of the porous material (rock) sample, measured in m
is the volume of the solid grains (containing no void space), measured in m
Hence, if the porous material is fully saturated with water at atmospheric pressure, it contains a volume of of water at atmospheric pressure.
However, in the general situation, the volume of water that can be extracted from a volume of porous material is often not because of the following.
The porous material may not be fully saturated with water, meaning that the water volume is less than
Even if fully saturated, the water may be under high pressure, and since it is slightly compressible and the rock grains are slightly compressible, the water volume at atmospheric pressure may be greater than . Similarly, because of the compressibilities, can depend on porepressure, so the that is measured at atmospheric pressure in a lab may be different than the in-situ .
Due to unsaturated physics — namely capillarity and relative permeability — it may be impossible to extract all the water without resorting to extreme pressures, temperatures or chemical stimulants.
All these phenomena are handled automatically by PorousFlow, so porosity is usually employed in PorousFlow models. Porosity is also a natural quantity to use in multi-phase flow scenarios. In most groundwater models, porosity or its cousin, the BiotModulus, are simply specified as fixed quantities for each hydrologic unit. In more elaborate models, porosity is specified at atmospheric pressure and temperature, at zero strain, and in-situ amount of mineralisation, and a PorousFlowPorosity Material handles any changes due to changes in porepressure, temperature, deformation and dissolution/precipitation.
However, for the 3 reasons stated above, storativity, specific yield and specific storage are often used instead of porosity in groundwater models, reports and literature. These are defined as follows.
Storativity or storage coefficient
Here
is the storativity, which is sometimes called the storage coefficient, which is dimensionless
is the specific storage, measured in m
is the aquifer thickness, measured in m
is the specific yield, which is dimensionless
Specific storage
(4) The symbol has been used instead of equality to emphasise that various assumptions have been made concerning the interaction between water pressure and solid-mechanical stresses and strains. In this formula:
is the specific storage, measured in m
is the water density (approximately 1000kg.m) and is the acceleration due to gravity (approximately 10m.s), so Pa.m.
is the compressibility of the rock grains (the porous material without any void space), which is the reciprocal of its bulk modulus, measured in Pa. Typical values are Pa for clay, Pa for sandstone, Pa for shale, Pa for limestone and granite.
is the compressibility of water (approximately Pa), measured in Pa.
Specific yield
Specific yield, , is the volume of water that can be drained by gravity from a volume of porous material. It is therefore close to the porosity, . However, due to unsaturated physics, it is typically a little less than .
Groundwater equations for confined situations
If you dislike equations as much as many groundwater hydrogeologists feel free to read the following note and skip to the next section!
There are many different versions of the groundwater flow equations, and even more numerical implementations of them. Each produces slightly different solutions. The differences are tiny compared with the uncertainties in input parameters such as an aquifer's hydraulic conductivity.
Groundwater models often consider only fully-saturated scenarios such as confined aquifers, and do not consider the vadose zone. This section contains a number of versions of the physical equations governing the groundwater flow in such situations.
These equations are sometimes used to model unsaturated flow (as in unconfined aquifers) however, strictly speaking their physics does not accurately describe such flows. Using the equations below for both saturated and unsaturated zones means your model will be simple and will run rapidly. The porepressure will be positive in fully-saturated regions, and negative in unsaturated regions. If required, you will need to invent your own relationship between saturation and the negative porepressure. Alternately, you may use a PorousFlowUnsaturated Action instead of a PorousFlow BasicTHM to accurately simulate models with both saturated and unsaturated zones, at the cost of providing extra information (the capillary and relative permeability functions) and computational speed.
Presentation 1
Groundwater flow is often modelled by the equation (5) Here:
is time, measured in s
is the porosity, which is dimensionless. It may vary spatially, corresponding to different aquifers and aquitards, for instance.
is the water density, measured in kg.m
is a vector of spatial derivatives, , measured in m. The index is 1, 2 or 3, depending on the spatial direction. If the situation is 2D then the sums in Eq. (5) run over 1 and 2 only.
is the permeability tensor, , measured in m. Usually this is transversely isotropic:
is the viscosity of water, measured in Pa.s
is the porepressure
is the acceleration due to gravity, measured in kg.m.s. If the axis points upwards, this is (i.e., points in the direction of falling).
Presentation 2
Eq. (5) may be extended to the situation in which the porous-material (rock) grains are compressible and may suffer volumetric expansion or compaction due to changes in groundwater porepressure. In this case the equation reads (6) In this equation is the volumetric strain of the porous material.
Presentation 3
As described here, various standard assumptions allow the left-hand side of Eq. (6) to be re-written: (7) In this equation is the so-called Biot modulus, measured in Pa: (8) Here:
is the porosity
is the compressibility of water (approximately Pa), measured in Pa.
is the compressibility of the rock grains (the porous material without any void space), which is the reciprocal of its bulk modulus, measured in Pa.
is the Biot coefficient, which is dimensionless, and describes the impact of the groundwater on solid-mechanical stresses. is usually chosen within the range . If then the groundwater has maximum impact on the solid-mechanical stresses.
In groundwater scenarios, the Biot modulus may be considered time-independent and only dependent on the porous material via material-dependent and .
Presentation 4
Assuming the density of water has very little spatial variation, it may be divided out from Eq. (9) leaving a diffusion equation: (10)
Presentation 5
Notice the similarity of Eq. (4) and Eq. (8). To reasonable accuracy (11) Hence Eq. (9) may be written (12)
Presentation 6
Using the definition of hydraulic head, Eq. (1), and hydraulic conductivity, Eq. (3) in Eq. (12), and making further assumptions based on the almost-incompressibility of water yields the usual groundwater flow equation (13)
Simulating a pumping test
To illustrate the basic setup of a PorousFlow groundwater model, an aquifer pumping test. This example forms part of the PorousFlow QA test suite and is documented from that perspective here.
Place a constant volumetric sink of strength (m.m.s) acting along an infinite line in an isotropic 3D medium. The situation is therefore two dimensional. Assuming that groundwater head obeys Eq. (13), Theis provided the solution for the head (14) which is frequently used in the groundwater literature. "Ei" is the exponential integral function, with values for small (where is the Euler number ), and for large .
A sketch of the modelled situation is shown in Figure 3.
Figure 3: Sketch of a pumping test. A well extracts water at a constant rate from an aquifer with initial porepressure 20MPa, which causes drawdown, as shown by the blue curve.
The MOOSE input file may be broken down into a number of sections, as follows.
The mesh
The situation is two dimensional, but radially symmetric, so can be simulated using a one-dimensional mesh, with a single radial coordinate:
[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
type = GeneratedMesh
dim = 1
nx = 20
xmax = 100
bias_x = 1.05
coord_type = RZ
[]The PorousFlowDictator
All PorousFlow input files must include a PorousFlowDictator. This specifies the number of fluid phases (1 for groundwater models), number of fluid components (1 for groundwater models without tracers) and various other things that enable consistency checking throughout the input file. Virtually all PorousFlow objects require a PorousFlowDictator input parameter. Hence, to shorten input files, virtually all PorousFlow input files contain the PorousFlowDictator information in their GlobalParams block:
[GlobalParams<<<{"href": "../../syntax/GlobalParams/index.html"}>>>]
PorousFlowDictator = dictator
[]This simply means you do not have to type PorousFlowDictator = dictator in all the other PorousFlow objects in your input files.
The physics
In order to replicate Eq. (14), the physics needs to be described by Eq. (13). This may be achieved by employing a PorousFlowBasicTHM Action:
[PorousFlowBasicTHM<<<{"href": "../../syntax/PorousFlowBasicTHM/index.html"}>>>]
dictator_name<<<{"description": "The name of the dictator user object that is created by this Action"}>>> = dictator
add_darcy_aux<<<{"description": "Add AuxVariables that record Darcy velocity"}>>> = false
fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = simple_fluid
gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
multiply_by_density<<<{"description": "If true, then the Kernels for fluid flow are multiplied by the fluid density. If false, this multiplication is not performed, which means the problem linearises, but that care must be taken when using other PorousFlow objects."}>>> = false
porepressure<<<{"description": "The name of the porepressure variable"}>>> = pp
[]At this stage, it is probably a good idea for readers to skim-read the description of the PorousFlowBasicTHM Action, as it will be used extensively in this page.
The Action implements Eq. (9). The Action also implements coupling with solid mechanics and heat flow, which are not of interest in simple groundwater models. To ignore these couplings, simply do not specify the coupling_type input parameter, since it defaults to hydro. Note the flags:
gravity = '0 0 0'
multiply_by_density = false
which mean that Eq. (9) becomes (in 1D, without the and indices) This is of the same form as Eq. (13) if the relationships mentioned above are used (with Pa.m and Pa.s), namely:
Hence, the action is implementing the desired equation, assuming the parameters and are chosen appropriately.
Some additional remarks are:
The Action also implements Eq. (10) when the flag
multiply_by_density = false.The Action also implements Eq. (12) when the Biot Modulus (see below) is set to as in Eq. (11)
The Action does not have a "hydraulic head" input parameter. Instead, it assumes that modellers are using "porepressure", users must supply the name of the porepressure variable using the
porepressure = XXXXinput (hereXXXXis the name of the porepressure variable).In the above block
add_darcy_aux = falseis used simply because the QA test does not involve checking the Darcy-velocity directions. Usually you will want to omit this line and PorousFlow will automatically add Darcy-velocity information to your output.
Variables
In this input file, and most groundwater input files, there is just one Variable (unknown to solve for). It is called pp here, although it may be interpreted as hydraulic head using the scaling :
[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
[pp]
initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 20E6
[]
[]Equation of state for water
PorousFlow requires users to specify their "equation of state". These are functional relationships that provide density and viscosity (and thermal conductivity, enthalpy, etc) as functions of porepressure and temperature. In virtually all groundwater models a SimpleFluidProperties object is sufficient:
[FluidProperties<<<{"href": "../../syntax/FluidProperties/index.html"}>>>]
[simple_fluid]
type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../source/fluidproperties/SimpleFluidProperties.html"}>>>
viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 0.001
[]
[]Here:
the
viscosityis specified exactly, to make the relationship clearer in this example.the name
simple_fluidis provided to the PorousFlowBasicTHM, above
Material properties
The permeability and Biot Modulus need to be specified. This is always done via a Materials block:
[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
[porosity]
type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../source/materials/PorousFlowPorosityConst.html"}>>>
porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material). This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.05
[]
[biot_mod]
type = PorousFlowConstantBiotModulus<<<{"description": "Computes the Biot Modulus, which is assumed to be constant for all time. Sometimes 1 / BiotModulus is called storativity", "href": "../../source/materials/PorousFlowConstantBiotModulus.html"}>>>
fluid_bulk_modulus<<<{"description": "Fluid bulk modulus"}>>> = 2E9
biot_coefficient<<<{"description": "Biot coefficient"}>>> = 1.0
[]
[permeability]
type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-14 0 0 0 1E-14 0 0 0 1E-14'
[]
[]In this case, setting porosity = 0.05, fluid_bulk_modulus = 2E9 and biot_coefficient = 1.0 means that the BiotModulus and the specific storage are The permeability is isotropic:
The produced water flux
This is implemented by a PorousFlowSquarePulsePointSource:
[DiracKernels<<<{"href": "../../syntax/DiracKernels/index.html"}>>>]
[sink]
type = PorousFlowSquarePulsePointSource<<<{"description": "Point source (or sink) that adds (removes) fluid at a constant mass flux rate for times between the specified start and end times.", "href": "../../source/dirackernels/PorousFlowSquarePulsePointSource.html"}>>>
point<<<{"description": "The x,y,z coordinates of the point source (sink)"}>>> = '0 0 0'
mass_flux<<<{"description": "The mass flux at this point in kg/s (positive is flux in, negative is flux out)"}>>> = -0.16E-3 # recall this is a volumetric flux because multiply_by_density = false in the Action, so this corresponds to a mass_flux of 0.16 kg/s/m because density=1000
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
[]
[]In this 2D case with multiply_by_density = false, this means m of water is removed from every 1m of aquifer height every second. If multiply_by_density = true and the model was 2D then the mass_flux in a PorousFlowSquarePulsePointSource would be measured in kg.s.m, and if the model was 3D the units would be kg.s.
The entire input file
The remainder of the input file involves the usual MOOSE objects, such as an Executioner to define the end time and time-stepping, as well as a Preconditioning block, and some output objects. The entire input file reads:
# Theis problem: Flow to single sink using BasicTHM
# SinglePhase
# RZ mesh
[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
type = GeneratedMesh
dim = 1
nx = 20
xmax = 100
bias_x = 1.05
coord_type = RZ
[]
[GlobalParams<<<{"href": "../../syntax/GlobalParams/index.html"}>>>]
PorousFlowDictator = dictator
[]
[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
[pp]
initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 20E6
[]
[]
[PorousFlowBasicTHM<<<{"href": "../../syntax/PorousFlowBasicTHM/index.html"}>>>]
dictator_name<<<{"description": "The name of the dictator user object that is created by this Action"}>>> = dictator
add_darcy_aux<<<{"description": "Add AuxVariables that record Darcy velocity"}>>> = false
fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = simple_fluid
gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 0'
multiply_by_density<<<{"description": "If true, then the Kernels for fluid flow are multiplied by the fluid density. If false, this multiplication is not performed, which means the problem linearises, but that care must be taken when using other PorousFlow objects."}>>> = false
porepressure<<<{"description": "The name of the porepressure variable"}>>> = pp
[]
[FluidProperties<<<{"href": "../../syntax/FluidProperties/index.html"}>>>]
[simple_fluid]
type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../source/fluidproperties/SimpleFluidProperties.html"}>>>
viscosity<<<{"description": "Constant dynamic viscosity (Pa.s)"}>>> = 0.001
[]
[]
[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
[porosity]
type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../source/materials/PorousFlowPorosityConst.html"}>>>
porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material). This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.05
[]
[biot_mod]
type = PorousFlowConstantBiotModulus<<<{"description": "Computes the Biot Modulus, which is assumed to be constant for all time. Sometimes 1 / BiotModulus is called storativity", "href": "../../source/materials/PorousFlowConstantBiotModulus.html"}>>>
fluid_bulk_modulus<<<{"description": "Fluid bulk modulus"}>>> = 2E9
biot_coefficient<<<{"description": "Biot coefficient"}>>> = 1.0
[]
[permeability]
type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-14 0 0 0 1E-14 0 0 0 1E-14'
[]
[]
[DiracKernels<<<{"href": "../../syntax/DiracKernels/index.html"}>>>]
[sink]
type = PorousFlowSquarePulsePointSource<<<{"description": "Point source (or sink) that adds (removes) fluid at a constant mass flux rate for times between the specified start and end times.", "href": "../../source/dirackernels/PorousFlowSquarePulsePointSource.html"}>>>
point<<<{"description": "The x,y,z coordinates of the point source (sink)"}>>> = '0 0 0'
mass_flux<<<{"description": "The mass flux at this point in kg/s (positive is flux in, negative is flux out)"}>>> = -0.16E-3 # recall this is a volumetric flux because multiply_by_density = false in the Action, so this corresponds to a mass_flux of 0.16 kg/s/m because density=1000
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
[]
[]
[VectorPostprocessors<<<{"href": "../../syntax/VectorPostprocessors/index.html"}>>>]
[pp]
type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../source/vectorpostprocessors/LineValueSampler.html"}>>>
num_points<<<{"description": "The number of points to sample along the line"}>>> = 25
start_point<<<{"description": "The beginning of the line"}>>> = '0 0 0'
end_point<<<{"description": "The ending of the line"}>>> = '100 0 0'
sort_by<<<{"description": "What to sort the samples by"}>>> = x
variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = pp
[]
[]
[Preconditioning<<<{"href": "../../syntax/Preconditioning/index.html"}>>>]
[smp]
type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
[]
[]
[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
type = Transient
solve_type = Newton
dt = 200
end_time = 1E3
nl_abs_tol = 1e-10
[]
[Outputs<<<{"href": "../../syntax/Outputs/index.html"}>>>]
perf_graph<<<{"description": "Enable printing of the performance graph to the screen (Console)"}>>> = true
[csv]
type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../source/outputs/CSV.html"}>>>
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = final
[]
[]The results
MOOSE agrees with the Theis solution:
Figure 4: Results of the aquifer pumping test in comparison with the Theis solution.
Abstracting from a groundwater system
In this section, a model of the groundwater system shown in Figure 5 is studied.
Figure 5: The mesh and properties used in the abstraction example.
The system consists of:
An upper aquifer of thickness 10m with porosity, , and isotropic permeability m. With Biot coefficient and water bulk compliance Pa, Eq. (8) implies the Biot modulus is GPa, and Eq. (11) implies the specific storage is m. This number is also produced by Eq. (4) if (Eq. (4) does not contain the Biot coefficient: it approximates ). Using Eq. (3), the hydraulic conductivity is m.day.
An aquitard of thickness 10m with porosity , horizontal permeability m and vertical permeability m. Using the same manipulations as in the preceding paragraph, this yields Biot modulus GPa, specific storage m, horizontal permeability mm.day and mm.day.
A lower aquifer of thickness 10m with the same physical parameters as the upper aquifer.
This system is assumed to exist in isolation from any other aquifers and aquitards. The roof of the upper aquifer is 70m below the water table.
A borehole abstracts water from either the lower aquifer or the aquitard (due to a mistake in the drilling or wellbore completion) and a model is built to investigate the drawdown in the upper aquifer due to these two scenarios.
The finite-element mesh
Creating the mesh is the most difficult aspect of numerical modelling in almost all models, be they groundwater models, structural-mechanical models, models of fluid turbulence, electromagnetic models, etc. For instance, the groundwater model in Herron et al. (2018) took about 3 months to construct since it contained over 100 underground mine workings, many aquifers and aquitards, streams and rivers and alluvium, zones of different evapotranspiration, the ocean, thousands of boreholes, etc. For groundwater modelling, only conceptualisation and parameterisation are more difficult, but these also impact the mesh, for example, there is no need to distinguish between different aquifers if their hydraulic conductivities are the same, up to measurement accuracy.
In this case, the mesh is very simple and may be created using MOOSE's inbuilt commands with the result shown in Figure 5:
[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
[basic_mesh]
type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 3
xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = -50
xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = 50
nx<<<{"description": "Number of elements in the X direction"}>>> = 20
ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = -25
ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = 25
ny<<<{"description": "Number of elements in the Y direction"}>>> = 10
zmin<<<{"description": "Lower Z Coordinate of the generated mesh"}>>> = -100
zmax<<<{"description": "Upper Z Coordinate of the generated mesh"}>>> = -70
nz<<<{"description": "Number of elements in the Z direction"}>>> = 3
[]
[lower_aquifer]
type = SubdomainBoundingBoxGenerator<<<{"description": "Changes the subdomain ID of elements either (XOR) inside or outside the specified box to the specified ID.", "href": "../../source/meshgenerators/SubdomainBoundingBoxGenerator.html"}>>>
input<<<{"description": "The mesh we want to modify"}>>> = basic_mesh
block_id<<<{"description": "Subdomain id to set for inside/outside the bounding box"}>>> = 1
block_name<<<{"description": "Subdomain name to set for inside/outside the bounding box (optional)"}>>> = lower_aquifer
bottom_left<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '-1000 -500 -100'
top_right<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '1000 500 -90'
[]
[aquitard]
type = SubdomainBoundingBoxGenerator<<<{"description": "Changes the subdomain ID of elements either (XOR) inside or outside the specified box to the specified ID.", "href": "../../source/meshgenerators/SubdomainBoundingBoxGenerator.html"}>>>
input<<<{"description": "The mesh we want to modify"}>>> = lower_aquifer
block_id<<<{"description": "Subdomain id to set for inside/outside the bounding box"}>>> = 2
block_name<<<{"description": "Subdomain name to set for inside/outside the bounding box (optional)"}>>> = aquitard
bottom_left<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '-1000 -500 -90'
top_right<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '1000 500 -80'
[]
[upper_aquifer]
type = SubdomainBoundingBoxGenerator<<<{"description": "Changes the subdomain ID of elements either (XOR) inside or outside the specified box to the specified ID.", "href": "../../source/meshgenerators/SubdomainBoundingBoxGenerator.html"}>>>
input<<<{"description": "The mesh we want to modify"}>>> = aquitard
block_id<<<{"description": "Subdomain id to set for inside/outside the bounding box"}>>> = 3
block_name<<<{"description": "Subdomain name to set for inside/outside the bounding box (optional)"}>>> = upper_aquifer
bottom_left<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '-1000 -500 -80'
top_right<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '1000 500 -70'
[]
[]Equation of state for water
PorousFlow requires users to specify their "equation of state". These are functional relationships that provide density and viscosity (and thermal conductivity, enthalpy, etc) as functions of porepressure and temperature. In virtually all groundwater models a SimpleFluidProperties object is sufficient. In this case, the fluid bulk modulus is chosen very large and the thermal expansion coefficient chosen to be zero so that the groundwater density is very close to 1000kg.m so that Pa.m.
[FluidProperties<<<{"href": "../../syntax/FluidProperties/index.html"}>>>]
[simple_fluid]
type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../source/fluidproperties/SimpleFluidProperties.html"}>>>
# the following mean that density = 1000 * exp(P / 1E15) ~ 1000
thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 1E15
[]
[]PorousFlowDictator specification
Most objects in PorousFlow models require the PorousFlowDictator to be specified, so this is most efficiently accomplished through the GlobalParams block:
[GlobalParams<<<{"href": "../../syntax/GlobalParams/index.html"}>>>]
PorousFlowDictator = dictator
[]Variables
The only MOOSE variable (unknown) in this situation is the porepressure.
[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
[pp]
[]
[]The physics
It is assumed that groundwater obeys Eq. (10). This physics is implemented in PorousFlow by the PorousFlowBasicTHM Action with multiply_by_density = false:
[PorousFlowBasicTHM<<<{"href": "../../syntax/PorousFlowBasicTHM/index.html"}>>>]
fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = simple_fluid
gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 -10'
porepressure<<<{"description": "The name of the porepressure variable"}>>> = pp
multiply_by_density<<<{"description": "If true, then the Kernels for fluid flow are multiplied by the fluid density. If false, this multiplication is not performed, which means the problem linearises, but that care must be taken when using other PorousFlow objects."}>>> = false
[]The physics is of the same form as Eq. (13) if the relationships mentioned above are used (with Pa.m and Pa.s), namely:
As mentioned in the section above, the Action does not have a "hydraulic head" input parameter. Instead, it assumes that modellers are using "porepressure", users must supply the name of the porepressure variable using the porepressure = XXXX input (here XXXX is the name of the porepressure variable).
In-situ conditions: hydraulic gradient
Assume that:
That the porepressure is related to hydraulic head via Eq. (1) with
The average hydraulic head in the upper aquifer is 10m and there is a hydraulic gradient of from east to west in the upper aquifer
The hydraulic head in the lower aquifer is 20m. This means the model is not at hydrostatic equilibrium which could be due to recharge into the lower aquifer from a hilly region outside the model
The hydraulic head in the aquitard is a linear combination of these two
The following functions then describe the in-situ conditions:
[upper_aquifer_head]
type = ParsedFunction
expression = '10 + x / 200'
[]
[lower_aquifer_head]
type = ParsedFunction
expression = '20'
[]
[insitu_head]
type = ParsedFunction
symbol_values = 'lower_aquifer_head upper_aquifer_head'
symbol_names = 'low up'
expression = 'if(z <= -90, low, if(z >= -80, up, (up * (z + 90) - low * (z + 80)) / (10.0)))'
[]
[insitu_pp]
type = ParsedFunction
symbol_values = 'insitu_head'
symbol_names = 'h'
expression = '(h - z) * 1E4'
[]Initial and boundary conditions
Assume that the initial conditions are as specified by the insitu_pp Function describe above. Assume also that porepressure is fixed on the model's boundaries. This corresponds to a ready source of water that flows into the model from regions outside the model. The relevant input-file blocks are:
[ICs<<<{"href": "../../syntax/ICs/index.html"}>>>]
[pp]
type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = pp
function<<<{"description": "The initial condition function."}>>> = insitu_pp
[]
[][BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
[pp]
type = FunctionDirichletBC<<<{"description": "Imposes the essential boundary condition $u=g(t,\\vec{x})$, where $g$ is a (possibly) time and space-dependent MOOSE Function.", "href": "../../source/bcs/FunctionDirichletBC.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
function<<<{"description": "The forcing function."}>>> = insitu_pp
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'left right top bottom front back'
[]
[]Prescribing material properties: porosity and permeability
The material properties are prescribed using the following block
[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
[porosity_aquifers]
type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../source/materials/PorousFlowPorosityConst.html"}>>>
porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material). This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.05
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 'upper_aquifer lower_aquifer'
[]
[porosity_aquitard]
type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../source/materials/PorousFlowPorosityConst.html"}>>>
porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material). This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.2
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = aquitard
[]
[biot_mod]
type = PorousFlowConstantBiotModulus<<<{"description": "Computes the Biot Modulus, which is assumed to be constant for all time. Sometimes 1 / BiotModulus is called storativity", "href": "../../source/materials/PorousFlowConstantBiotModulus.html"}>>>
fluid_bulk_modulus<<<{"description": "Fluid bulk modulus"}>>> = 2E9
biot_coefficient<<<{"description": "Biot coefficient"}>>> = 1.0
[]
[permeability_aquifers]
type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 'upper_aquifer lower_aquifer'
[]
[permeability_aquitard]
type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-16 0 0 0 1E-16 0 0 0 1E-17'
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = aquitard
[]
[]Note the aquitard has anisotropic permeability: m while m.
The abstraction well
An abstraction well may be modelled using a PorousFlow sink. The reader is encouraged to study the documentation on the PorousFlow sinks page in order to understand the following block, since boreholes are frequently used for injection or production in groundwater models. A PorousFlowPolyLineSink is used in this example:
[DiracKernels<<<{"href": "../../syntax/DiracKernels/index.html"}>>>]
[sink]
type = PorousFlowPolyLineSink<<<{"description": "Approximates a polyline sink by using a number of point sinks with given weighting whose positions are read from a file. NOTE: if you are using PorousFlowPorosity that depends on volumetric strain, you should set strain_at_nearest_qp=true in your GlobalParams, to ensure the nodal Porosity Material uses the volumetric strain at the Dirac quadpoints, and can therefore be computed", "href": "../../source/dirackernels/PorousFlowPolyLineSink.html"}>>>
SumQuantityUO<<<{"description": "User Object of type=PorousFlowSumQuantity in which to place the total outflow from the line sink for each time step."}>>> = volume_extracted
point_file<<<{"description": "The file containing the coordinates of the points and their weightings that approximate the line sink. The physical meaning of the weightings depend on the scenario, eg, they may be borehole radii. Each line in the file must contain a space-separated weight and coordinate, viz r x y z. For boreholes, the last point in the file is defined as the borehole bottom, where the borehole pressure is bottom_pressure. If your file contains just one point, you must also specify the line_length and line_direction parameters. Note that you will get segementation faults if your points do not lie within your mesh!"}>>> = ex01.bh_lower
line_length<<<{"description": "Line length. Note this is only used if there is only one point in the point_file."}>>> = 10
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
# following produces a flux of 0 m^3(water)/m(borehole length)/s if porepressure = 0, and a flux of 1 m^3/m/s if porepressure = 1E9
p_or_t_vals<<<{"description": "Tuple of pressure (or temperature) values. Must be monotonically increasing."}>>> = '0 1E9'
fluxes<<<{"description": "Tuple of flux values (measured in kg.m^-1.s^-1 if no 'use_*' are employed). These flux values are multiplied by the line-segment length to achieve a flux in kg.s^-1. A piecewise-linear fit is performed to the (p_or_t_vals,flux) pairs to obtain the flux at any arbitrary pressure (or temperature). If a quad-point pressure is less than the first pressure value, the first flux value is used. If quad-point pressure exceeds the final pressure value, the final flux value is used. This flux is OUT of the medium: hence positive values of flux means this will be a SINK, while negative values indicate this flux will be a SOURCE."}>>> = '0 1'
[]
[]1 0 0 -95
The ex01.bh_lower file contains just one line, and the sink sets line_length = 10 and line_direction = '0 0 1' (by default) so this corresponds to a vertical borehole of length 10m. The lines
p_or_t_vals = '0 1E9'
fluxes = '0 1'
mean that the borehole will extract a flux of 0m(water)/m(borehole)/s if the porepressure is zero or lower, and a flux of 1m(water)/m/s if the porepressure is 1GPa or higher. In between these two porepressures, the flux varies linearly. Hence, this will cause porepressure to reduce to close to zero before an equilibrium is found between the borehole extraction rate and the groundwater flux through the rock. If multiply_by_density = true was used in the PorousFlowBasicTHM, then the units of fluxes would be kg(water)/m(borehole)/s.
If a fixed flux was required instead, then fluxes = '1 1' could be used, or another of the PorousFlow sinks could be employed.
The appropriate PorousFlow sink must be used for your model, because they represent significantly different physics, so produce significantly different results. This example uses a PorousFlowPolyLineSink for a borehole, while virtually all real-life models would use a PorousFlowPeacemanBorehole instead. A rule of thumb is to use a PorousFlowPeacemanBorehole for boreholes, and a PorousFlowPolyLineSink for rivers. This is discussed further towards the bottom of this page.
It is common to refine the mesh around boreholes:
a refinement in the plan mesh allows the extent of drawdown cones to be more accurately predicted;
a refinement in the vertical mesh so that the screened section is captured by at least one finite element is recommended as it allows drawdown on vertical sections to be accurately simulated.
Such refinements are not mandatory. They are not performed in this example because only MOOSE's simple inbuilt meshing functionality is used. (Using the PorousFlowPolyLineSink just so happens to be fine in this model because the single Dirac point coincides with a node in the plane.)
Recording the drawdown
AuxVariables and a VectorPostprocessor are used to record the head change throughout the model:
[AuxVariables<<<{"href": "../../syntax/AuxVariables/index.html"}>>>]
[insitu_head]
[]
[head_change]
[]
[][AuxKernels<<<{"href": "../../syntax/AuxKernels/index.html"}>>>]
[insitu_head]
type = FunctionAux<<<{"description": "Auxiliary Kernel that creates and updates a field variable by sampling a function through space and time.", "href": "../../source/auxkernels/FunctionAux.html"}>>>
variable<<<{"description": "The name of the variable that this object applies to"}>>> = insitu_head
function<<<{"description": "The function to use as the value"}>>> = insitu_head
[]
[head_change]
type = ParsedAux<<<{"description": "Sets a field variable value to the evaluation of a parsed expression.", "href": "../../source/auxkernels/ParsedAux.html"}>>>
coupled_variables<<<{"description": "Vector of coupled variable names"}>>> = 'pp insitu_head'
use_xyzt<<<{"description": "Make coordinate (x,y,z) and time (t) variables available in the function expression."}>>> = true
expression<<<{"description": "Parsed function expression to compute"}>>> = 'pp / 1E4 + z - insitu_head'
variable<<<{"description": "The name of the variable that this object applies to"}>>> = head_change
[]
[][VectorPostprocessors<<<{"href": "../../syntax/VectorPostprocessors/index.html"}>>>]
[drawdown]
type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../source/vectorpostprocessors/LineValueSampler.html"}>>>
variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = head_change
start_point<<<{"description": "The beginning of the line"}>>> = '-50 0 -75'
end_point<<<{"description": "The ending of the line"}>>> = '50 0 -75'
num_points<<<{"description": "The number of points to sample along the line"}>>> = 101
sort_by<<<{"description": "What to sort the samples by"}>>> = x
[]
[]The entire input file
For reference, the entire input file is:
# Groundwater extraction example.
# System consists of two confined aquifers separated by an aquitard
# There is a hydraulic gradient in the upper aquifer
# A well extracts water from the lower aquifer, and the impact on the upper aquifer is observed
# In the center of the model, the roof of the upper aquifer sits 70m below the local water table
[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
[basic_mesh]
type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 3
xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = -50
xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = 50
nx<<<{"description": "Number of elements in the X direction"}>>> = 20
ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = -25
ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = 25
ny<<<{"description": "Number of elements in the Y direction"}>>> = 10
zmin<<<{"description": "Lower Z Coordinate of the generated mesh"}>>> = -100
zmax<<<{"description": "Upper Z Coordinate of the generated mesh"}>>> = -70
nz<<<{"description": "Number of elements in the Z direction"}>>> = 3
[]
[lower_aquifer]
type = SubdomainBoundingBoxGenerator<<<{"description": "Changes the subdomain ID of elements either (XOR) inside or outside the specified box to the specified ID.", "href": "../../source/meshgenerators/SubdomainBoundingBoxGenerator.html"}>>>
input<<<{"description": "The mesh we want to modify"}>>> = basic_mesh
block_id<<<{"description": "Subdomain id to set for inside/outside the bounding box"}>>> = 1
block_name<<<{"description": "Subdomain name to set for inside/outside the bounding box (optional)"}>>> = lower_aquifer
bottom_left<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '-1000 -500 -100'
top_right<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '1000 500 -90'
[]
[aquitard]
type = SubdomainBoundingBoxGenerator<<<{"description": "Changes the subdomain ID of elements either (XOR) inside or outside the specified box to the specified ID.", "href": "../../source/meshgenerators/SubdomainBoundingBoxGenerator.html"}>>>
input<<<{"description": "The mesh we want to modify"}>>> = lower_aquifer
block_id<<<{"description": "Subdomain id to set for inside/outside the bounding box"}>>> = 2
block_name<<<{"description": "Subdomain name to set for inside/outside the bounding box (optional)"}>>> = aquitard
bottom_left<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '-1000 -500 -90'
top_right<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '1000 500 -80'
[]
[upper_aquifer]
type = SubdomainBoundingBoxGenerator<<<{"description": "Changes the subdomain ID of elements either (XOR) inside or outside the specified box to the specified ID.", "href": "../../source/meshgenerators/SubdomainBoundingBoxGenerator.html"}>>>
input<<<{"description": "The mesh we want to modify"}>>> = aquitard
block_id<<<{"description": "Subdomain id to set for inside/outside the bounding box"}>>> = 3
block_name<<<{"description": "Subdomain name to set for inside/outside the bounding box (optional)"}>>> = upper_aquifer
bottom_left<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '-1000 -500 -80'
top_right<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '1000 500 -70'
[]
[]
[GlobalParams<<<{"href": "../../syntax/GlobalParams/index.html"}>>>]
PorousFlowDictator = dictator
[]
[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
[pp]
[]
[]
[ICs<<<{"href": "../../syntax/ICs/index.html"}>>>]
[pp]
type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = pp
function<<<{"description": "The initial condition function."}>>> = insitu_pp
[]
[]
[BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
[pp]
type = FunctionDirichletBC<<<{"description": "Imposes the essential boundary condition $u=g(t,\\vec{x})$, where $g$ is a (possibly) time and space-dependent MOOSE Function.", "href": "../../source/bcs/FunctionDirichletBC.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
function<<<{"description": "The forcing function."}>>> = insitu_pp
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'left right top bottom front back'
[]
[]
[Functions<<<{"href": "../../syntax/Functions/index.html"}>>>]
[upper_aquifer_head]
type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../source/functions/MooseParsedFunction.html"}>>>
expression<<<{"description": "The user defined function."}>>> = '10 + x / 200'
[]
[lower_aquifer_head]
type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../source/functions/MooseParsedFunction.html"}>>>
expression<<<{"description": "The user defined function."}>>> = '20'
[]
[insitu_head]
type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../source/functions/MooseParsedFunction.html"}>>>
symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'lower_aquifer_head upper_aquifer_head'
symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'low up'
expression<<<{"description": "The user defined function."}>>> = 'if(z <= -90, low, if(z >= -80, up, (up * (z + 90) - low * (z + 80)) / (10.0)))'
[]
[insitu_pp]
type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../source/functions/MooseParsedFunction.html"}>>>
symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'insitu_head'
symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'h'
expression<<<{"description": "The user defined function."}>>> = '(h - z) * 1E4'
[]
[l_rate]
type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../source/functions/MooseParsedFunction.html"}>>>
symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'm3_produced dt'
symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'm3_produced dt'
expression<<<{"description": "The user defined function."}>>> = '1000 * m3_produced / dt'
[]
[]
[AuxVariables<<<{"href": "../../syntax/AuxVariables/index.html"}>>>]
[insitu_head]
[]
[head_change]
[]
[]
[AuxKernels<<<{"href": "../../syntax/AuxKernels/index.html"}>>>]
[insitu_head]
type = FunctionAux<<<{"description": "Auxiliary Kernel that creates and updates a field variable by sampling a function through space and time.", "href": "../../source/auxkernels/FunctionAux.html"}>>>
variable<<<{"description": "The name of the variable that this object applies to"}>>> = insitu_head
function<<<{"description": "The function to use as the value"}>>> = insitu_head
[]
[head_change]
type = ParsedAux<<<{"description": "Sets a field variable value to the evaluation of a parsed expression.", "href": "../../source/auxkernels/ParsedAux.html"}>>>
coupled_variables<<<{"description": "Vector of coupled variable names"}>>> = 'pp insitu_head'
use_xyzt<<<{"description": "Make coordinate (x,y,z) and time (t) variables available in the function expression."}>>> = true
expression<<<{"description": "Parsed function expression to compute"}>>> = 'pp / 1E4 + z - insitu_head'
variable<<<{"description": "The name of the variable that this object applies to"}>>> = head_change
[]
[]
[Postprocessors<<<{"href": "../../syntax/Postprocessors/index.html"}>>>]
[m3_produced]
type = PorousFlowPlotQuantity<<<{"description": "Extracts the value from the PorousFlowSumQuantity UserObject", "href": "../../source/postprocessors/PorousFlowPlotQuantity.html"}>>>
uo<<<{"description": "PorousFlowSumQuantity user object name that holds the required information"}>>> = volume_extracted
outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = 'none'
[]
[dt]
type = TimestepSize<<<{"description": "Reports the timestep size", "href": "../../source/postprocessors/TimestepSize.html"}>>>
outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = 'none'
[]
[l_per_s]
type = FunctionValuePostprocessor<<<{"description": "Computes the value of a supplied function at a single point (scalable)", "href": "../../source/postprocessors/FunctionValuePostprocessor.html"}>>>
function<<<{"description": "The function which supplies the postprocessor value."}>>> = l_rate
[]
[]
[VectorPostprocessors<<<{"href": "../../syntax/VectorPostprocessors/index.html"}>>>]
[drawdown]
type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../source/vectorpostprocessors/LineValueSampler.html"}>>>
variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = head_change
start_point<<<{"description": "The beginning of the line"}>>> = '-50 0 -75'
end_point<<<{"description": "The ending of the line"}>>> = '50 0 -75'
num_points<<<{"description": "The number of points to sample along the line"}>>> = 101
sort_by<<<{"description": "What to sort the samples by"}>>> = x
[]
[]
[PorousFlowBasicTHM<<<{"href": "../../syntax/PorousFlowBasicTHM/index.html"}>>>]
fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = simple_fluid
gravity<<<{"description": "Gravitational acceleration vector downwards (m/s^2)"}>>> = '0 0 -10'
porepressure<<<{"description": "The name of the porepressure variable"}>>> = pp
multiply_by_density<<<{"description": "If true, then the Kernels for fluid flow are multiplied by the fluid density. If false, this multiplication is not performed, which means the problem linearises, but that care must be taken when using other PorousFlow objects."}>>> = false
[]
[FluidProperties<<<{"href": "../../syntax/FluidProperties/index.html"}>>>]
[simple_fluid]
type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../source/fluidproperties/SimpleFluidProperties.html"}>>>
# the following mean that density = 1000 * exp(P / 1E15) ~ 1000
thermal_expansion<<<{"description": "Constant coefficient of thermal expansion (1/K)"}>>> = 0
bulk_modulus<<<{"description": "Constant bulk modulus (Pa)"}>>> = 1E15
[]
[]
[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
[porosity_aquifers]
type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../source/materials/PorousFlowPorosityConst.html"}>>>
porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material). This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.05
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 'upper_aquifer lower_aquifer'
[]
[porosity_aquitard]
type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../source/materials/PorousFlowPorosityConst.html"}>>>
porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material). This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.2
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = aquitard
[]
[biot_mod]
type = PorousFlowConstantBiotModulus<<<{"description": "Computes the Biot Modulus, which is assumed to be constant for all time. Sometimes 1 / BiotModulus is called storativity", "href": "../../source/materials/PorousFlowConstantBiotModulus.html"}>>>
fluid_bulk_modulus<<<{"description": "Fluid bulk modulus"}>>> = 2E9
biot_coefficient<<<{"description": "Biot coefficient"}>>> = 1.0
[]
[permeability_aquifers]
type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-12 0 0 0 1E-12 0 0 0 1E-12'
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 'upper_aquifer lower_aquifer'
[]
[permeability_aquitard]
type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-16 0 0 0 1E-16 0 0 0 1E-17'
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = aquitard
[]
[]
[DiracKernels<<<{"href": "../../syntax/DiracKernels/index.html"}>>>]
[sink]
type = PorousFlowPolyLineSink<<<{"description": "Approximates a polyline sink by using a number of point sinks with given weighting whose positions are read from a file. NOTE: if you are using PorousFlowPorosity that depends on volumetric strain, you should set strain_at_nearest_qp=true in your GlobalParams, to ensure the nodal Porosity Material uses the volumetric strain at the Dirac quadpoints, and can therefore be computed", "href": "../../source/dirackernels/PorousFlowPolyLineSink.html"}>>>
SumQuantityUO<<<{"description": "User Object of type=PorousFlowSumQuantity in which to place the total outflow from the line sink for each time step."}>>> = volume_extracted
point_file<<<{"description": "The file containing the coordinates of the points and their weightings that approximate the line sink. The physical meaning of the weightings depend on the scenario, eg, they may be borehole radii. Each line in the file must contain a space-separated weight and coordinate, viz r x y z. For boreholes, the last point in the file is defined as the borehole bottom, where the borehole pressure is bottom_pressure. If your file contains just one point, you must also specify the line_length and line_direction parameters. Note that you will get segementation faults if your points do not lie within your mesh!"}>>> = ex01.bh_lower
line_length<<<{"description": "Line length. Note this is only used if there is only one point in the point_file."}>>> = 10
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
# following produces a flux of 0 m^3(water)/m(borehole length)/s if porepressure = 0, and a flux of 1 m^3/m/s if porepressure = 1E9
p_or_t_vals<<<{"description": "Tuple of pressure (or temperature) values. Must be monotonically increasing."}>>> = '0 1E9'
fluxes<<<{"description": "Tuple of flux values (measured in kg.m^-1.s^-1 if no 'use_*' are employed). These flux values are multiplied by the line-segment length to achieve a flux in kg.s^-1. A piecewise-linear fit is performed to the (p_or_t_vals,flux) pairs to obtain the flux at any arbitrary pressure (or temperature). If a quad-point pressure is less than the first pressure value, the first flux value is used. If quad-point pressure exceeds the final pressure value, the final flux value is used. This flux is OUT of the medium: hence positive values of flux means this will be a SINK, while negative values indicate this flux will be a SOURCE."}>>> = '0 1'
[]
[]
[UserObjects<<<{"href": "../../syntax/UserObjects/index.html"}>>>]
[volume_extracted]
type = PorousFlowSumQuantity<<<{"description": "Records total mass flowing into a borehole", "href": "../../source/userobjects/PorousFlowSumQuantity.html"}>>>
[]
[]
[Preconditioning<<<{"href": "../../syntax/Preconditioning/index.html"}>>>]
[smp]
type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
[]
[]
[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
type = Transient
solve_type = Newton
[TimeStepper<<<{"href": "../../syntax/Executioner/TimeStepper/index.html"}>>>]
type = SolutionTimeAdaptiveDT
dt = 1.1E5
[]
end_time = 3.456E5 # 4 days
nl_abs_tol = 1E-13
[]
[Outputs<<<{"href": "../../syntax/Outputs/index.html"}>>>]
[csv]
type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../source/outputs/CSV.html"}>>>
file_base<<<{"description": "The desired solution output name without an extension. If not provided, MOOSE sets it with Outputs/file_base when available. Otherwise, MOOSE uses input file name and this object name for a master input or uses master file_base, the subapp name and this object name for a subapp input to set it."}>>> = ex01_lower_extraction
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = final
[]
[]Results
Water is abstracted from the lower aquifer when point_file = ex01.bh_lower in the PorousFlowPolyLineSink. This production rate is approximately l.s, and the upper aquifer experiences very little drawdown. Water is abstracted from the aquitard when point_file = ex01.bh_aquitard in the PorousFlowPolyLineSink, which results in production rate of about l.s and significant drawdown in the upper aquifer. The drawdown cone reaches steady state after about 4 days, and is shown in Figure 6.
Figure 6: Drawdown in the upper aquifer when water is abstracted from the lower aquifer (blue line) or the aquitard (red line).
The results are mesh-dependent, as expected. A finer mesh produces more drawdown near to the well since it is more able to resolve spatial structure there.
The severe abstraction from the aquitard causes negative hydraulic head in the upper aquifer. This does not mean the porepressure is negative in this case. However, it is easy to build a model in which the porepressure becomes negative. This corresponds physically to desaturation, and to properly model such situations a PorousFlowUnsaturated Action should be used instead of a PorousFlowBasicTHM Action. In many cases, the extra computational cost of modelling unsaturated physics, and the difficulty of appropriately parameterising the model, mean that groundwater modellers choose to ignore partially-saturated physics and simply employ the equivalent of PorousFlowBasicTHM regardless of unsaturated zones. In PorousFlow it is relatively straightforward to swap between PorousFlowBasicTHM and PorousFlowUnsaturated to investigate the impact of including partially-saturated physics.
Baseflow, ET, recharge, unsaturated flow: impact of groundwater abstraction on baseflow to a river
In this example, the impact of groundwater abstraction on baseflow to a river is modelled. The hydrostratigraphy is shown in Figure 7. It consists of:
an upper, unconfined aquitard with topography of varying elevation representing a hill and valley
an aquitard that outcrops in the valley region
a lower aquifer that also outcrops in the valley region
The lower confined aquifer may be recharged by a perennial river that is assumed to exist at the valley's base. Or, it may provide baseflow to the river. Whether leakage or baseflow occurs depends on the groundwater porepressure: the position of the water table with respect to the river stage height.
A horizontal abstraction borehole is drilled into the lower aquifer below the river. The question to answer is: how does the groundwater abstraction influence the groundwater-river interaction? (Such horizontal boreholes are used in gas exploitation, but the main reason it is used here is to illustrate boreholes don't need to be vertical.)
To answer this question, the parameters defining the system will be described in the context of the PorousFlow input file, then the steady-state virgin system will be solved, then the borehole introduced and the impact assessed.
Figure 7: The hydrostratigraphy and the mesh used to explore the impact of groundwater abstraction on baseflow to the river. The unconfined aquifer is shown in brown-yellow, the aquitard in red and the confined aquifer in blue. The river is shown as a bright green line. Vertical exaggeration equals 2 in this figure. The mesh was created by a closed-source script built by CSIRO for these types of problems.
The mesh
The mesh in this case was created by an external program that defined the hydrogeology and the mesh, but not the boundaries that are required for this model. Hence, MOOSE's in-built MeshModifiers are used to identify the appropriate boundaries:
[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
[basic_mesh]
# mesh create by external program: lies within -500<=x<=500 and -200<=y<=200, with varying z
type = FileMeshGenerator<<<{"description": "Read a mesh from a file.", "href": "../../source/meshgenerators/FileMeshGenerator.html"}>>>
file<<<{"description": "The filename to read."}>>> = ex02_mesh.e
[]
[name_blocks]
type = RenameBlockGenerator<<<{"description": "Changes the block IDs and/or block names for a given set of blocks defined by either block ID or block name. The changes are independent of ordering. The merging of blocks is supported.", "href": "../../source/meshgenerators/RenameBlockGenerator.html"}>>>
input<<<{"description": "The mesh we want to modify"}>>> = basic_mesh
old_block<<<{"description": "Elements with these block ID(s)/name(s) will be given the new block information specified in 'new_block'"}>>> = '2 3 4'
new_block<<<{"description": "The new block ID(s)/name(s) to be given by the elements defined in 'old_block'."}>>> = 'bot_aquifer aquitard top_aquifer'
[]
[zmax]
type = SideSetsFromNormalsGenerator<<<{"description": "Adds a new named sideset to the mesh for all faces matching the specified normal.", "href": "../../source/meshgenerators/SideSetsFromNormalsGenerator.html"}>>>
input<<<{"description": "The mesh we want to modify"}>>> = name_blocks
normal_tol<<<{"description": "If normal is supplied then faces are only added if face_normal.normal_hat >= 1 - normal_tol, where normal_hat = normal/|normal|"}>>> = 0.1
new_boundary<<<{"description": "The list of boundary names to create on the supplied subdomain"}>>> = zmax
normals<<<{"description": "A list of normals for which to start painting sidesets"}>>> = '0 0 1'
[]
[xmin_bot_aquifer]
type = ParsedGenerateSideset<<<{"description": "A MeshGenerator that adds element sides to a sideset if the centroid of the side satisfies the `combinatorial_geometry` expression.", "href": "../../source/meshgenerators/ParsedGenerateSideset.html"}>>>
input<<<{"description": "The mesh we want to modify"}>>> = zmax
included_subdomains<<<{"description": "A set of subdomain names or ids whose sides will be included in the new sidesets. A side is only added if the subdomain id of the corresponding element is in this set."}>>> = 2
normal<<<{"description": "If supplied, only faces with normal equal to this, up to normal_tol, will be added to the sidesets specified"}>>> = '-1 0 0'
combinatorial_geometry<<<{"description": "Function expression encoding a combinatorial geometry"}>>> = 'x <= -500.0'
new_sideset_name<<<{"description": "The name of the new sideset"}>>> = xmin_bot_aquifer
[]
[]The physics
In this case, unsaturated flow and saturated flow are used, in order to use appropriate physics for the unconfined aquifer. Hence, instead of a PorousFlowBasicTHM Action, a PorousFlowUnsaturated action is used instead:
[PorousFlowUnsaturated<<<{"href": "../../syntax/PorousFlowUnsaturated/index.html"}>>>]
fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = simple_fluid
porepressure<<<{"description": "The name of the porepressure variable"}>>> = pp
[]Readers are encouraged to read the documentation for the PorousFlowUnsaturated action if they need to model unsaturated flow. Usually, the above block would contain information concerning the relative permeability and the capillary pressure functions, sourced from site experiments, but in this case just the default values are used.
Note the important facet of PorousFlowUnsaturated: it does not have a multiply_by_density flag. Hence, all quantities are measured in mass, not volume. For instance, the recharge quantities are measured in kg.s, not m.s.
Porosity and permeability
The porosity is assumed to be 0.05 everywhere, while the permeabilities are assumed to be:
aquifers: m in the horizontal direction, m in the vertical direction
aquitard: m in the horizontal direction, m in the vertical direction
These properties are specified via:
[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
[porosity_everywhere]
type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../source/materials/PorousFlowPorosityConst.html"}>>>
porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material). This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.05
[]
[permeability_aquifers]
type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 'top_aquifer bot_aquifer'
permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-12 0 0 0 1E-12 0 0 0 1E-13'
[]
[permeability_aquitard]
type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = aquitard
permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-16 0 0 0 1E-16 0 0 0 1E-17'
[]
[]Initial conditions
It is likely that the steady-state solution will be approximately hydrostatic. Therefore, to help PorousFlow converge, the model is initialized with the hydrostatic condition. This is non-trivial here because, for fixed elevation (), the depth varies spatially (with and ). The initial mesh, however, has depth information inserted into it, as the variable cosflow_depth. Hence, it may be extracted using a SolutionUserObject:
[initial_mesh]
type = SolutionUserObject
execute_on = INITIAL
mesh = ex02_mesh.e
timestep = LATEST
system_variables = cosflow_depth
[]and then used in a Function:
[initial_pp]
type = SolutionFunction
scale_factor = 1E4
from_variable = cosflow_depth
solution = initial_mesh
[]before finally being used in an initial condition:
[ICs<<<{"href": "../../syntax/ICs/index.html"}>>>]
[pp]
type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = pp
function<<<{"description": "The initial condition function."}>>> = initial_pp
[]
[]This kind of manipulation — using a SolutionUserObject and a SolutionFunction — is rather common in groundwater modelling.
Boundary conditions
The boundary conditions are the most complicated aspect of this input file.
Rainfall recharge
It is assumed that rainfall recharges the groundwater system at the topography with a constant rate of 0.1mm.day, uniformly over the entire topography. This quantity needs to be converted to a mass flux: Hence the BC:
[rainfall_recharge]
type = PorousFlowSink
boundary = zmax
variable = pp
flux_function = -1E-6 # recharge of 0.1mm/day = 1E-4m3/m2/day = 0.1kg/m2/day ~ 1E-6kg/m2/s
[]In most groundwater models, rainfall recharge would vary spatially and temporally, since different surficial expressions and vegetation cover and seasonal patterns impact rainfall recharge. This is easily incorporated through the flux_function or other parameters to the PorousFlowSink.
Evapotranspiration (ET)
A PorousFlowHalfCubicSink is used to model ET:
[evapotranspiration]
type = PorousFlowHalfCubicSink
boundary = zmax
variable = pp
center = 0.0
cutoff = -5E4 # roots of depth 5m. 5m of water = 5E4 Pa
use_mobility = true
fluid_phase = 0
# Assume pan evaporation of 4mm/day = 4E-3m3/m2/day = 4kg/m2/day ~ 4E-5kg/m2/s
# Assume that if permeability was 1E-10m^2 and water table at topography then ET acts as pan strength
# Because use_mobility = true, then 4E-5 = maximum_flux = max * perm * density / visc = max * 1E-4, so max = 40
max = 40
This block requires some explanation, but the reader is strongly encouraged to read the PorousFlow documentation on the PorousFlowHalfCubicSink and PorousFlow boundary conditions in general for further information.
It is assumed that:
ET is zero if the groundwater head is more than 5m below the topography, corresponding to a maximum root depth of plants of 5m.
ET is maximum if the groundwater table is at, or above, the topography.
It is common to make these types of assumptions, with the rooting depth being dependent on vegetation type and/or surficial geology. In this case, it is assumed that the 5m cutoff is independent of time and does not vary spatially. In the above input-file block, more sophisticated PorousFlow groundwater models would make use of the flux_function and PT_shift and make cutoff be a more complicated function to incorporate spatio-temporal variation. The cubic in the PorousFlowHalfCubicSink interpolates smoothly between the maximum ET and zero, but users may build their own ET functions of porepressure or saturation if data is available, employing a PorousFlowPiecewiseLinearSink.
Furthermore, it is assumed that the pan evaporation is 4mm.day, or, expressed in mass units:
It is assumed that if the surficial geology has permeability m and the water table is at the topography, then ET removes groundwater at this rate. Because of the flag use_mobility = true, the max parameter obeys:
where kg.m is the density of water, and Pa.s is its viscosity. Hence, .
The river
The river is modelled using a PorousFlowPolyLineSink:
[DiracKernels<<<{"href": "../../syntax/DiracKernels/index.html"}>>>]
[river]
type = PorousFlowPolyLineSink<<<{"description": "Approximates a polyline sink by using a number of point sinks with given weighting whose positions are read from a file. NOTE: if you are using PorousFlowPorosity that depends on volumetric strain, you should set strain_at_nearest_qp=true in your GlobalParams, to ensure the nodal Porosity Material uses the volumetric strain at the Dirac quadpoints, and can therefore be computed", "href": "../../source/dirackernels/PorousFlowPolyLineSink.html"}>>>
SumQuantityUO<<<{"description": "User Object of type=PorousFlowSumQuantity in which to place the total outflow from the line sink for each time step."}>>> = baseflow
point_file<<<{"description": "The file containing the coordinates of the points and their weightings that approximate the line sink. The physical meaning of the weightings depend on the scenario, eg, they may be borehole radii. Each line in the file must contain a space-separated weight and coordinate, viz r x y z. For boreholes, the last point in the file is defined as the borehole bottom, where the borehole pressure is bottom_pressure. If your file contains just one point, you must also specify the line_length and line_direction parameters. Note that you will get segementation faults if your points do not lie within your mesh!"}>>> = ex02_river.bh
# Assume a perennial river.
# Assume the river has an incision depth of 1m and a stage height of 1.5m, and these are constant in time and uniform over the whole model. Hence, if groundwater head is 0.5m (5000Pa) there will be no baseflow and leakage.
p_or_t_vals<<<{"description": "Tuple of pressure (or temperature) values. Must be monotonically increasing."}>>> = '-999995000 5000 1000005000'
# Assume the riverbed conductance, k_zz*density*river_segment_length*river_width/riverbed_thickness/viscosity = 1E-6*river_segment_length kg/Pa/s
fluxes<<<{"description": "Tuple of flux values (measured in kg.m^-1.s^-1 if no 'use_*' are employed). These flux values are multiplied by the line-segment length to achieve a flux in kg.s^-1. A piecewise-linear fit is performed to the (p_or_t_vals,flux) pairs to obtain the flux at any arbitrary pressure (or temperature). If a quad-point pressure is less than the first pressure value, the first flux value is used. If quad-point pressure exceeds the final pressure value, the final flux value is used. This flux is OUT of the medium: hence positive values of flux means this will be a SINK, while negative values indicate this flux will be a SOURCE."}>>> = '-1E3 0 1E3'
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
[]
[]This object has already been encountered above, but readers are encouraged to read the detailed documentation on the PorousFlowPolyLineSink and line sources in general.
In this case, it is assumed that the river is perennial and has:
an incision depth of 1m, and
a stage height of 1.5m, and
Therefore, its surface sits 0.5m above the topography. It is assumed that this holds for all simulated times and over the whole river. More sophisticated PorousFlow models can modify the inputs of the PorousFlowPolyLineSink in order to achieve spatio-temporal variation in these parameters. In the case at hand, if the porepressure at the topography is 5kPa (corresponding to 0.5m head) there is zero baseflow or leakage. Assume that the flux satisfies
where kPa. This is implemented by setting
p_or_t_vals = '-Pbig+5000 5000 Pbig+5000'
fluxes = '-M 0 M'
where Pbig is chosen large enough so that it exceeds any conceivable porepressure in the simulation, and . In the case at hand, Pbig = 1E9.
As explained in the sinks page, riverbed conductance, C (measured in kg.Pa.s), is often specified for rivers in groundwater models, and . In the case at hand, using a river width of m and a riverbed thickness of m and a vertical permeability of m, this means , and .
The entire steady-state input file
This reads:
# Steady-state groundwater model. See groundwater_models.md for a detailed description
[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
[basic_mesh]
# mesh create by external program: lies within -500<=x<=500 and -200<=y<=200, with varying z
type = FileMeshGenerator<<<{"description": "Read a mesh from a file.", "href": "../../source/meshgenerators/FileMeshGenerator.html"}>>>
file<<<{"description": "The filename to read."}>>> = ex02_mesh.e
[]
[name_blocks]
type = RenameBlockGenerator<<<{"description": "Changes the block IDs and/or block names for a given set of blocks defined by either block ID or block name. The changes are independent of ordering. The merging of blocks is supported.", "href": "../../source/meshgenerators/RenameBlockGenerator.html"}>>>
input<<<{"description": "The mesh we want to modify"}>>> = basic_mesh
old_block<<<{"description": "Elements with these block ID(s)/name(s) will be given the new block information specified in 'new_block'"}>>> = '2 3 4'
new_block<<<{"description": "The new block ID(s)/name(s) to be given by the elements defined in 'old_block'."}>>> = 'bot_aquifer aquitard top_aquifer'
[]
[zmax]
type = SideSetsFromNormalsGenerator<<<{"description": "Adds a new named sideset to the mesh for all faces matching the specified normal.", "href": "../../source/meshgenerators/SideSetsFromNormalsGenerator.html"}>>>
input<<<{"description": "The mesh we want to modify"}>>> = name_blocks
normal_tol<<<{"description": "If normal is supplied then faces are only added if face_normal.normal_hat >= 1 - normal_tol, where normal_hat = normal/|normal|"}>>> = 0.1
new_boundary<<<{"description": "The list of boundary names to create on the supplied subdomain"}>>> = zmax
normals<<<{"description": "A list of normals for which to start painting sidesets"}>>> = '0 0 1'
[]
[xmin_bot_aquifer]
type = ParsedGenerateSideset<<<{"description": "A MeshGenerator that adds element sides to a sideset if the centroid of the side satisfies the `combinatorial_geometry` expression.", "href": "../../source/meshgenerators/ParsedGenerateSideset.html"}>>>
input<<<{"description": "The mesh we want to modify"}>>> = zmax
included_subdomains<<<{"description": "A set of subdomain names or ids whose sides will be included in the new sidesets. A side is only added if the subdomain id of the corresponding element is in this set."}>>> = 2
normal<<<{"description": "If supplied, only faces with normal equal to this, up to normal_tol, will be added to the sidesets specified"}>>> = '-1 0 0'
combinatorial_geometry<<<{"description": "Function expression encoding a combinatorial geometry"}>>> = 'x <= -500.0'
new_sideset_name<<<{"description": "The name of the new sideset"}>>> = xmin_bot_aquifer
[]
[]
[GlobalParams<<<{"href": "../../syntax/GlobalParams/index.html"}>>>]
PorousFlowDictator = dictator
[]
[Variables<<<{"href": "../../syntax/Variables/index.html"}>>>]
[pp]
[]
[]
[ICs<<<{"href": "../../syntax/ICs/index.html"}>>>]
[pp]
type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = pp
function<<<{"description": "The initial condition function."}>>> = initial_pp
[]
[]
[BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
[rainfall_recharge]
type = PorousFlowSink<<<{"description": "Applies a flux sink to a boundary.", "href": "../../source/bcs/PorousFlowSink.html"}>>>
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = zmax
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
flux_function<<<{"description": "The flux. The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE. The functional form is useful for spatially or temporally varying sinks. Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)"}>>> = -1E-6 # recharge of 0.1mm/day = 1E-4m3/m2/day = 0.1kg/m2/day ~ 1E-6kg/m2/s
[]
[evapotranspiration]
type = PorousFlowHalfCubicSink<<<{"description": "Applies a flux sink to a boundary. The base flux defined by PorousFlowSink is multiplied by a cubic.", "href": "../../source/bcs/PorousFlowHalfCubicSink.html"}>>>
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = zmax
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
center<<<{"description": "Center of the cubic flux multiplier (measured in Pa (or K for temperature BCs))."}>>> = 0.0
cutoff<<<{"description": "Cutoff of the cubic (measured in Pa (or K for temperature BCs)). This needs to be less than zero."}>>> = -5E4 # roots of depth 5m. 5m of water = 5E4 Pa
use_mobility<<<{"description": "If true, then fluxes are multiplied by (density*permeability_nn/viscosity), where the '_nn' indicates the component normal to the boundary. In this case bare_flux is measured in Pa.m^-1. This can be used in conjunction with other use_*"}>>> = true
fluid_phase<<<{"description": "If supplied, then this BC will potentially be a function of fluid pressure, and you can use mass_fraction_component, use_mobility, use_relperm, use_enthalpy and use_energy. If not supplied, then this BC can only be a function of temperature"}>>> = 0
# Assume pan evaporation of 4mm/day = 4E-3m3/m2/day = 4kg/m2/day ~ 4E-5kg/m2/s
# Assume that if permeability was 1E-10m^2 and water table at topography then ET acts as pan strength
# Because use_mobility = true, then 4E-5 = maximum_flux = max * perm * density / visc = max * 1E-4, so max = 40
max<<<{"description": "Maximum of the cubic flux multiplier. Denote x = porepressure - center (or in the case of a heat flux with no fluid, x = temperature - center). Then Flux out is multiplied by (max/cutoff^3)*(2x + cutoff)(x - cutoff)^2 for cutoff < x < 0. Flux out is multiplied by max for x >= 0. Flux out is multiplied by 0 for x <= cutoff."}>>> = 40
[]
[]
[DiracKernels<<<{"href": "../../syntax/DiracKernels/index.html"}>>>]
[river]
type = PorousFlowPolyLineSink<<<{"description": "Approximates a polyline sink by using a number of point sinks with given weighting whose positions are read from a file. NOTE: if you are using PorousFlowPorosity that depends on volumetric strain, you should set strain_at_nearest_qp=true in your GlobalParams, to ensure the nodal Porosity Material uses the volumetric strain at the Dirac quadpoints, and can therefore be computed", "href": "../../source/dirackernels/PorousFlowPolyLineSink.html"}>>>
SumQuantityUO<<<{"description": "User Object of type=PorousFlowSumQuantity in which to place the total outflow from the line sink for each time step."}>>> = baseflow
point_file<<<{"description": "The file containing the coordinates of the points and their weightings that approximate the line sink. The physical meaning of the weightings depend on the scenario, eg, they may be borehole radii. Each line in the file must contain a space-separated weight and coordinate, viz r x y z. For boreholes, the last point in the file is defined as the borehole bottom, where the borehole pressure is bottom_pressure. If your file contains just one point, you must also specify the line_length and line_direction parameters. Note that you will get segementation faults if your points do not lie within your mesh!"}>>> = ex02_river.bh
# Assume a perennial river.
# Assume the river has an incision depth of 1m and a stage height of 1.5m, and these are constant in time and uniform over the whole model. Hence, if groundwater head is 0.5m (5000Pa) there will be no baseflow and leakage.
p_or_t_vals<<<{"description": "Tuple of pressure (or temperature) values. Must be monotonically increasing."}>>> = '-999995000 5000 1000005000'
# Assume the riverbed conductance, k_zz*density*river_segment_length*river_width/riverbed_thickness/viscosity = 1E-6*river_segment_length kg/Pa/s
fluxes<<<{"description": "Tuple of flux values (measured in kg.m^-1.s^-1 if no 'use_*' are employed). These flux values are multiplied by the line-segment length to achieve a flux in kg.s^-1. A piecewise-linear fit is performed to the (p_or_t_vals,flux) pairs to obtain the flux at any arbitrary pressure (or temperature). If a quad-point pressure is less than the first pressure value, the first flux value is used. If quad-point pressure exceeds the final pressure value, the final flux value is used. This flux is OUT of the medium: hence positive values of flux means this will be a SINK, while negative values indicate this flux will be a SOURCE."}>>> = '-1E3 0 1E3'
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pp
[]
[]
[Functions<<<{"href": "../../syntax/Functions/index.html"}>>>]
[initial_pp]
type = SolutionFunction<<<{"description": "Function for reading a solution from file.", "href": "../../source/functions/SolutionFunction.html"}>>>
scale_factor<<<{"description": "Scale factor (a) to be applied to the solution (x): ax+b, where b is the 'add_factor'"}>>> = 1E4
from_variable<<<{"description": "The name of the variable in the file that is to be extracted"}>>> = cosflow_depth
solution<<<{"description": "The SolutionUserObject to extract data from."}>>> = initial_mesh
[]
[baseflow_rate]
type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../source/functions/MooseParsedFunction.html"}>>>
symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'baseflow_kg dt'
symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'baseflow_kg dt'
expression<<<{"description": "The user defined function."}>>> = 'baseflow_kg / dt * 24.0 * 3600.0 / 400.0'
[]
[]
[PorousFlowUnsaturated<<<{"href": "../../syntax/PorousFlowUnsaturated/index.html"}>>>]
fp<<<{"description": "The name of the user object for fluid properties. Only needed if fluid_properties_type = PorousFlowSingleComponentFluid"}>>> = simple_fluid
porepressure<<<{"description": "The name of the porepressure variable"}>>> = pp
[]
[FluidProperties<<<{"href": "../../syntax/FluidProperties/index.html"}>>>]
[simple_fluid]
type = SimpleFluidProperties<<<{"description": "Fluid properties for a simple fluid with a constant bulk density", "href": "../../source/fluidproperties/SimpleFluidProperties.html"}>>>
[]
[]
[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
[porosity_everywhere]
type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../source/materials/PorousFlowPorosityConst.html"}>>>
porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material). This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.05
[]
[permeability_aquifers]
type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 'top_aquifer bot_aquifer'
permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-12 0 0 0 1E-12 0 0 0 1E-13'
[]
[permeability_aquitard]
type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../source/materials/PorousFlowPermeabilityConst.html"}>>>
block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = aquitard
permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1E-16 0 0 0 1E-16 0 0 0 1E-17'
[]
[]
[UserObjects<<<{"href": "../../syntax/UserObjects/index.html"}>>>]
[initial_mesh]
type = SolutionUserObject<<<{"description": "Reads a variable from a mesh in one simulation to another", "href": "../../source/userobjects/SolutionUserObject.html"}>>>
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = INITIAL
mesh<<<{"description": "The name of the mesh file (must be xda/xdr or exodusII file)."}>>> = ex02_mesh.e
timestep<<<{"description": "Index of the single timestep used or \"LATEST\" for the last timestep (exodusII only). If not supplied, time interpolation will occur."}>>> = LATEST
system_variables<<<{"description": "The name of the nodal and elemental variables from the file you want to use for values"}>>> = cosflow_depth
[]
[baseflow]
type = PorousFlowSumQuantity<<<{"description": "Records total mass flowing into a borehole", "href": "../../source/userobjects/PorousFlowSumQuantity.html"}>>>
[]
[]
[Postprocessors<<<{"href": "../../syntax/Postprocessors/index.html"}>>>]
[baseflow_kg]
type = PorousFlowPlotQuantity<<<{"description": "Extracts the value from the PorousFlowSumQuantity UserObject", "href": "../../source/postprocessors/PorousFlowPlotQuantity.html"}>>>
uo<<<{"description": "PorousFlowSumQuantity user object name that holds the required information"}>>> = baseflow
outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = 'none'
[]
[dt]
type = TimestepSize<<<{"description": "Reports the timestep size", "href": "../../source/postprocessors/TimestepSize.html"}>>>
outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = 'none'
[]
[baseflow_l_per_m_per_day]
type = FunctionValuePostprocessor<<<{"description": "Computes the value of a supplied function at a single point (scalable)", "href": "../../source/postprocessors/FunctionValuePostprocessor.html"}>>>
function<<<{"description": "The function which supplies the postprocessor value."}>>> = baseflow_rate
indirect_dependencies<<<{"description": "If the evaluated function depends on other postprocessors they must be listed here to ensure proper dependency resolution"}>>> = 'baseflow_kg dt'
[]
[]
[Preconditioning<<<{"href": "../../syntax/Preconditioning/index.html"}>>>]
[smp]
type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
# following 2 lines are not mandatory, but illustrate a popular preconditioner choice in groundwater models
petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-pc_type -sub_pc_type -pc_asm_overlap'
petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = ' asm ilu 2 '
[]
[]
[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
type = Transient
solve_type = Newton
dt = 1E6
[TimeStepper<<<{"href": "../../syntax/Executioner/TimeStepper/index.html"}>>>]
type = FunctionDT
function = 'max(1E6, t)'
[]
end_time = 1E12
nl_abs_tol = 1E-13
[]
[Outputs<<<{"href": "../../syntax/Outputs/index.html"}>>>]
print_linear_residuals<<<{"description": "Enable printing of linear residuals to the screen (Console)"}>>> = false
[ex]
type = Exodus<<<{"description": "Object for output data in the Exodus format", "href": "../../source/outputs/Exodus.html"}>>>
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = final
[]
[csv]
type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../source/outputs/CSV.html"}>>>
[]
[]Steady-state results
At steady-state, the baseflow to the river is litres(water)/m(river-length)/day, and the saturation is shown in Figure 8.
Figure 8: The saturation at steady-state.
Abstracting groundwater
The only modifications to the steady-state input file needed when modelling the abstraction are the reading of the steady-state solution from the steady-state output result, and the inclusion of the abstraction borehole.
Reading the steady-state solution
The steady-state solution is used to initialise the transient simulation. The mesh is read:
[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
[from_steady_state]
type = FileMeshGenerator<<<{"description": "Read a mesh from a file.", "href": "../../source/meshgenerators/FileMeshGenerator.html"}>>>
file<<<{"description": "The filename to read."}>>> = gold/ex02_steady_state_ex.e
[]
[]A SolutionUserObject is employed to read the steady-state solution:
[steady_state_solution]
type = SolutionUserObject
execute_on = INITIAL
mesh = gold/ex02_steady_state_ex.e
timestep = LATEST
system_variables = pp
[]and the pp Variable is initialized:
[ICs<<<{"href": "../../syntax/ICs/index.html"}>>>]
[pp]
type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = pp
function<<<{"description": "The initial condition function."}>>> = steady_state_pp
[]
[]The abstraction borehole
A PorousFlowPeacemanBorehole is used to model the borehole:
[horizontal_borehole]
type = PorousFlowPeacemanBorehole
SumQuantityUO = abstraction
bottom_p_or_t = -1E5
unit_weight = '0 0 -1E4'
character = 1.0
point_file = ex02.bh
variable = pp
[]with point file representing a long horizontal borehole with radius 0.2m at an elevation of -90m that sits almost directly underneath the river (the river sits at , while the borehole is at ):
0.2 240 -200 -90
0.2 240 -190 -90
0.2 240 -180 -90
0.2 240 -170 -90
0.2 240 -160 -90
0.2 240 -150 -90
0.2 240 -140 -90
0.2 240 -130 -90
0.2 240 -120 -90
0.2 240 -110 -90
0.2 240 -100 -90
0.2 240 -90 -90
0.2 240 -80 -90
0.2 240 -70 -90
0.2 240 -60 -90
0.2 240 -50 -90
0.2 240 -40 -90
0.2 240 -30 -90
0.2 240 -20 -90
0.2 240 -10 -90
0.2 240 0 -90
0.2 240 10 -90
0.2 240 20 -90
0.2 240 30 -90
0.2 240 40 -90
0.2 240 50 -90
0.2 240 60 -90
0.2 240 70 -90
0.2 240 80 -90
0.2 240 90 -90
0.2 240 100 -90
0.2 240 110 -90
0.2 240 120 -90
0.2 240 130 -90
0.2 240 140 -90
0.2 240 150 -90
0.2 240 160 -90
0.2 240 170 -90
0.2 240 180 -90
0.2 240 190 -90
0.2 240 200 -90
The value bottom_p_or_t = -1E5 means that the pump at the borehole's bottom point is able to suck at -1atmosphere of pressure. (In this case, the borehole is horizontal, so its "bottom point" is meaningless, but in general it is the first point in the point_file, as explained in the documentation.)
It is vital that readers understand that a PorousFlowPeacemanBorehole uses an approximation of the groundwater solution around a small borehole in order to accurately represent the small borehole in a coarse mesh. The small borehole need not coincide with the finite-element nodes or elements, since the Peaceman formulation automatically takes the mesh into consideration. In this case, it means that the porepressure at any node around the borehole will not be Pa, even though bottom_p_or_t = -1E5. This is in stark contrast to a PorousFlowPolyLineSink, which can easily force the porepressure to be a fixed value at the borehole if the magnitude of fluxes is set high enough. The difference is quantified in a section below. Almost always:
Peaceman boreholes should be used to represent boreholes
PolyLine sinks should be used to represent surface features such as streams and the interaction with the atmosphere, where the stream or atmosphere is sufficiently "huge" that it does actually fix porepressure at the point of interest.
Results
The borehole extracts very little groundwater — only a paltry 7.6litre/day — which is not surprising given its small diameter and the "tightness" of the aquifer. The porepressure in the proximity of the river is only reduced by around 0.2Pa, as shown in Figure 9. This reduces the baseflow from 11.09litre/m/day to 11.07litre/m/day.
Figure 9: The porepressure change due to groundwater abstraction.
Groundwater modelling therefore suggests that the borehole will have negligible impact on the environment, but will also produce virtually no water.
Using a PolyLineSink
It is a mistake to use a PorousFlowPolyLineSink to represent the borehole in this case, such as
[polyline_sink_borehole]
type = PorousFlowPolyLineSink
SumQuantityUO = abstraction
fluxes = '-0.4 0 0.4'
p_or_t_vals = '-1E8 0 1E8'
point_file = ex02.bh
variable = pp
If the fluxes parameters were chosen sufficiently large, this would fix the porepressure at finite-element nodes surrounding the borehole to the value of 0Pa. Since the finite-element mesh has a resolution of around , this could correspond to a colossal borehole with cross-sectional size (and length of 400m from the width of the model)!
In this case, the water extracted by the borehole is around 640litre/day, and the river becomes leaky with 4.5litre/m/day of river water entering the groundwater system (and flowing to the borehole). The porepressure change is much larger, as shown in Figure 10.
Figure 10: The porepressure change due to groundwater abstraction when erroneously using a PorousFlowPolyLineSink to represent the borehole.
Final words
This page has presented some very simple groundwater models that you can use as a framework to build your realistic models. PorousFlow is highly functional and has been used to solve very complicated groundwater models, but if you require extra functionality not already included, please ask on the MOOSE discussion list.
References
- N. F. Herron, L. Peeters, R. Crosbie, S. P. Marvanek, D. Pagendam, A. Ramage, P. K. Rachakonda, and A. Wilkins.
Groundwater numerical modelling for the Hunter subregion. Product 2.6.2 for the Hunter subregion from the Northern Sydney Basin Bioregional Assessment.
Technical Report, Department of the Environment and Energy, Bureau of Meteorology, CSIRO and Geoscience Australia, Canberra, 2018.
URL: http://data.bioregionalassessments.gov.au/product/NSB/HUN/2.6.2.[BibTeX]