Setup 3D steady-state finite difference model from scratch.
Specifying model dimensions
The 3D steady state FDM will be based on a regular grid consisting of rows an columns and layers...
Column widths and row heights are constant.
Layer thickness can vary on cell to cell basis.
Specify rectangular grid with the following dimensions:
$x$ and $y$: -1000 to 1000 by 25 meters intervals
$z$:
## Load libraries
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve # import with from to use its short name
import matplotlib.pyplot as plt # combines namespace of numpy and pyplot
## Setup of the model by specifying its dimensions
# specify a rectangular grid
x = np.arange(-1000., 1000., 25.)
y = np.arange(-1000., 1000., 25.) # backward, i.e. first row grid line has highest y
z = np.arange(-100., 0., 20.) # backward, i.e. from top to bottom
From these coordinates we obtain the number of cells along each axis and the cell sizes.
# As well as the number of cells along the three axes
Nx = len(x)-1
Ny = len(y)-1
Nz = len(z)-1
sz = (Nz,Ny,Nx) # the shape of the model
Nod = np.prod(sz) # total number of cells in the model
# From this we have the width of columns, rows and layers
dx = np.diff(x).reshape(1, 1, Nx)
dy = np.diff(y).reshape(1, Ny,1)
dz = np.abs(np.diff(z)).reshape(Nz, 1,1)
IBOUND array - Active, non-active, and prescribed head
Specify cells with prescribed head (fixed head) and inactive cells (zero-conductance). The IBOUND arrray (MODFLOW terminology) is an array with the same dimensions as the cell grid. Inactive cells are cells are not a part of the compuitation (represents impermeable boundaries, such as bed rock, etc.)..
IBOUND > 0, heads will be computed
IBOUND == 0, means cells are inactive
IBOUND < 0, prescribed heads
In this particular example we specify that the vertical $z$-$x$ plane at the last row of the model will have prescribed heads equal to zero.
IBOUND = np.ones(sz)
IBOUND[:,-1,:] = -1 # Last row of model heads are prescribed
IBOUND[:, 40:45, 20:70] = 0 # These cells are inactive
To figure out which cells are active, inactive or have a prescribed head, we can easily peform queries using simple algebreaic equations:
active = (IBOUND>0).reshape(Nod) # active is now a vector of booleans of length Nod
inact = (IBOUND==0).reshape(Nod) # dito for inact
fxhd = (IBOUND<0).reshape(Nod) # dito for fxhd
Cell conductancies: defining the ease of flow between adjacent cells
Flow resistance is defined for each cell in the model in 3 grid directions ($x$, $y$, and $z$). The flow resistance is represented using the hydraulic conductivity in 3D arrays $kx$, $ky$, and $kz$.
k = 10.0 # m/d uniform conductivity
kx = k * np.ones(sz) # [L/T] 3D kx array
ky = k * np.ones(sz) # [L/T] 3D ky array with same values as kx
kz = k * np.ones(sz) # [L/T] 3D kz array with same values as kx
Inactive cells are respresented by setting their resistance to infinite.
Compute the conductance between each pair of adjacent cells across their connecting faces. The conductance is just the reverse of the resistance of the two connected cells in $x$, $y$, and $z$. The resistance in each direction is the sum of the two half cells because the cells are in series in terms of the flow directions in $x$, $y$, and $z$.
# Conductances between adjacent cells
Cx = 1 / (Rx[:, :, :-1] + Rx[:, :,1:]) # [L2/T] in x-direction
Cy = 1 / (Ry[:, :-1,:] + Ry[:, 1:,:]) # idem in y-direction
Cz = 1 / (Rz[:-1,:,:] + Rz[1:,:,:]) # idem in z-direction
Setting up the system matrix - set of water balance equations
To be able to indentify adjacent cells we generate cell numbers in an array that has the size of the model grid:
NOD = np.arange(Nod).reshape(sz) # this is a full 3D array of node numbers (cell numbers)
Now we can easily find the neighbour cells by identifying the adjacent cells using their cells numbers. We do this by creating arrays where the cell ID's are shifted one place in each direction.
IE = NOD[:, :, 1:] # Numbers of the eastern neighbors of each cell
IW = NOD[:, :, :-1] # Same western neighbors
IN = NOD[:, :-1,:] # Same northern neighbors
IS = NOD[:, 1:,:] # Southern neighbors
IT = NOD[:-1,:,:] # Top neighbors
IB = NOD[1:,:,:] # Bottom neighbors
Notice that the shape of the $IE$ and $IW$ is the same as that of $C_x$, the size of $IN$ and $IS$ is the same as that of $C_y$ and the size of $IT$ and $IB$ is the same as that of $C_x$.