Basic Example of a MODFLOW Model Review, Simulation and Output Representation with Flopy - Tutorial

Flopy is a package of tools written in Python for MODFLOW groundwater flow model construction, simulation and output analysis. Flopy is build on top of well know and powerful Python packages as Numpy and works with Matplotlib and Pandas that allows to do a great amount of analysis with few lines of code. Several new capabilities in the water balance analysis can be done with Flopy bringing a better control to the modeler in terms of a more available and user friendly information of the inputs, outputs and discrepancies of the model. This tutorial shows the complete procedure to read, simulate and output analysis of a MODFLOW NWT model of a tunnel development with time. The tutorial include a discussion and review of the different tools available in Flopy and the interaction with QGIS.

You can download Flopy and access to more tutorial and documentation from its Github page:

https://github.com/modflowpy/flopy

 

Tutorial

Input files

You can download the input files from this link.

Code

This is the complete unformatted Python code used in the tutorial:

# coding: utf-8

# # Import and display packages

# In[ ]:

get_ipython().magic('matplotlib inline')
from IPython.display import clear_output, display
from __future__ import print_function
import time
import os
import sys
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import flopy

print(sys.version)
print('numpy version: {}'.format(np.__version__))
print('matplotlib version: {}'.format(mpl.__version__))
print('pandas version: {}'.format(pd.__version__))
print('flopy version: {}'.format(flopy.__version__))


# # Executables and model files

# In[ ]:

#route to the executable file
#in Windows
#modflownwt = '/MODFLOW-NWT_64.exe'
#in Ubuntu/MacOS
modflownwt = os.path.join(os.getcwd(), 'mfnwt')

# first lets load an existing model
model_ws = 'Model'
ml = flopy.modflow.Modflow.load("ModeloDrenajeInfraestructura.nam", model_ws=model_ws, verbose=False,
                               check=False, exe_name=modflownwt)


# In[ ]:

#display model spatial reference
ml.dis.sr


# In[ ]:

#display model executable name
ml.exe_name


# In[ ]:

# display model heading
ml.heading


# In[ ]:

# check modflow options in Flopy
ml.version_types


# # Display model discretization and hydraulic parameters

# In[ ]:

#layer thickness and interface elevation
ml.dis.plot();


# In[ ]:

#horizontal anisotropy, hydraulic conductivy, specific yield, specific storage, etc
ml.upw.plot(mflay=1)


# # Display tunel as drain on different stress periods

# In[ ]:

#ml.drn.plot(mflay=2, kper=1)
ml.drn.plot(mflay=2, kper=1)


# # Export model features to shapefile

# In[ ]:

#Remember to have pyshp installed
#pip install pyshp
ml.dis.export(os.path.join('Shp',"Tunel_dis.shp"))


# In[ ]:

ml.drn.export(os.path.join('Shp',"Tunel_drn.shp"))


# In[ ]:

ml.rch.export(os.path.join('Shp',"Tunel_rch.shp"))


# # Change model working forder and run model

# In[ ]:

#ml.external_path = "ref"
ml.model_ws = "Run"


# In[ ]:

ml.write_input()

ml.run_model()


# # Read water budget per time step as pandas dataframe

# In[ ]:

mfl = flopy.utils.MfListBudget(os.path.join(os.getcwd(),"Run/ModeloDrenajeInfraestructura.lst"))
df_flux, df_vol = mfl.get_dataframes(start_datetime="10-21-2015")
df_flux.head()


# In[ ]:

#plot water inflows with time
df_flux['DRAINS_OUT'].plot()


# In[ ]:

#graph of the water balance at the beginning of tunnel boring

df_flux.rename(columns={'DRAINS':'DRAINS_IN'},inplace=True)

groups = df_flux.groupby(lambda x:x.split('_')[-1],axis=1).groups
df_flux_in = df_flux.loc[:,groups["IN"]]
df_flux_in.columns = df_flux_in.columns.map(lambda x:x.split('_')[0])

df_flux_out = df_flux.loc[:,groups["OUT"]]
df_flux_out.columns = df_flux_out.columns.map(lambda x:x.split('_')[0])

df_flux_delta = df_flux_in - df_flux_out
df_flux_delta
df_flux_delta.iloc[1,:].plot(kind="bar",figsize=(10,10),grid=True, color='darkcyan');


# # Spatial water head representation 

# In[ ]:

# display model times
h = flopy.utils.formattedfile.FormattedHeadFile(os.path.join(os.getcwd(),"Run/ModeloDrenajeInfraestructura.fhd"),model=ml)
h.times


# In[ ]:

#plot specific time for specified layer
h.plot(totim=2.73312e+08, contour=True, mflay=2, grid=False, colorbar=True, 
       figsize=(10,10), masked_values=[-2.00000E+20,-1.00000E+20]);


# In[ ]:

#export heads as shapefile
h.to_shapefile(os.path.join('Shp',"Tunel_hds.shp"))


# In[ ]:

#create animation
f = plt.figure(figsize=(10,10))
ax = plt.subplot(1,1,1,aspect="equal")
for i,t in enumerate(h.times[:-6]):
    ax.set_title("totim:{0}".format(t))
    
    #ml.drn.plot()
    h.plot(totim=t,mflay=3,contour=True,
           grid=True,figsize=(10,10),
           masked_values=[-2.00000E+20,-1.00000E+20],
           axes=[ax])
    ax.set_ylim([6000,11000])
    ax.set_xlim([11000,16000])
    time.sleep(1)
    clear_output(True)
    display(f)
    ax.cla()


# In[ ]: