Modflow 6 has changed in many ways the traditional concepts of groundwater modeling; now a group of models can run under a "simulation" and these models can interchange flow and transport. Under this structure Flopy has implemented a tool to split a groundwater model into several models where their results can be reconstructed and compared to the original entire models. Capabilities and features like this make Flopy and Modflow 6 a great tool for developing advanced groundwater models that simulate complicated tasks or requirements to the groundwater flow regime.

Tutorial

Input data

You can download the input data from this link

Code

import os, sys, copy
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import flopy

from flopy.mf6.utils import Mf6Splitter
from flopy.plot import styles
from flopy.utils.geometry import LineString, Polygon

#load workspace and set new sim path
workspace = os.path.abspath('../Mf6Model/')
sim = flopy.mf6.MFSimulation.load('mfsim', sim_ws=workspace)
sim.set_sim_path('../Mf6Model/entire')
sim.exe_name = 'c://WRDAPP/mf6.exe'

loading simulation...
  loading simulation name file...
  loading tdis package...
  loading model gwf6...
    loading package dis...
    loading package ic...
WARNING: Block "options" is not a valid block name for file type ic.
    loading package npf...
    loading package obs...
    loading package sto...
    loading package oc...
    loading package ghb...
    loading package wel...
    loading package riv...
    loading package rch...
  loading solution package modflow...

#get model
gwf = sim.get_model()

#write simulation and run
sim.write_simulation()
success, buff = sim.run_simulation(silent=True)
assert success

writing simulation...
  writing simulation name file...
  writing simulation tdis package...
  writing solution package modflow...
  writing model modflow...
    writing model name file...
    writing package dis...
    writing package ic...
    writing package npf...
    writing package obs-1...
    writing package sto...
    writing package oc...
    writing package ghb-1...
    writing package wel-1...
    writing package riv-1...
    writing package rch-1...

#Get model extent
gwf.modelgrid.extent

(198748.1387899443, 215789.41221654008, 8792833.583466, 8809874.856892597)

#Set corner coordinates
Xmin = gwf.modelgrid.extent[0]
Xmax = gwf.modelgrid.extent[1]
Ymin = gwf.modelgrid.extent[2]
Ymax = gwf.modelgrid.extent[3]

# Plot the model results

# +
water_table = flopy.utils.postprocessing.get_water_table(
    gwf.output.head().get_data()
)
heads = gwf.output.head().get_alldata()[-1]
hmin, hmax = water_table.min(), water_table.max()
contours = np.arange(0, 10, 1)
hmin, hmax

(5.077448514558371e-06, 3.353759226075293)

# +
with styles.USGSMap():
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot()
    ax.set_xlim(Xmin, Xmax)
    ax.set_ylim(Ymin, Ymax)
    ax.set_aspect("equal")
    pmv = flopy.plot.PlotMapView(modelgrid=gwf.modelgrid, ax=ax)
    h = pmv.plot_array(heads, vmin=hmin, vmax=hmax)
    c = pmv.contour_array(
        water_table,
        levels=contours,
        colors="cyan",
        linewidths=1.5,
        linestyles=":",
    )
    plt.clabel(c, fontsize=8)
    pmv.plot_inactive(alpha=0.1)
    plt.colorbar(h, ax=ax, shrink=0.5)
    plt.show()
# -

# ### Split the watershed model
#
# Build a splitting array and split this model into many models for parallel modflow runs

nrow = gwf.modelgrid.nrow
ncol = gwf.modelgrid.ncol 

nrow_blocks, ncol_blocks = 2, 4
row_inc, col_inc = int(nrow / nrow_blocks), int(ncol / ncol_blocks)
row_inc, col_inc

# +
icnt = 0
row_blocks = [icnt]
for i in range(nrow_blocks):
    icnt += row_inc
    row_blocks.append(icnt)
if row_blocks[-1] < nrow:
    row_blocks[-1] = nrow
row_blocks
# -

# +
icnt = 0
col_blocks = [icnt]
for i in range(ncol_blocks):
    icnt += col_inc
    col_blocks.append(icnt)
if col_blocks[-1] < ncol:
    col_blocks[-1] = ncol
col_blocks
# -

# +
mask = np.zeros((nrow, ncol), dtype=int)
# -

# +
# create masking array
ival = 1
model_row_col_offset = {}
for idx in range(len(row_blocks) - 1):
    for jdx in range(len(col_blocks) - 1):
        mask[
            row_blocks[idx] : row_blocks[idx + 1],
            col_blocks[jdx] : col_blocks[jdx + 1],
        ] = ival
        model_row_col_offset[ival - 1] = (row_blocks[idx], col_blocks[jdx])
        # increment model number
        ival += 1
# -

# +
fig = plt.figure()
plt.imshow(mask, aspect=.2)
plt.show()
# -

# ### Now split the model into many models using `Mf6Splitter()`

mfsplit = Mf6Splitter(sim)
new_sim = mfsplit.split_model(mask)

new_sim.set_sim_path('../Mf6Model/splitted')
new_sim.exe_name='c://WRDAPP/mf6.exe'
new_sim.write_simulation()
success, buff = new_sim.run_simulation(silent=True)
assert success
# -

writing simulation...
  writing simulation name file...

# ### Reassemble the heads to the original model shape for plotting
#
# Create a dictionary of model number : heads and use the `reconstruct_array()` method to get a numpy array that is the original shape of the unsplit model.

model_names = list(new_sim.model_names)
head_dict = {}
for modelname in model_names:
    mnum = int(modelname.split("_")[-1])
    head = new_sim.get_model(modelname).output.head().get_alldata()[-1]
    head_dict[mnum] = head

ra_heads = mfsplit.reconstruct_array(head_dict)
ra_watertable = flopy.utils.postprocessing.get_water_table(ra_heads)

# +
with styles.USGSMap():
    fig, axs = plt.subplots(nrows=3, figsize=(12, 24))
    diff = ra_heads - heads
    hv = [ra_heads, heads, diff]
    titles = ["Multiple models", "Single model", "Multiple - single"]
    for idx, ax in enumerate(axs):
        ax.set_aspect("equal")
        ax.set_title(titles[idx])

        if idx < 2:
            levels = contours
            vmin = hmin
            vmax = hmax
        else:
            levels = None
            vmin = None
            vmax = None

        pmv = flopy.plot.PlotMapView(modelgrid=gwf.modelgrid, ax=ax, layer=0)
        h = pmv.plot_array(hv[idx], vmin=vmin, vmax=vmax)
        if levels is not None:
            c = pmv.contour_array(
                hv[idx],
                levels=levels,
                colors="white",
                linewidths=0.75,
                linestyles=":",
            )
            plt.clabel(c, fontsize=8)
        pmv.plot_inactive()
        plt.colorbar(h, ax=ax, shrink=0.5)
    plt.show()

# -