MODFLOW

MODFLOW is the U.S. Geological Survey modular finite-difference flow model, which is a computer code that solves the groundwater flow equation. The program is used by hydrogeologists to simulate the flow of groundwater through aquifers. The source code is free public domain software, written primarily in Fortran, and can compile and run on Microsoft Windows or Unix-like operating systems. The compiled program is a command line tool that reads text file inputs that define the model grid, boundary conditions, and simulation time frame. Extensions beyond the command line program are several actively developed commercial and non-commercial graphical user interfaces.

Groundwater flow equation

The governing partial differential equation for a confined aquifer used in MODFLOW is:

∂

∂

x

[

K

x

x

∂

h

∂

x

]

+

∂

∂

y

[

K

y

y

∂

h

∂

y

]

+

∂

∂

z

[

K

z

z

∂

h

∂

z

]

+

W

=

S

S

∂

h

∂

t

{\displaystyle {\frac {\partial }{\partial x}}\left[K_{xx}{\frac {\partial h}{\partial x}}\right]+{\frac {\partial }{\partial y}}\left[K_{yy}{\frac {\partial h}{\partial y}}\right]+{\frac {\partial }{\partial z}}\left[K_{zz}{\frac {\partial h}{\partial z}}\right]+W=S_{S}{\frac {\partial h}{\partial t}}}

where

K

x

x

{\displaystyle K_{xx}}

,

K

y

y

{\displaystyle K_{yy}}

and

K

z

z

{\displaystyle K_{zz}}

are the values of hydraulic conductivity along the x, y, and z coordinate axes (L/T)

h

{\displaystyle h}

is the potentiometric head (L)

W

{\displaystyle W}

is a volumetric flux per unit volume representing sources and/or sinks of water, where negative values are extractions, and positive values are injections (T−1)

S

S

{\displaystyle S_{S}}

is the specific storage of the porous material (L−1); and

t

{\displaystyle t\,}

is time (T)

Finite difference

The finite difference form of the partial differential in a discretized aquifer domain (represented using rows, columns and layers) is:

C

R

i

,

j

−

1

2

,

k

(

h

i

,

j

−

1

,

k

m

−

h

i

,

j

,

k

m

)

+

C

R

i

,

j

+

1

2

,

k

(

h

i

,

j

+

1

,

k

m

−

h

i

,

j

,

k

m

)

+

C

C

i

−

1

2

,

j

,

k

(

h

i

−

1

,

j

,

k

m

−

h

i

,

j

,

k

m

)

+

C

C

i

+

1

2

,

j

,

k

(

h

i

+

1

,

j

,

k

m

−

h

i

,

j

,

k

m

)

+

C

V

i

,

j

,

k

−

1

2

(

h

i

,

j

,

k

−

1

m

−

h

i

,

j

,

k

m

)

+

C

V

i

,

j

,

k

+

1

2

(

h

i

,

j

,

k

+

1

m

−

h

i

,

j

,

k

m

)

+

P

i

,

j

,

k

h

i

,

j

,

k

m

+

Q

i

,

j

,

k

=

S

S

i

,

j

,

k

(

Δ

r

j

Δ

c

i

Δ

v

k

)

h

i

,

j

,

k

m

−

h

i

,

j

,

k

m

−

1

t

m

−

t

m

−

1

{\displaystyle {\begin{aligned}&{\mathit {CR}}_{i,j-{\tfrac {1}{2}},k}\left(h_{i,j-1,k}^{m}-h_{i,j,k}^{m}\right)+{\mathit {CR}}_{i,j+{\tfrac {1}{2}},k}\left(h_{i,j+1,k}^{m}-h_{i,j,k}^{m}\right)+\\&{\mathit {CC}}_{i-{\tfrac {1}{2}},j,k}\left(h_{i-1,j,k}^{m}-h_{i,j,k}^{m}\right)+{\mathit {CC}}_{i+{\tfrac {1}{2}},j,k}\left(h_{i+1,j,k}^{m}-h_{i,j,k}^{m}\right)+\\&{\mathit {CV}}_{i,j,k-{\tfrac {1}{2}}}\left(h_{i,j,k-1}^{m}-h_{i,j,k}^{m}\right)+{\mathit {CV}}_{i,j,k+{\tfrac {1}{2}}}\left(h_{i,j,k+1}^{m}-h_{i,j,k}^{m}\right)+\\&P_{i,j,k}\,h_{i,j,k}^{m}+Q_{i,j,k}={\mathit {SS}}_{i,j,k}\left(\Delta r_{j}\Delta c_{i}\Delta v_{k}\right){\frac {h_{i,j,k}^{m}-h_{i,j,k}^{m-1}}{t^{m}-t^{m-1}}}\end{aligned}}}

where

h

i

,

j

,

k

m

{\displaystyle h_{i,j,k}^{m}\,}

is the hydraulic head at cell i,j,k at time step m

CV, CR and CC are the hydraulic conductances, or branch conductances between node i,j,k and a neighboring node

P

i

,

j

,

k

{\displaystyle P_{i,j,k}\,}

is the sum of coefficients of head from source and sink terms

Q

i

,

j

,

k

{\displaystyle Q_{i,j,k}\,}

is the sum of constants from source and sink terms, where

Q

i

,

j

,

k

<

0.0

{\displaystyle Q_{i,j,k}<0.0\,}

is flow out of the groundwater system (such as pumping) and

Q

i

,

j

,

k

>

0.0

{\displaystyle Q_{i,j,k}>0.0\,}

is flow in (such as injection)

S

S

i

,

j

,

k

{\displaystyle {\mathit {SS}}_{i,j,k}\,}

is the specific storage

Δ

r

j

Δ

c

i

Δ

v

k

{\displaystyle \Delta r_{j}\Delta c_{i}\Delta v_{k}\,}

are the dimensions of cell i,j,k, which, when multiplied, represent the volume of the cell; and

t

m

{\displaystyle t^{m}\,}

is the time at time step m

This equation is formulated into a system of equations to be solved as:

C

V

i

,

j

,

k

−

1

2

h

i

,

j

,

k

−

1

m

+

C

C

i

−

1

2

,

j

,

k

h

i

−

1

,

j

,

k

m

+

C

R

i

,

j

−

1

2

,

k

h

i

,

j

−

1

,

k

m

+

(

−

C

V

i

,

j

,

k

−

1

2

−

C

C

i

−

1

2

,

j

,

k

−

C

R

i

,

j

−

1

2

,

k

−

C

R

i

,

j

+

1

2

,

k

−

C

C

i

+

1

2

,

j

,

k

−

C

V

i

,

j

,

k

+

1

2

+

H

C

O

F

i

,

j

,

k

)

h

i

,

j

,

k

m

+

C

R

i

,

j

+

1

2

,

k

h

i

,

j

+

1

,

k

m

+

C

C

i

+

1

2

,

j

,

k

h

i

+

1

,

j

,

k

m

+

C

V

i

,

j

,

k

+

1

2

h

i

,

j

,

k

+

1

m

=

R

H

S

i

,

j

,

k

{\displaystyle {\begin{aligned}&{\mathit {CV}}_{i,j,k-{\tfrac {1}{2}}}h_{i,j,k-1}^{m}+{\mathit {CC}}_{i-{\tfrac {1}{2}},j,k}h_{i-1,j,k}^{m}+{\mathit {CR}}_{i,j-{\tfrac {1}{2}},k}h_{i,j-1,k}^{m}\\&+\left(-{\mathit {CV}}_{i,j,k-{\tfrac {1}{2}}}-{\mathit {CC}}_{i-{\tfrac {1}{2}},j,k}-{\mathit {CR}}_{i,j-{\tfrac {1}{2}},k}-{\mathit {CR}}_{i,j+{\tfrac {1}{2}},k}-{\mathit {CC}}_{i+{\tfrac {1}{2}},j,k}-{\mathit {CV}}_{i,j,k+{\tfrac {1}{2}}}+{\mathit {HCOF}}_{i,j,k}\right)h_{i,j,k}^{m}\\&+{\mathit {CR}}_{i,j+{\tfrac {1}{2}},k}h_{i,j+1,k}^{m}+{\mathit {CC}}_{i+{\tfrac {1}{2}},j,k}h_{i+1,j,k}^{m}+{\mathit {CV}}_{i,j,k+{\tfrac {1}{2}}}h_{i,j,k+1}^{m}={\mathit {RHS}}_{i,j,k}\end{aligned}}}

where

H

C

O

F

i

,

j

,

k

=

P

i

,

j

,

k

−

S

S

i

,

j

,

k

Δ

r

j

Δ

c

i

Δ

k

t

m

−

t

m

−

1

R

H

S

i

,

j

,

k

=

−

Q

i

,

j

,

k

−

S

S

i

,

j

,

k

Δ

r

j

Δ

c

i

Δ

v

k

h

i

,

j

,

k

m

−

1

t

m

−

t

m

−

1

{\displaystyle {\begin{aligned}{\mathit {HCOF}}_{i,j,k}&=P_{i,j,k}-{\frac {{\mathit {SS}}_{i,j,k}\Delta r_{j}\Delta c_{i}\Delta _{k}}{t^{m}-t^{m-1}}}\\{\mathit {RHS}}_{i,j,k}&=-Q_{i,j,k}-{\mathit {SS}}_{i,j,k}\Delta r_{j}\Delta c_{i}\Delta v_{k}{\frac {h_{i,j,k}^{m-1}}{t^{m}-t^{m-1}}}\end{aligned}}}

or in matrix form as:

A

h

=

q

{\displaystyle A\mathbf {h} =\mathbf {q} }

where

A is a matrix of the coefficients of head for all active nodes in the grid

h

{\displaystyle \mathbf {h} }

is a vector of head values at the end of time step m for all nodes in the grid; and

q

{\displaystyle \mathbf {q} }

is a vector of the constant terms, RHS, for all nodes of the grid.

Limitations

The water must have a constant density, dynamic viscosity (and consequently temperature) throughout the modelling domain (SEAWAT is a modified version of MODFLOW which is designed for density-dependent groundwater flow and transport)

The principal components of anisotropy of the hydraulic conductivity used in MODFLOW is displayed on the right. This tensor does not allow non-orthogonal anisotropies, as could be expected from flow in fractures. Horizontal anisotropy for an entire layer can be represented by the coefficient "TRPY" (Data Item 3 Page 153).

Versions

"Modular Model"

The USGS throughout the 1970s had developed several hundred models, written in different dialects of FORTRAN. At the time, it was common practice to rewrite a new model to fit the need of a new groundwater scenario.