Note

This page was generated from a jupyter notebook.

Coupling a Landlab groundwater with a Mesa agent-based model

This notebook shows a toy example of how one might couple a simple groundwater model (Landlab’s GroundwaterDupuitPercolator, by Litwin et al. (2020)) with an agent-based model (ABM) written using the Mesa Agent-Based Modeling (ABM) package.

The purpose of this tutorial is to demonstrate the technical aspects of creating an integrated Landlab-Mesa model. The example is deliberately very simple in terms of the processes and interactions represented, and not meant to be a realistic portrayal of water-resources decision making. But the example does show how one might build a more sophisticated and interesting model using these basic ingredients.

(Greg Tucker, November 2021; created from earlier notebook example used in May 2020 workshop)

Running the groundwater model

The following section simply illustrates how to create a groundwater model using the GroundwaterDupuitPercolator component.

Imports:

import matplotlib.pyplot as plt
from tqdm.notebook import trange

from landlab import RasterModelGrid
from landlab.components import GroundwaterDupuitPercolator

Matplotlib is building the font cache; this may take a moment.

Set parameters:

base_depth = 22.0  # depth of aquifer base below ground level, m
initial_water_table_depth = 2.0  # starting depth to water table, m
dx = 100.0  # cell width, m
pumping_rate = 0.001  # pumping rate, m3/s
well_locations = [800, 1200]
K = 0.001  # hydraulic conductivity, (m/s)
n = 0.2  # porosity, (-)
dt = 3600.0  # time-step duration, s
background_recharge = 0.1 / (3600 * 24 * 365.25)  # recharge rate from infiltration, m/s

Create a grid and add fields:

# Raster grid with closed boundaries
# boundaries = {'top': 'closed','bottom': 'closed','right':'closed','left':'closed'}
grid = RasterModelGrid((41, 41), xy_spacing=dx)  # , bc=boundaries)

# Topographic elevation field (meters)
elev = grid.add_zeros("topographic__elevation", at="node")

# Field for the elevation of the top of an impermeable geologic unit that forms
# the base of the aquifer (meters)
grid.at_node["aquifer_base__elevation"] = elev - base_depth

# Field for the elevation of the water table (meters)
grid.at_node["water_table__elevation"] = elev - initial_water_table_depth

# Field for the groundwater recharge rate (meters per second)
recharge = grid.add_full("recharge__rate", background_recharge, at="node")

# pumping rate, in terms of recharge
recharge[well_locations] -= pumping_rate / (dx * dx)

Instantiate the component (note use of an array/field instead of a scalar constant for recharge_rate):

gdp = GroundwaterDupuitPercolator(
    grid,
    hydraulic_conductivity=K,
    porosity=n,
    recharge_rate=recharge,
    regularization_f=0.01,
)

Define a couple of handy functions to run the model for a day or a year:

def run_for_one_day(gdp, dt):
    num_iter = int(3600.0 * 24 / dt)
    for _ in range(num_iter):
        gdp.run_one_step(dt)

def run_for_one_year(gdp, dt):
    num_iter = int(365.25 * 3600.0 * 24 / dt)
    for _ in range(num_iter):
        gdp.run_one_step(dt)

Run for a year and plot the water table:

for _ in trange(365):
    run_for_one_day(gdp, dt)

grid.imshow("water_table__elevation", colorbar_label="Water table elevation (m)")

Aside: calculating a pumping rate in terms of recharge

The pumping rate at a particular grid cell (in volume per time, representing pumping from a well at that location) needs to be given in terms of a recharge rate (depth of water equivalent per time) in a given grid cell. Suppose for example you’re pumping 16 gallons/minute (horrible units of course). That equates to:

16 gal/min x 0.00378541 m3/gal x (1/60) min/sec =

Qp = 16.0 * 0.00378541 / 60.0
print(Qp)

0.0010094426666666667

…equals about 0.001 m\(^3\)/s. That’s \(Q_p\). The corresponding negative recharge in a cell of dimensions \(\Delta x\) by \(\Delta x\) would be

\(R_p = Q_p / \Delta x^2\)

Rp = Qp / (dx * dx)
print(Rp)

1.0094426666666667e-07

A very simple ABM with farmers who drill wells into the aquifer

For the sake of illustration, our ABM will be extremely simple. There are \(N\) farmers, at random locations, who each pump at a rate \(Q_p\) as long as the water table lies above the depth of their well, \(d_w\). Once the water table drops below their well, the well runs dry and they switch from crops to pasture.

Check that Mesa is installed

For the next step, we must verify that Mesa is available. If it is not, use one of the installation commands below to install, then re-start the kernel (Kernel => Restart) and continue.

try:
    from mesa import Model
except ModuleNotFoundError:
    print(
        """
Mesa needs to be installed in order to run this notebook.

Normally Mesa should be pre-installed alongside the Landlab notebook collection.
But it appears that Mesa is not already installed on the system on which you are
running this notebook. You can install Mesa from a command prompt using either:

`conda install -c conda-forge mesa`

or

`pip install mesa`
    """
    )
    raise

Could not import SolaraViz. If you need it, install with 'pip install --pre mesa[viz]'

Defining the ABM

In Mesa, an ABM is created using a class for each Agent and a class for the Model. Here’s the Agent class (a Farmer). Farmers have a grid location and an attribute: whether they are actively pumping their well or not. They also have a well depth: the depth to the bottom of their well. Their action consists of checking whether their well is wet or dry; if wet, they will pump, and if dry, they will not.

from mesa import Agent
from mesa.space import MultiGrid
from mesa.time import RandomActivation

class FarmerAgent(Agent):
    """An agent who pumps from a well if it's not dry."""

    def __init__(self, unique_id, model, well_depth=5.0):
        # super().__init__(unique_id, model)
        super().__init__(model)
        self.pumping = True
        self.well_depth = well_depth

    def step(self):
        x, y = self.pos
        print(f"Farmer {self.unique_id}, ({x}, {y})")
        print(f"    Depth to the water table: {self.model.wt_depth_2d[x, y]}")
        print(f"    Depth to the bottom of the well: {self.well_depth}")
        if self.model.wt_depth_2d[x, y] >= self.well_depth:  # well is dry
            print("    Well is dry.")
            self.pumping = False
        else:
            print("    Well is pumping.")
            self.pumping = True

Next, define the model class. The model will take as a parameter a reference to a 2D array (with the same dimensions as the grid) that contains the depth to water table at each grid location. This allows the Farmer agents to check whether their well has run dry.

import random


class FarmerModel(Model):
    """A model with several agents on a grid."""

    def __init__(self, N, width, height, well_depth, depth_to_water_table, seed=None):
        super().__init__()

        self.num_agents = N
        self.grid = MultiGrid(width, height, True)
        self.depth_to_water_table = depth_to_water_table
        self.schedule = RandomActivation(self)
        self.random = random.Random(seed)

        # Create agents
        for i in range(self.num_agents):
            a = FarmerAgent(i, self, well_depth)
            self.schedule.add(a)
            # Add the agent to a random grid cell (excluding the perimeter)
            x = self.random.randrange(self.grid.width - 2) + 1
            y = self.random.randrange(self.grid.width - 2) + 1
            self.grid.place_agent(a, (x, y))

    def step(self):
        self.wt_depth_2d = self.depth_to_water_table.reshape(
            (self.grid.width, self.grid.height)
        )
        self.schedule.step()

Setting up the Landlab grid, fields, and groundwater simulator

base_depth = 22.0  # depth of aquifer base below ground level, m
initial_water_table_depth = 2.8  # starting depth to water table, m
dx = 100.0  # cell width, m
pumping_rate = 0.004  # pumping rate, m3/s
well_depth = 3  # well depth, m
background_recharge = 0.002 / (365.25 * 24 * 3600)  # recharge rate, m/s
K = 0.001  # hydraulic conductivity, (m/s)
n = 0.2  # porosity, (-)
dt = 3600.0  # time-step duration, s
num_agents = 12  # number of farmer agents
run_duration_yrs = 5  # run duration in years

grid = RasterModelGrid((41, 41), xy_spacing=dx)

elev = grid.add_zeros("topographic__elevation", at="node")

grid.at_node["aquifer_base__elevation"] = elev - base_depth
grid.at_node["water_table__elevation"] = elev - initial_water_table_depth
grid.add_full("water_table__depth_below_ground", initial_water_table_depth, at="node")
recharge = grid.add_full("recharge__rate", background_recharge, at="node")

# pumping rate, in terms of recharge
recharge[well_locations] -= pumping_rate / (dx * dx)

gdp = GroundwaterDupuitPercolator(
    grid,
    hydraulic_conductivity=K,
    porosity=n,
    recharge_rate=recharge,
    regularization_f=0.01,
)

Set up the Farmer model

farmer_model = FarmerModel(
    num_agents,
    grid.number_of_node_columns,
    grid.number_of_node_rows,
    well_depth,
    # depth_to_wt.reshape(grid.shape),
    grid.at_node["water_table__depth_below_ground"].reshape(grid.shape),
    seed=1945,
)

/var/folders/2s/h6hvv9ps03xgz_krkkstvq_r0000gn/T/ipykernel_6669/2212081909.py:13: DeprecationWarning: The time module and all its Schedulers are deprecated and will be removed in Mesa 3.1. They can be replaced with AgentSet functionality. See the migration guide for details. https://mesa.readthedocs.io/latest/migration_guide.html#time-and-schedulers
  self.schedule = RandomActivation(self)

Check the spatial distribution of wells:

import numpy as np


def get_well_count(model):
    well_count = np.zeros(grid.shape, dtype=int)
    pumping_well_count = np.zeros(grid.shape, dtype=int)
    for cell in model.grid.coord_iter():
        cell_content, (x, y) = cell
        well_count[x][y] = len(cell_content)
        for agent in cell_content:
            if agent.pumping:
                pumping_well_count[x][y] += 1
    return well_count, pumping_well_count


well_count, p_well_count = get_well_count(farmer_model)
grid.imshow(well_count)

Set the initial recharge field

recharge[:] = -(pumping_rate / (dx * dx)) * p_well_count.flatten()
grid.imshow(-recharge * 3600 * 24, colorbar_label="Pumping rate (m/day)")

Run the model

depth_to_wt = grid.at_node["water_table__depth_below_ground"]

for i in trange(run_duration_yrs):
    # Run the groundwater simulator for one year
    run_for_one_year(gdp, dt)

    # Update the depth to water table
    depth_to_wt[:] = (
        grid.at_node["topographic__elevation"] - grid.at_node["water_table__elevation"]
    )

    # Run the farmer model
    farmer_model.step()

    # Count the number of pumping wells
    well_count, pumping_well_count = get_well_count(farmer_model)
    total_pumping_wells = np.sum(pumping_well_count)
    print(f"In year {i + 1} there are {total_pumping_wells} pumping wells")
    print(f" and the greatest depth to water table is {np.amax(depth_to_wt)} meters.")

    # Update the recharge field according to current pumping rate
    recharge[:] = (
        background_recharge - (pumping_rate / (dx * dx)) * pumping_well_count.flatten()
    )
    print(f"Total recharge: {np.sum(recharge)}")
    print("")

    plt.figure()
    grid.imshow("water_table__elevation")
    plt.show()

Farmer 10, (36, 6)
    Depth to the water table: 2.9532171973868735
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 3, (19, 24)
    Depth to the water table: 3.125042438222149
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 1, (27, 20)
    Depth to the water table: 3.108003201164557
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 9, (17, 39)
    Depth to the water table: 2.887000015621961
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 11, (25, 12)
    Depth to the water table: 3.116970742578278
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 5, (14, 16)
    Depth to the water table: 3.1509947334725688
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 12, (36, 16)
    Depth to the water table: 2.9843739795064295
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 4, (25, 6)
    Depth to the water table: 3.037622566192745
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 6, (31, 37)
    Depth to the water table: 2.9320254253580877
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 7, (14, 9)
    Depth to the water table: 3.0812406798846723
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 2, (16, 20)
    Depth to the water table: 3.1586006697665425
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 8, (8, 13)
    Depth to the water table: 3.056278550336227
    Depth to the bottom of the well: 3
    Well is dry.
In year 1 there are 4 pumping wells
 and the greatest depth to water table is 3.1586006697665425 meters.
Total recharge: -1.4934646487692335e-06

Farmer 9, (17, 39)
    Depth to the water table: 2.8769211378838335
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 11, (25, 12)
    Depth to the water table: 2.824266525495407
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 2, (16, 20)
    Depth to the water table: 2.815141441130109
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 4, (25, 6)
    Depth to the water table: 2.8153339866056015
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 6, (31, 37)
    Depth to the water table: 2.9138457953766945
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 10, (36, 6)
    Depth to the water table: 2.9231182726368097
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 8, (8, 13)
    Depth to the water table: 2.805415000400476
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 7, (14, 9)
    Depth to the water table: 2.8079596555990136
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 3, (19, 24)
    Depth to the water table: 2.818472450469848
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 1, (27, 20)
    Depth to the water table: 2.8282651014507287
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 5, (14, 16)
    Depth to the water table: 2.811962766912295
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 12, (36, 16)
    Depth to the water table: 2.929668376739542
    Depth to the bottom of the well: 3
    Well is pumping.
In year 2 there are 12 pumping wells
 and the greatest depth to water table is 2.929668376739542 meters.
Total recharge: -4.693464648769233e-06

Farmer 11, (25, 12)
    Depth to the water table: 3.1143755941466367
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 7, (14, 9)
    Depth to the water table: 3.0790560260349658
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 1, (27, 20)
    Depth to the water table: 3.1051900454472623
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 4, (25, 6)
    Depth to the water table: 3.0359216454975524
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 5, (14, 16)
    Depth to the water table: 3.1481775464319455
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 9, (17, 39)
    Depth to the water table: 2.886638832506346
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 12, (36, 16)
    Depth to the water table: 2.9831172831320423
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 8, (8, 13)
    Depth to the water table: 3.054297728905631
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 10, (36, 6)
    Depth to the water table: 2.952407294852545
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 3, (19, 24)
    Depth to the water table: 3.1220281602523556
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 6, (31, 37)
    Depth to the water table: 2.9312215240038277
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 2, (16, 20)
    Depth to the water table: 3.155575445940869
    Depth to the bottom of the well: 3
    Well is dry.
In year 3 there are 4 pumping wells
 and the greatest depth to water table is 3.155575445940869 meters.
Total recharge: -1.4934646487692335e-06

Farmer 5, (14, 16)
    Depth to the water table: 2.811890486026062
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 10, (36, 6)
    Depth to the water table: 2.9231062348845462
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 3, (19, 24)
    Depth to the water table: 2.818391536671431
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 12, (36, 16)
    Depth to the water table: 2.929643146222773
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 6, (31, 37)
    Depth to the water table: 2.9138327971284497
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 2, (16, 20)
    Depth to the water table: 2.8150602977777908
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 8, (8, 13)
    Depth to the water table: 2.8053722758515605
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 9, (17, 39)
    Depth to the water table: 2.876914619030277
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 7, (14, 9)
    Depth to the water table: 2.807910311977299
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 11, (25, 12)
    Depth to the water table: 2.8242027101142426
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 4, (25, 6)
    Depth to the water table: 2.81529819732798
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 1, (27, 20)
    Depth to the water table: 2.8281922819535943
    Depth to the bottom of the well: 3
    Well is pumping.
In year 4 there are 12 pumping wells
 and the greatest depth to water table is 2.929643146222773 meters.
Total recharge: -4.693464648769233e-06

Farmer 2, (16, 20)
    Depth to the water table: 3.1555734391319206
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 8, (8, 13)
    Depth to the water table: 3.0542966780127614
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 1, (27, 20)
    Depth to the water table: 3.1051882524011702
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 6, (31, 37)
    Depth to the water table: 2.93122120939519
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 7, (14, 9)
    Depth to the water table: 3.0790548104511437
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 12, (36, 16)
    Depth to the water table: 2.9831166700784735
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 4, (25, 6)
    Depth to the water table: 3.0359207669681787
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 11, (25, 12)
    Depth to the water table: 3.114374020719641
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 10, (36, 6)
    Depth to the water table: 2.9524070028017135
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 5, (14, 16)
    Depth to the water table: 3.14817575877451
    Depth to the bottom of the well: 3
    Well is dry.
Farmer 9, (17, 39)
    Depth to the water table: 2.8866386746886192
    Depth to the bottom of the well: 3
    Well is pumping.
Farmer 3, (19, 24)
    Depth to the water table: 3.12202616458503
    Depth to the bottom of the well: 3
    Well is dry.
In year 5 there are 4 pumping wells
 and the greatest depth to water table is 3.1555734391319206 meters.
Total recharge: -1.4934646487692335e-06

# Display the area of water table that lies below the well depth
too_deep = (
    grid.at_node["topographic__elevation"] - grid.at_node["water_table__elevation"]
    > well_depth
)

grid.imshow(too_deep)

---

Generated by nbsphinx from a Jupyter notebook.