MODFLOWObservedSimulatedHeadComparisonPlotwithModelMuseandFlopy3.jpg

Groundwater modeling with MODFLOW and other codes are defined as inverse modeling where the aquifer parameters can be calculated from the comparison of model results with observed data. This comparison process is time consuming, employs acceptance criteria and trend analysis of the boundary condition influence.

There are tools for the comparison of observed and simulated heads in Model Muse and custom charts can be done with few lines in Python. This tutorial cover the whole procedure to create a simulated / observed plot in Python from the results of a MODFLOW model run on Model Muse. The study case is over a regional model with more than 100 piezometers. The tutorial creates a graph with a colorbar and exports it as a JPG file.

Tutorial

Python code

This is the code in Python for the tutorial:

# Import the required libraries
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# Setup the model folder and model name
modPath = '../Model/'
modName = 'Modelo1'

# Import observation results as a Pandas dataframe
hobOut = pd.read_csv(modPath + modName + '.hob_out',delimiter='\s+')
hobOut.head()

	SIMULATED EQUIVALENT	OBSERVED VALUE	OBSERVATION NAME
0	3741.871338	3758.0	SD_04_1
1	3744.384277	3763.0	SD_04_2
2	3748.123779	3767.0	SD_04_3
3	3765.178711	3776.0	SD_04_4
4	3761.951904	3784.0	SD_05_5

# Column names are renamed. There is one more column on the dataframe from the difference in between observed - simulated
hobOut.rename(columns={'SIMULATED EQUIVALENT':'Simulated Head',
                      'OBSERVED VALUE':'Observed Head',
                      'OBSERVATION NAME':'Piezometer'},inplace=True)
hobOut.set_index('Piezometer', inplace=True)
hobOut['Residual'] = hobOut['Simulated Head'] - hobOut['Observed Head']
hobOut.head()

	Simulated Head	Observed Head	Residual
Piezometer
SD_04_1	3741.871338	3758.0	-16.128662
SD_04_2	3744.384277	3763.0	-18.615723
SD_04_3	3748.123779	3767.0	-18.876221
SD_04_4	3765.178711	3776.0	-10.821289
SD_05_5	3761.951904	3784.0	-22.048096

#Piezometers that went dry are filtered from the plot, because there is no mean to compare a simulated dry piezometer with
#a observed head
hobOutFil = hobOut[hobOut['Simulated Head']>-10000]
hobOutFil.head()

	Simulated Head	Observed Head	Residual
Piezometer
SD_04_1	3741.871338	3758.0	-16.128662
SD_04_2	3744.384277	3763.0	-18.615723
SD_04_3	3748.123779	3767.0	-18.876221
SD_04_4	3765.178711	3776.0	-10.821289
SD_05_5	3761.951904	3784.0	-22.048096

#Create a basic scatter plot with Matplotlib
#A identity line was created based on the minumum and maximum observed value
#Points markers are colored by the residual and a residual colorbar is added to the figure
fig = plt.figure(figsize=(10,8))
x = np.linspace(hobOutFil['Observed Head'].min(), hobOutFil['Observed Head'].max(), 100)
plt.plot(x, x, linestyle='dashed')
plt.scatter(hobOutFil['Observed Head'],hobOutFil['Simulated Head'], marker='o', c=hobOutFil['Residual'])

cbar = plt.colorbar()
cbar.set_label('Residual (m)', fontsize=14)

plt.grid()
plt.xlabel('Observed Head (m)', fontsize=14)
plt.ylabel('Simulated Head (m)', fontsize=14)
fig.tight_layout()
fig.savefig('../Figures/ObservedSimulatedHeadComparisonPlot.png')