How to insert and read Observation Points (OBS6) in Modflow 6 with Model Muse and Flopy – Tutorial

Hydrate

Modflow 6 has a new approach in setting up observation points and it’s essentially different to the previous versions. The OBS6 package works not only with heads and drawdowns but also with flows, so it’s also possible to calibrate the model against baseflow or any other recorded flow from a boundary condition directly. We have created an applied case of the implementation of piezometers in a hillslope groundwater flow model in Modflow 6 and Model Muse. The tutorial covers all the steps related to the implementation of the observed points in Model Muse as well as the comparison among simulated and observed heads through scripts in Flopy.

Tutorial

Code

This is the code for simulated vs observed heads representation:

#import the require packages
import flopy
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
#define model name and model path
modelWs = '../Model/'
modelName = 'HillslopeModel_v1'
namFile = modelName+'.nam'
#open the observation file out heads
obsOut = flopy.utils.observationfile.Mf6Obs(modelWs+modelName+'.ob_gw_out_head.csv',isBinary=False)
#create a dataframe for the observation points
totalObsDf = obsOut.get_dataframe().transpose()
totalObsDf = totalObsDf.rename(columns={totalObsDf.keys()[0]:'SimHead'})
totalObsDf = totalObsDf.drop('totim',axis=0)
totalObsDf['ObsHead'] = 0
totalObsDf['Lat'] = 0
totalObsDf['Lon'] = 0
totalObsDf.head()
SimHead ObsHead Lat Lon
HD_OBS_474820122185601 1.045374e+02 0 0 0
HD_OBS_474851122165301 1.084544e+02 0 0 0
HD_OBS_474907122194901 -1.000000e+30 0 0 0
HD_OBS_474909122164401 1.096796e+02 0 0 0
HD_OBS_474923122195601 -1.000000e+30 0 0 0
#open the observed data file
obsDf = pd.read_csv('../Txt/HobDf.csv',index_col=0)
obsDf.head()
SiteNo Latitude Longitude LatLonDatum SurfaceElevation_m SiteCode ScreenElevation_m WTElevation_m
0 12126800 47.825556 -122.254167 NAD83 96.6216 Obs_12126800 NaN NaN
1 12128100 47.884461 -122.299211 NAD83 124.9680 Obs_12128100 NaN NaN
2 12128130 47.909753 -122.297894 NAD83 97.5360 Obs_12128130 NaN NaN
3 474758122192101 47.822041 -122.298186 NAD83 119.9810 Obs_474758122192101 78.2234 NaN
4 474820122185601 47.814541 -122.306520 NAD83 107.7890 Obs_474820122185601 39.2090 88.1294
#look for the observed elevation and coordinates
for index, row in totalObsDf.iterrows():
    for index2, row2 in obsDf.iterrows():
        if index[7:] == str(row2.SiteNo):
            totalObsDf.loc[index,'ObsHead'] = row2.WTElevation_m
            totalObsDf.loc[index,'Lat'] = row2.Latitude
            totalObsDf.loc[index,'Lon'] = row2.Longitude
#compute the residual and save csv file
totalObsDf = totalObsDf[totalObsDf['SimHead']>0]
totalObsDf['Residual'] = totalObsDf['SimHead'] - totalObsDf['ObsHead']
totalObsDf.to_csv('../Txt/HobDf_Out.csv')
totalObsDf.head()
SimHead ObsHead Lat Lon Residual
HD_OBS_474820122185601 104.537359 88.129400 47.814541 -122.306520 16.407959
HD_OBS_474851122165301 108.454441 115.713800 47.813986 -122.282630 -7.259359
HD_OBS_474909122164401 109.679565 113.885000 47.818986 -122.280130 -4.205435
HD_OBS_474953122170201 108.180467 123.665778 47.830930 -122.286520 -15.485311
HD_OBS_475044122155601 108.793058 80.357000 47.845374 -122.266797 28.436058
totalObsDf
SimHead ObsHead Lat Lon Residual
HD_OBS_474820122185601 104.537359 88.129400 47.814541 -122.306520 16.407959
HD_OBS_474851122165301 108.454441 115.713800 47.813986 -122.282630 -7.259359
HD_OBS_474909122164401 109.679565 113.885000 47.818986 -122.280130 -4.205435
HD_OBS_474953122170201 108.180467 123.665778 47.830930 -122.286520 -15.485311
HD_OBS_475044122155601 108.793058 80.357000 47.845374 -122.266797 28.436058
HD_OBS_475104122191901 55.070670 69.384200 47.850930 -122.323189 -14.313530
HD_OBS_475106122170101 104.033737 114.494600 47.851485 -122.284854 -10.460863
HD_OBS_475120122162401 108.820731 123.029000 47.855374 -122.274576 -14.208269
HD_OBS_475150122192101 40.718962 41.590660 47.863707 -122.323745 -0.871699
HD_OBS_475328122181901 65.566638 43.476200 47.890929 -122.306523 22.090438
HD_OBS_475334122180401 72.714256 72.279800 47.892596 -122.302357 0.434456
#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(totalObsDf['SimHead'].min()-5, totalObsDf['SimHead'].max()+5, 100)
plt.plot(x, x, linestyle='dashed')
plt.scatter(totalObsDf['ObsHead'],totalObsDf['SimHead'], marker='o', c=totalObsDf['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()

Input files

You can download the input files from this link.

 

—————
TRANSCOM ISPFree Sigma HSE Email
Level 6 Domain Names are FREE – Register Now.

Related Posts