Navigation:  Numerical Modeling Workflow - Finite Difference Grids > Define Properties >

Background on Flow and Transport Properties

 

The following sections present an overview of the property zone parameters required for flow and transport models in Visual MODFLOW Flex.

A steady-state groundwater flow model requires conductivity and initial heads property values for each active grid cell in order to run a flow simulation, while transient flow models also require storage properties for each active grid cell. Similarly, a transport model requires various model and species parameter values for each active grid cell in order to run a transport simulation. Upon creating a Visual MODFLOW Flex project, the default flow and transport parameter values are assigned to every grid cell in the model domain. This will ensure the model has the minimum data required to run a simulation. However, in most situations, the flow and transport properties will not be uniform throughout the entire model domain, and it will be necessary to assign different property values to different areas of the model.

Heterogeneous property values are supported by Visual MODFLOW Flex using either Constant Value Property Zones, or Distributed Value Property Zones. These two different approaches are described below.

 

Constant Value Property Zones

The Constant Value Property Zones approach is the most simple and straight-forward, and can be used for all model properties supported by Visual MODFLOW Flex. Different model properties are accommodated by grouping grid cells sharing the same property values into “property zones”. Each property zone will (normally) contain a unique set of property values, and is represented by a different grid cell color.

The Constant Value Property Zones approach requires the development of a conceptual model, whereby each hydrostratigraphic unit of the model is assigned a uniform set of property values. For example, consider an aquifer where there is pumping test data and slug test data indicating a range of horizontal conductivity values from 1x10-4 cm/s to 5x10-4 cm/s at different locations within the aquifer. One conceptual approach would assign a uniform Kx and Ky value of 2.5x10-4 cm/s to the entire aquifer. This value would be adjusted up or down for calibration purposes within the range of values reported. If a reasonable calibration cannot be achieved using this conceptual model, it may be necessary to sub-divide this region into several zones to accommodate local irregularities in the flow pattern. However, almost all modeling textbooks strongly recommend starting out simple at first, getting as close a solution as possible, and then making the model more complex, if necessary.

 

Distributed Value Property Zones

The Distributed Value Property Zones approach is available for flow properties (conductivity, storage, and initial Heads), and most transport properties (e.g. initial concentrations, longitudinal dispersivity, sorption parameters). This approach is a little more complicated than the constant value property zone approach because it involves linking a property zone to one or more parameter distribution arrays containing data interpolated from scattered observation points. When a property zone is linked to a distribution array, the property values assigned to each grid cell within that zone are calculated by multiplying the zone parameter value with the corresponding value from the parameter distribution array. If the grid spacing from the model does not match the grid spacing from the distribution array, a bivariate interpolation scheme is used to calculate the appropriate parameter value at the center of the model grid cell using the four nearest data nodes in the parameter distribution array.

 

Flow Properties

 

Conductivity

Hydraulic conductivity is a property of both the porous media and the groundwater (or other fluid) passing through the pore space due to a pressure gradient based on Darcy's Law, an empirical relationship that relates the fluid flux to hydraulic gradient, with the conductivity as a constant of proportionality. In general:

where:

is the groundwater flux [L/T]

is the hydraulic conductivity [L/T]

i is the hydraulic gradient, which is the change in pressure head over a distance , given by over infinitesimal distances or by over finite distances

 

It should be noted that Darcy's Law is limited to conditions where flow is laminar (slow and viscous), which is applicable in most groundwater systems.

Darcy's Law can be extended to three dimensional cases as follows:

where:

and are the X, Y, or Z principal directions

is the hydraulic conductivity tensor with individual terms for flow in the direction due to a gradient in the direction

 

In most versions of MODFLOW and as currently implemented in Flex, the off-diagonal terms of the tensor are assumed to be negligible and the equation reduces to:

 

The conductivity parameters may be defined on a cell-by-cell basis using constant property values and/or distributed property values. When importing or assigning the conductivity property zones, Visual MODFLOW Flex will require valid data for Kx, Ky, and Kz.

Kx is defined as the hydraulic conductivity along the X-axis direction

Ky is defined as the hydraulic conductivity along the Y-axis direction

Kz is defined as the hydraulic conductivity along the Z-axis direction

 

Anisotropy

The reason Visual MODFLOW Flex prompts for both Kx and Ky is because there are two options for defining the horizontal anisotropy of the Conductivity property zones:

Anisotropy by layer

Anisotropy as specified

 

Please Note: The anisotropy option is set in Translation settings. (see Anisotropy for more details).

 

If the 'Anisotropy by layer' option is used, the Kx value will determine the conductivity in the X-direction, and the specified horizontal anisotropy ratio (Ky/Kx) for each layer will be used to calculate the Ky value for each grid cell. Dy default, the horizontal anisotropy is set to 1.

If the 'Anisotropy as specified' option is used, the model will use the Kx and Ky values defined for each property zone.

Storage

Storage parameters are used in transient simulations, particle tracking, and solute transport:

Specific Storage (Ss) is defined as the volume of water that a unit volume of aquifer releases from storage under a unit decline in hydraulic head due to aquifer compaction and water expansion. Visual MODFLOW Flex determines the primary storage coefficient (sf1) for MODFLOW based on the user defined Specific Storage parameter. The primary storage coefficient is calculated by Visual MODFLOW Flex to be equal to the specific storage multiplied by the layer thickness (Specific Storage x thickness = Primary storage coefficient). Specific Storage is is only used in transient simulations.

Specific Yield (Sy) is known as the storage term for an unconfined aquifer. It is defined as the volume of water that an unconfined aquifer releases from storage per unit surface area per unit decline in the water table. For sand and gravel aquifers, specific yield is generally equal to the porosity. MODFLOW uses Ss or Sy depending on the layer type assigned by the user (please refer to "Layer Type Settings"). For an unconfined layer, MODFLOW uses Sy to determine storage volumes. For a confined layer, Ss is used. For a variable layer, MODFLOW will check the head value of the cell to determine if it is confined or not. If you do not have measured parameter values for Sc and Sy it is recommended that you refer to literature values as a default. Specific Yield is only used in transient simulations.

Effective Porosity (Eff. Por) is the pore space through which flow actually occurs, and is used by the particle tracking algorithms (MODPATH/MOD-PATH3DU) to determine the average linear groundwater velocities for use in time-dependent capture zones and time markers along pathlines. This term is not used for MODFLOW simulations.

Total Porosity (Tot. Por) is the percentage of the rock or soil that is void of material, and is used by MT3D/RT3D/SEAWAT (by default) to determine the chemical reaction coefficients, and for calculating the average linear groundwater flow velocity in the particle tracking solution schemes. A different porosity is used for MT3D/RT3D/SEAWAT than for particle tracking because the transport engines account for additional transport and reactive processes, such as dispersion. The total porosity term is not used for MODFLOW simulations.

These storage parameters may be defined on a cell-by-cell basis using constant property values and/or distributed property values. When importing or assigning the storage property zones, Visual MODFLOW Flex will require valid data for each of the above parameters.

Please Note: MODFLOW-SURFACT supports the time-varying material properties (TMP1) package (as an add-on) where flow properties such as conductivities and storage parameters may change over time, as discussed here.

 

Initial Heads

In order to start solving the flow simulation, flow engines require an initial “guess” for the head values in the model. A good initial guess for the starting heads of the simulation can reduce the required run time significantly, while a poor initial guess can increase run times or lead to stability and/or convergence issues. The Initial Head values are also used to calculate drawdown values, as measured by the difference between the starting head and the calculated head.

 

Vadose Zone

The vadose zone object is only available when an unsaturated flow model (MODFLOW-2005 with the unsaturated zone flow (UZF) package or MODFLOW-SURFACT) is selected in the Modeling Objectives workflow step.

The vadose zone, or unsaturated zone, extends from the land surface to the water table. The vadose zone typically has a lower hydraulic conductivity because some of the pore space is filled with air, and the soil moisture in the vadose zone only travels through the wetted cross-section of the pore space. The relative proportion of air to water in the pores can vary, and consequently the hydraulic properties of the porous media can also vary.

In general, soil moisture in the vadose zone is under tension, and it travels through the vadose zone due to total soil-moisture potential (the sum of the elevation potential, and the pressure head). The pressure head is a function of the volumetric water content of the soil, and depends upon whether the soil has previously undergone wetting or drying. The relationship between pressure head and volumetric water content for a particular soil is known as a soil-water characteristic curve. Using the experimental soil-water characteristic curve and the soil hydrologic parameters, the movement of water in the vadose zone can be determined using the BCF4 package in MODFLOW-SURFACT or the UZF package in MODFLOW 2005. A sample index of soil properties is given below:

 

Table: Sample Index of Soil Hydraulic Properties by Soil Texture

Soil Type

Sand
[%]

Clay
[%]

Bulk Density

[g/cm3]

Field Capacity

[cm3/cm3]*

Wilting Point

[cm3/cm3]**

Porosity Fraction

[-]

Saturated Hydraulic Conductivity

[cm/s]

Slope of Retention Curve
[log space] ***

Van Genuchten
Alpha parameter

[cm-1] ****

sand (s)

94.83

2.27

1.49

0.08

0.03

0.43

0.0107

4.1

0.145

loamy sand (ls)

85.23

6.53

1.52

0.15

0.06

0.42

0.003

3.99

0.124

sandy loam (sl)

69.28

12.48

1.57

0.21

0.09

0.4

0.0015

4.84

0.075

silt loam (sil)

19.28

17.11

1.42

0.32

0.12

0.46

0.0011

3.79

0.02

silt (si)

4.5

8.3

1.28

0.28

0.08

0.52

0.0023

3.05

0.016

loam (l)

41

20.69

1.49

0.29

0.14

0.43

0.0005

5.3

0.036

sandy clay loam (scl)

60.97

26.33

1.6

0.27

0.17

0.39

0.0007

8.66

0.059

silty clay loam (sicl)

9.04

33.05

1.38

0.36

0.21

0.48

0.0015

7.48

0.01

clay loam (cl)

30.08

33.46

1.43

0.34

0.21

0.46

0.0005

8.02

0.019

sandy clay (sc)

50.32

39.3

1.57

0.31

0.23

0.41

0.0003

13

0.027

silty clay (sic)

8.18

44.58

1.35

0.37

0.25

0.49

0.0008

9.76

0.005

clay (c)

24.71

52.46

1.39

0.36

0.27

0.47

0.0009

12.28

 

Source: "Average hydraulic conductivity properties of ARS soil texture classes", draft dated February, 2000 by J. Schaake. This expanded the work of others and included a total of 2128 soil samples.

 

* Field capacity is the fractional water content at 1/3 bar tension.

 

** Wilting point is the fractional water content at 15 bar tension;

 

*** Used in Campbell’s equation. ref: Coby et al. (1984)

 

**** Average values of the van Genuchten soil parameters obtained by experimental means (Carsel and Parrish, 1988) and also reported on page 180 in Contaminant Hydrogeology book by C. W. Fetter (1999). You may also want to consider the ROSETTA code for estimating the unsaturated soil hydraulic properties which can be downloaded from the US Salinity Laboratory website: https://data.nal.usda.gov/dataset/rosetta.

 

 

Figure: Soil texture classes triangle.

 

Source: Agriculture and Food Canada, 2017 (http://sis.agr.gc.ca/cansis/taxa/cssc3/chpt17.html)

 

Unsaturated Flow - UZF Package

In order to solve the variably saturated flow problem using the UZF package, it is necessary to specify the relationship of both the relative permeability () and total pressure head () as functions of the water phase saturation ().  Flow through and storage in the unsaturated zone is estimated in the UZF package using a kinematic wave approximation of Richard's equation using the method of characteristics.  The approach includes the assumptions that flow in the unsaturated zone is driven by gravity, that negative potential gradients are ignored, and that flow is vertical through a homogenous column of soil.  The Brooks-Corey function (Brooks-Corey, 1966) is used to estimate the relationship between relative hydraulic conductivity and water saturation:

where:

Krw is the relative hydraulic conductivity [L/T]

ε is the Brooks-Corey exponent [-]

Se is the effective water saturation [-]

θ is the current water saturation, which is a function of pressure head [-]

θr is the residual water saturation [-]

 

The parameters are entered as follows:

EPS - (Brooks-Corey exponent ε) controls the air-entry pressure - the smaller the α value, the larger the capillary fringe above the water table.  Gravels and sands have a large α while clayey soils have a small α.

THTi - (initial water saturation θi ) initial water content in the vadose zone.

THTs - (residual water saturation θr) depends upon several factors of the soil including packing. Residual saturation values are generally low for sands and gravels and high for clayey soils. Residual saturation typically ranges from 0.01 to 0.4.

 

In addition to the property values described above, the UZF also includes transient boundary condition parameters for potential infiltration rates, potential evapotranspiration rates, extinction water depth and water content, as described in the boundary conditions theory section.

 

Unsaturated Flow - MODFLOW-SURFACT

The BCF4 package in MODFLOW SURFACT can solve for variably saturated flow in the vadose zone in three dimensions using the Brooks-Corey equation described above or using the Van Genuchten (1980) equation for relative permeability:

where:

is an empirical parameter [dimensionless]

 

The relationship of the pressure head () [L] to relative saturation is described by van Genuchten (1980) using the following equation:

where:

is an empirical parameter [1/L]

is an empirical parameter [-] and is related to empirical parameter as follows:

is the capillary pressure [L] and is related to pressure head () and elevation (z) by

 

Depending on the saturation model selected, the primary inputs for simulating groundwater flow in the vadose zone using MODFLOW-SURFACT are the van Genuchten parameters (α, β, and θr ) and the Brooks-Corey parameter (ε).  The Brooks-Corey and van Genuchten empirical parameters can be measured in the laboratory for a given soil by plotting a soil-water characteristic curve.

The parameters are entered as follows:

VANAL - (van Genuchten α parameter) controls the air-entry pressure - the smaller the α value, the larger the capillary fringe above the water table.  Gravels and sands have a large α while clayey soils have a small α.

VANBT - (van Genuchten β parameter) controls how rapidly the saturation drops from unity to residual, and is an indicator of the grain size distribution. For larger values of β, the saturation drops rapidly with head (uniform grain size distribution), while for smaller values of β, the drop is more gradual (grain-size, and therefore pore size is more varied). β is greater than unity, but typically can be from 1.3 to 6, depending upon the soil.

VANSR - (residual water saturation θr) depends upon several factors of the soil including packing. Residual saturation values are generally low for sands and gravels and high for clayey soils. Residual saturation typically ranges from 0.01 to 0.4.

Brook - (Brooks-Corey ε)  is the Brooks-Corey empirical exponent.  It must be greater than zero (default value is 0.5).

 

Unsaturated Flow - Vapor Flow MODFLOW-SURFACT

The BCF4 package in SURFACT can also be used to simulate vapor flow through the vadose zone.  Pressure head in the air phase (Pa) and total head in the air phase (hap) are related as follows:

where:

ρw is density of water [L3/T]

ρa is the density of air [L3/T]

g is the gravitational acceleration constant [L/T2]

h is the total head in the air phase [L]

z is the elevation [L]

 

Soil vapor flow options for MODFLOW-SURFACT simulations are set using the advanced settings for the BCF4 package in the Translation Settings workflow step. These are initially based on the project defaults for physical constants in the Model Defaults tab (found under File / Project Defaults / Model Defaults) and include the following:

RHOWPw) is density of water [L3/T]

PHOAPa) is the density of air [L3/T]

VISWw) is the viscosity of water [M/LT]

VISGa) is the viscosity of air [M/LT]

COMPWATw) compressibility of water [LT2/M]

COMPAIRa) compressibility of air [LT2/M]

ATMGP (Patm) is the standard atmospheric pressure [M/LT2]

GRAV (g) is the gravitational acceleration constant [L/T2]

 

Transport Properties

If transport has been enabled at the Define Modeling Objectives step, then you may chose modeling objectives for the transport simulations, which include selections for the retardation and reaction models. If a transport engine has been selected then one or more of the following properties must be defined.

Initial Concentration

In many cases, the historical conditions of the site are unknown, and the contaminant source has been removed or remediated. In such cases, the groundwater contamination is still present and the mass transport simulation must be run forward in time, starting from the existing conditions, to predict the potential downstream impacts. The Initial Concentration properties define the existing conditions of each chemical species at the start of the simulation period.

Initial Concentrations must be defined for each chemical species that you have defined in the Define Modeling Objectives; the default value is 0.  Furthermore, if you have selected a dual or multi-domain model, then initial concentrations must be defined for each species in each immobile domain.

 

Retardation Models

Visual MODFLOW Flex supports a variety of retardation and sorption models for solute transport.

 

Retardation

Retardation is a transport process that slows solute migration.

where:

is the retarded flow velocity of species i [L/T]

is the average linear groundwater flow velocity [L/T]

is the retardation coefficient of species i [-]

 

There are four sorption models for calculating the retardation, as described below:

 

Linear Sorption:

In linear sorption, the sorbed concentration is assumed to be directly proportional to the dissolved concentration in the adjacent pore space:

where:

is the sorbed concentration of species i [M/M] [-]

is the dissolved concentration of species i [M/L3]

is the distribution coefficient of species i [L3/M]

 

The retardation coefficient for linear sorption is calculated as:

where:

is the retardation coefficient of species i [-]

is the dry bulk density of the soil [M/L3]

is the effective porosity of the soil [L3/L3] [-]

is the distribution coefficient of species i [L3/M]

 

When the Linear Sorption model is included in a Transport Simulation, Visual MODFLOW Flex requires the following parameters:

Species Parameters

Kd is a species parameter for the distribution coefficient [1/(M/L3)]

 

Model Parameters

RHOB is the bulk density of the porous media  [M/L3].

 

If a dual domain or multi domain model objective is specified (see below), then one or more of the following parameters may also be added:

Species Parameters

KdIm (or KdIm1) is the species parameter for distribution in the (first) immobile domain  [1/(M/L3)]

KdIm2 is the species parameter for distribution in the (second) immobile domain  [1/(M/L3)]

KdIm3 is the species parameter for distribution in the (third) immobile domain  [1/(M/L3)]

 

Model Parameters

RhoIm (or RHOBim1) is the bulk density of the porous media in the first immobile domain [M/L3]

RHOBim2 is the bulk density of the porous media in the second immobile domain [M/L3]

RHOBim3 is the bulk density of the porous media in the third immobile domain [M/L3]

 

Freundlich Isotherm

The Freundlich Isotherm simulates sorption using an exponential function:

where:

is the sorbed concentration of species i [M/M] [-]

is the dissolved concentration of species i [M/L3]

is the empirical Freundlich constant of species i [(L3/M)a]

is the empirical Freundlich exponent of species i [-]  

 

The retardation coefficient for the Freundlich isotherm is calculated as:

where:

is the retardation coefficient of species i [-]

is the dry bulk density of the soil [M/L3]

is the effective porosity of the soil [L3/L3] [-]

is the empirical Freundlich constant of species i [(L3/M)a]

is the empirical Freundlich exponent of species i [-]  

 

When the Freundlich Isotherm model is included in a Transport Simulation, Visual MODFLOW Flex requires the following parameters:

Species Parameters (one value per cell per species)

Kf is the Freundlich constant [1/(M/L3)a]

a Freundlich exponent [-]

 

Note: MODFLOW-SURFACT also includes support for specifiying a different set of Freundlich isotherm parameters for the immobile domain:
 

KfIm is the Freundlich constant in the immobile domain [1/(M/L3)a]

aim is the Freundlich exponent in the immobile domain [-]

 

Model Parameters (one value per cell)

RHOB is the bulk density of the porous media  [M/L3]

 

Langmuir Isotherm

The Langmuir Isotherm simulates sorption as a non-linear empirical function:

where:

is the sorbed concentration of species i [M/M] [-]

is the dissolved concentration of species i [M/L3]

is the total concentration of sorption sites available for species i [M/M] [-]

is the empirical Langmuir constant of species i [1/(M/L3)]

 

The retardation coefficient for Langmuir isotherm is calculated as:

where:

is the retardation coefficient of species i [-]

is the dry bulk density of the soil [M/L3]

is the effective porosity of the soil [L3/L3] [-]

is the total concentration of sorption sites available for species i [M/M] [-]

is the empirical Langmuir constant of species i [1/(M/L3)]

 

When the Freundlich Isotherm model is included in a Transport Simulation, Visual MODFLOW Flex requires the following parameters:

Species Parameters (one value per cell per species)

Kl is the empirical Langmuir constant [1/(M/L3)]

S is the total concentration of sorption sites available [M/M] [-]

 

Model Parameters (one value per cell)

RHOB is the bulk density of the porous media  [M/L3]

 

First Order Kinetic Sorption

The sorption models described above all assume that solute mass transfer between the dissolved and sorbed phases is instantaneous and in equilibrium. In some cases, this assumption may not be valid and sorption may be represented as a first-order reversible kinetic reaction:

where:

is the sorbed concentration of species i [M/M] [-]

is the dissolved concentration of species i [M/L3]

is the dry bulk density of the soil [M/L3]

is the first-order kinetic mass transfer rate between the dissolved and sorbed phase for species i [1/T]

is the distribution coefficient of species i [L3/M]

 

The retardation rate cannot be solved directly in the case of first order kinetic sorption. Instead, the above equation must incorporated into the governing transport equation and the sorbed/dissolved phase concentrations must be solved simultaneously. As the value of the mass transfer rate () increases, the system will tend towards the linear sorption model and as values of decrease, the system will tend towards negligible sorption.

 

Dual Domain and Multi Domain Mass Transfer

Heterogeneous groundwater systems may be conceptualized as a dual or multi domain system. In such cases, the conceptualized system consists of a primary mobile domain where advection is a dominant transport mechanism and one or more adjacent immobile domains where advection is negligible relative to transport in the mobile domain. Solute mass transfer between the mobile and immobile domain(s) can be approximated as a first-order concentration gradient-driven diffusion-based mass balance for the jth immobile domain is calculated as:

where:

is the dissolved concentration of species i in the mobile domain [M/L3]

is the dissolved concentration of species i in the jth immobile domain [M/L3]

is time [T]

is the first-order mass transfer rate between the mobile domain and the jth immobile domain for species i [1/T]

is the porosity of the jth immobile domain. For dual-domain systems is typically assumed to be the difference between total porosity () and effective porosity () [L3/L3] [-]

 

Mass transfer between immobile domains is assumed to be negligible and is not included in the formulation.

If linear sorption is included in the modeling objectives, the mass balance for the jth immobile domain is calculated as:

where:

is the fraction of the jth immobile domain [-]

is the bulk density of the jth immobile domain; in most cases is assumed to be the same as the bulk density of the mobile fraction [M/L3]

is the linear distribution coefficient for species i in the jth immobile domain [1/(M/L3)]

 

The implementations in MT3D-MS and MODFLO-6 assume that and as such, it is not parameterized.  The MODFLOW-SURFACT formulation of the transport equation defines the fraction of the immobile domain as the difference between the total unit volume of aquifer () and the fraction of the mobile domain , that is and is included as a parameter (PHIF).

Similar equations to the one above can be derived for non-linear sorption isotherms, as described in the MODFLOW-SURFACT Manual (HydroGeoLogic Inc., 1996), which also supports the Freundlich Isotherm as part of a dual domain formulation. However, for the other transport engines supported by Visual MODFLOW Flex that implement a dual domain formulation (e.g. MT3D-MS and MODFLOW-6) are more restrictive; if a dual or multi-domain approach is included in the modeling objectives:

MT3D-MS supports models with linear sorption in the mobile domain only

MODFLOW-6 supports models with Freundlich isotherm or Langmuir isotherm sorption if and only if sorption is not simulated in the immobile domain(s). However, linear sorption may be simulated in some, all, or none of the mobile and immobile domain(s).

 

It should also be noted that MT3D-MS and MODFLOW-SURFACT only support dual domain formulations (i.e. a mobile domain and a single immobile domain. MODFLOW-6 allows support for an unspecified number of immobile domains; however, the Visual MODFLOW Flex interface limits the number of immobile domains to three (as supported by the modeling objective entitled "Quadruple Domain with Complex Sorption").

Decay reactions may also be included in the above formulation. For more details, please refer to the documentation for the relevant transport engine:

MT3D-MS: Zheng, C., and Wang, P., 1999. MT3DMS, A Modular Three-Dimensional Multi-species Transport Model for Simulation of Advection, Dispersion, and Chemical Reactions of Contaminants in Groundwater Systems: Documentation and User’s Guide. U. S. Army Corps of Engineers, U. S. Army Engineer Research and Development Center, Vicksburg, Mississippi, SERDP-99-1.

 

MODFLOW-6: Langevin, C.D., Provost, A.M., Panday, S., and Hughes, J.D., 2022. Documentation for the MODFLOW 6 Groundwater Transport Model: U.S. Geological Survey Techniques and Methods, book 6, chap. A61, 56 p., https://doi.org/10.3133/tm6A61.

 

MODFLOW-SURFACT: HydroGeoLogic, Inc., 1996. MODFLOW-SURFACT ver. 2.2 User’s manual. A three dimensional fully integrated finite difference code for simulating Fluid flow and Transport of contaminant in saturated-unsaturated porous media. Herndon, VA 20170, USA.

 

When a dual domain or multi domain retardation model is included in a Transport Simulation, Visual MODFLOW Flex may require one or more of the parameters following may be required:

Species Parameters

Kd is the species parameter for the distribution coefficient in the mobile domain [1/(M/L3)]

Kl is the species parameter for the Langmuir is constant in the mobile domain [1/(M/L3)]

S  is the species parameter for the total concentration of sorption sites available in the mobile domain [-]

Kf is the species parameter for the Freundlich isotherm constant in the mobile domain [1/(M/L3)a]

a is the species parameter for the Freundlich exponent in the mobile domain [-]

 

K_mass (or K_mass1) is the mass transfer rate between the mobile domain and the (first) immobile domain [1/T]

K_mass2 is the mass transfer rate between the mobile domain and the (second) immobile domain [1/T]

K_mass3) is the mass transfer rate between the mobile domain and the (third) immobile domain [1/T]

 

KdIm (or KdIm1) is the species parameter for the distribution coefficient in the (first) immobile domain  [1/(M/L3)]

KdIm2 is the species parameter for the distribution coefficient in the (second) immobile domain  [1/(M/L3)]

KdIm3 is the species parameter for the distribution coefficient in the (third) immobile domain  [1/(M/L3)]

 

Model Parameters

RhobIm (or RHOBim1) is the bulk density of the porous media in the first immobile domain [M/L3]

RHOBim2 is the bulk density of the porous media in the second immobile domain [M/L3]

RHOBim3 is the bulk density of the porous media in the third immobile domain [M/L3]

 

PhiF is the volumetric fraction of the mobile domain relative to the model domain [-]

PhiIm (or PhiIm1) is the porosity of the (first) immobile domain [-]

PhiIm2 is the porosity of the second immobile domain [-]

PhiIm3 is the porosity of the third immobile domain [-]

 

Notes:

oAdditional species parameters may be required if a non-trivial reaction model is selected (see below).

oSome species and model parameters have different names depending on which soprtion model is specified (e.g. K_mass vs K_mass1). These have the same physical meaning; this first name listed is for the specific dual-domain sorption models, whereas the second name is for the dual-, triple-, and quadruple-domain with complex sorption options supported by MODFLOW-6, which necessitate a numeric suffix to identify/differentiate the indices of the immobile domain(s).

 

Reaction Models

Visual MODFLOW Flex supports a variety of reaction for solute transport:

 

First-Order Irreversible Decay

In first order irreversible decay, the rate of change of solute concentration is proportional to the concentration at a given time.  The mass balance term for first order decay as a function of time, (assuming no other processes), is given by:

where:

is the concentration of species i at time t [M/L3]

is the initial concentration of species i [M/L3]

is the effective porosity of the soil [L3/L3] [-]

is the first order decay rate of species i [1/T]

is elapsed time [T]

 

Assuming the effective porosity does not change over time, can be dropped from both sides of the equation and the general solution for the concentration as a function of time is given by:

 

First-order decay rate coefficients can be derived from half-life values, which are more commonly available. The half-life is the time for the constituent concentration to reach half of its initial concentration (i.e. ). The first order decay rate is related to the half-life as follows:

where:

is the first order decay rate of species i [1/T]

is the half-life of species i [T]

 

The first order decay may be applied to the constituent in supported phases (dissolved or sorbed) and/or domains (mobile or immobile), as supported by the various transport engines. Please refer to the descriptions of Modeling Objectives and Species Parameters for details.

When first order decay (or complex Model Reactions) is included in a Transport Simulation, Visual MODFLOW Flex may require one or more of the parameters following may be required:

K_mobile is the first-order decay rate coefficient for dissolved constituents in both the mobile (and immobile domain, if applicable) [1/T]

 

K_sorbed is the first-order decay rate coefficient for sorbed constituents in both the mobile (and immobile domain, if applicable) [1/T]

 

K1 is the first-order decay rate coefficient for dissolved constituents in the mobile domain [1/T]

 

K1_sorbed is the first-order decay rate coefficient for sorbed constituents in the mobile domain [1/T]

 

K1_im1 is the first-order decay rate coefficient for dissolved consitituents in the first immobile domain [1/T]

 

K1_im1_sorbed is the first-order decay rate coefficient for sorbed consitituents in the first immobile domain [1/T]

 

K1_im2 is the first-order decay rate coefficient for dissolved consitituents in the second immobile domain [1/T]

 

K1_im2_sorbed is the first-order decay rate coefficient for sorbed consitituents in the second immobile domain [1/T]

 

K1_im3 is the first-order decay rate coefficient for dissolved consitituents in the third immobile domain [1/T]

 

K1_im3_sorbed is the first-order decay rate coefficient for sorbed consitituents in the third immobile domain [1/T]

 

Zeroth-Order Irreversible Decay

In zeroth order irreversible decay, the rate of change of solute concentration is constant. Zeroth order decay is most commonly used in situations where bio-geochemical decay (or production) does not depend on solute concentration or when simulating groundwater age, where the groundwater age where the is mathematically analogous to a solute (Goode, 1996).  The mass balance term for zeroth order decay as a function of time, (assuming no other processes), is given by:

where:

is the concentration of species i at time t [M/L3]

is the effective porosity of the soil [L3/L3] [-]

is the zeroth order decay rate of species i [1/T]

is elapsed time [T]

 

A negative value indicates a loss in mass concentration, which is commonly used models of bio-geochemical decay.  A positive value of the zeroth order decay rate indicates a gain in mass concentration, which may be used for simulating a daughter product in some bio-geochemical reactions or for simulating groundwater age as the groundwater ages at a constant rate, (unless the some of the groundwater is accelerated to relativistic speeds, which is generally improbable and also apt to cause significant problems). Assuming the effective porosity does not change over time, can be dropped from both sides of the equation and the general solution for the concentration as a function of time is given by:

   

When zeroth order decay (or complex Model Reactions) is included in a Transport Simulation, Visual MODFLOW Flex may require one or more of the parameters following may be required:

K_mobile is the zeroth-order decay rate coefficient for dissolved constituents in both the mobile (and immobile domain, if applicable) [1/T]

 

K_sorbed is the zeroth-order decay rate coefficient for sorbed constituents in both the mobile (and immobile domain, if applicable) [1/T]

 

K0 is the zeroth-order decay rate coefficient for dissolved constituents in the mobile domain [1/T]

 

K0_sorbed is the zeroth-order decay rate coefficient for sorbed constituents in the mobile domain [1/T]

 

K0_im1 is the zeroth-order decay rate coefficient for dissolved consitituents in the first immobile domain [1/T]

 

K0_im1_sorbed is the zeroth-order decay rate coefficient for sorbed consitituents in the first immobile domain [1/T]

 

K0_im2 is the zeroth-order decay rate coefficient for dissolved consitituents in the second immobile domain [1/T]

 

K0_im2_sorbed is the zeroth-order decay rate coefficient for sorbed consitituents in the second immobile domain [1/T]

 

K0_im3 is the zeroth-order decay rate coefficient for dissolved consitituents in the third immobile domain [1/T]

 

K0_im3_sorbed is the zeroth-order decay rate coefficient for sorbed consitituents in the third immobile domain [M/L]

 

 

RT3D Reaction Modules

RTD3 is an advanced transport engine that includes support for coupled reactions for several well known reaction types:

Instantaneous aerobic degradation of of BTEX;

Six-Species, First-Order, Rate-Limited, BTEX Degradation using Sequential Electron Acceptors;

Rate-Limited Sorption;

Double-Monod Model;

Sequential First-Order Decay; or

Aerobic/Anaerobic PCE/TCE Dechlorination

 

A description of these reactions is beyond the scope of the Visual MODFLOW Flex documentation as RT3D ; however, the reactions are well descbribed in the RT3D documentation by Clement (1997 and 2003).

 

Dispersion

Dispersion is a physical process that tends to ‘disperse’, or spread, the solute mass in the X, Y and Z directions along the advective path of the plume, and acts to reduce the solute concentration (but not overall mass). Dispersion combines two physical processes, mechanical dispersion () which is due to the tortuosity of the flowpaths of the groundwater as the macro-scale as the solute travels through the interconnected pores of the soil and by diffusion () which is caused by random Fickian motion that tends to spread solute mass away from areas of higher concentration. Dispersion is calculated using the equation:

where:

 

is the Dispersion Coefficient [L2/T]

is the Mechanical Dispersion Coefficient [L2/T]

is the diffusion coefficient [L2/T]

is the longitudinal dispersivity [L]

is the longitudinal velocity of flow along the plume migration pathway [L/T]

is the horizontal dispersivity [L]

is the horizontal velocity of flow along the plume migration pathway [L/T]

is the vertical dispersivity [L]

is the vertical velocity of flow along the plume migration pathway [L/T]

is the magnitude of the groundwater velocity [L/T]

 

MT3DMS, RT3D, SEAWAT, and SURFACT

MT3DMS, RT3D, SEAWAT, and SURFACT calculate the Dispersion Coefficient for the mass transport model using the following parameters:

Longitudinal Dispersivity () for each transport grid cell

Ratio of Horizontal to Longitudinal Dispersivity () for each layer

Ratio of Vertical to Longitudinal Dispersivity () for each layer

Molecular Diffusion Coefficient for each layer

 

Longitudinal Dispersion can be defined using the regular set of tools provided for most parameters (i.e. cell-by-cell, by layer/row/column, or using polylines/polygons/data objects). The Longitudinal Dispersion can also be defined on a layer-by-layer basis  by right-clicking the Longitudinal Dispersion object in the model tree and selecting the Dispersion Parameters option.

 

The following dialog will appear:

 

Longitudinal Dispersion can also be defined based on property zones by right-clicking the Longitudinal Dispersion object in the model tree and selecting the Edit Property Zones... option.  The absolute values of the horizontal and vertical dispersion parameters will be adjusted according to the property zone values.

 

For SEAWAT models, you can also specify a diffusion coefficient for each species using the MDCOEFF variable at the define modeling objectives step.  The values will only be used if the Multi-Diffusion option is set to true in the Advanced Settings for the DSP package in the Translation step.

 

MODFLOW-6

MODFLOW-6 includes a more complex conceptualization of the Mechanical Dispersion tensor (please see Langevin, et al, 2022) and calculates the Dispersion Coefficient for the mass transport model using the following parameters on a cell-by-cell basis:

ALH is the longitudinal dispersivity in the horizontal direction. If flow is strictly horizontal, then this is the longitudinal dispersivity that will be used. If flow is not strictly horizontal or strictly vertical, then the longitudinal dispersivity () is a function of both ALH and ALV.

 

ALV is the longitudinal dispersivity in vertical direction. If flow is strictly vertical, then this is the longitudinal dispsersivity value that will be used. If flow is not strictly horizontal or strictly vertical, then the longitudinal dispersivity () is a function of both ALH and ALV.

 

ATH1 is the transverse dispersivity in the horizontal direction. This is the transverse dispersivity value for the second ellipsoid axis. If flow is strictly horizontal and directed in the X direction, then this value controls spreading in the Y direction. If mechanical dispersion is represented (by specifying any dispersivity values) then this array is required.

 

ATH2 is the transverse dispersivity in the horizontal direction. This is the transverse dispersivity value for the third ellipsoid axis. If flow is strictly horizontal and directed in the X direction, then this value controls spreading in the Z direction. If this value is not specified and mechanical dispersion is represented, then this array is set equal to ATH1.

 

ATV is the transverse dispersivity when flow is in vertical direction. If flow is strictly vertical and directed in the Z direction, then this value controls spreading in the X and Y directions. If this value is not specified and mechanical dispersion is represented, then this array is set equal to ATH2.

 

 


Page url:http://www.waterloohydrogeologic.com/help/vmod-flex/index.html?vm_propertyzones_theory.htm