PCRaster Modflow Python

This chapter describes the usage of the PCRasterModflow extension within Python scrips. Using the extension requires the specification of the geographic extents and the creation of an extension object:

from pcraster import *

setclone("clone.map")
# Construct Modflow extension object
mf = initialise(clone())

This will initialise the data structures used in the extension. All operations will operate on the object mf.

The next step is the grid specification using the operations described in the DIS package. This must be done before any other package is defined. Afterwards packages can be defined in arbitrary order. For a Modflow simulation at least the DIS, BAS and BCF package must be specified.

The DIS, BAS, BCF and a solver package must be set in the initial section of a script. Stress packages (RIV, DRN, RCH and WEL) can be activated and modified in the dynamic section.

Note

One timestep in PCRaster represents one stress period in Modflow.

In the following non-spatial arguments are written capitalised.

Setting Modflow’s input data

These sections describe the provided operations for the input data of Modflow:

• Discretization (DIS) package
• Basic (BAS) package
• Block-Centered Flow (BCF) package
• River (RIV) package
• Recharge (RCH) package
• Well (WEL) package
• Drain (DRN) package
• General Head Boundary (GHB) package
• Solver packages (PCG, SOR, SIP, DSP)

Running Modflow

After specifying the required packages, Modflow can be called with

mf.run()

After a successful Modflow run the results of a stress period can be reported with the commands described in this section. If the run fails the end of the Modflow list file is displayed, the execution of the script ends.

In Python you can also specify a directory where the Modflow run will be executed, for example in 1:

mf.run("1")

Retrieving Modflow’s output

The head and boundary values are retrieved automatically and must not be set again for the next stress period. Applying the operations to layer specified as quasi-3D confining beds will result in an error.

head = mf.getHeads(LAYER)

head is a scalar PCRaster map containing the resulting head values of the layer layer. Cells converted to dry obtain the data type missing value.

river = mf.getRiverLeakage(LAYER)

river is a scalar PCRaster map containing the resulting river cell-by-cell flow values (in [\(L^3T^-1\)]) of the layer LAYER.

rch = mf.getRecharge(LAYER)

rch is a scalar PCRaster map containing the resulting recharge cell-by-cell flow values (in [\(L^3T^-1]\)) of the layer LAYER.

drn = mf.getDrain(LAYER)

drn is a scalar PCRaster map containing the resulting drain cell-by-cell flow values (in [\(L^3T^-1]\)) of the layer LAYER.

st = mf.getStorage(LAYER)

st is a scalar PCRaster map containing the resulting cell-by-cell storage terms of the layer LAYER.

ch = mf.getConstantHead(LAYER)

ch is a scalar PCRaster map containing the resulting cell-by-cell constant head flow terms of the layer LAYER.

rf = mf.getRightFace(LAYER)

rf is a scalar PCRaster map containing the resulting internal cell-by-cell flows (right) of the layer LAYER.

ff = mf.getFrontFace(LAYER)

ff is a scalar PCRaster map containing the resulting internal cell-by-cell flows (front) of the layer LAYER.

lf = mf.getLowerFace(LAYER)

lf is a scalar PCRaster map containing the resulting internal cell-by-cell flows (lower) of the layer LAYER.

gh = mf.getGeneralHeadLeakage(LAYER)

gh is a scalar PCRaster map containing the resulting internal cell-by-cell flows terms of the layer LAYER.