# Effective stress calculation from MODFLOW Groundwater Flow Model with Model Muse & Flopy – Tutorial

Effective stress theory was developed by Terzaghi in the 1920s. Based on our modeling experience we wanted to calculate the effective stress based on the results from a MODFLOW groundwater model. Finally, after 6 years from the first thought about it we came with a full deduction of the effective stress calculation based on the model geometry and an applied example for the effective stress calculation on a hillslope groundwater flow model.

The example model was developed in Modflow-Nwt and Model Muse, whereas the effective stress determination was done with scripts in Python and Hataripy (our fork of Flopy). The scripts can also generate 3D objects as VTK files of the model results, geometry and stresses that can be visualized in Paraview.

## Coding

This in the Python for the effective stress calculation and geometry generation:

``````import hataripy
import numpy as np
import matplotlib.pyplot as plt

``````#define model name, path, executables
modelWs = '../Model/'
modName = 'HillslopeModel_v1d_mf2005'
nwtPath = '../Exe/MODFLOW-NWT_64.exe'

check=False, exe_name=nwtPath)

# get a list of the model packages
mf.get_package_list()

['DIS', 'BAS6', 'OC', 'NWT', 'UPW', 'GHB', 'DRN', 'RCH']

[1.0]

(20, 57, 46)

#plot surface topography
plt.imshow(mf.modelgrid.top)

<matplotlib.image.AxesImage at 0x28f1f09ac70>``````

### Effective stress calculation

``````# calculate water table from heads
else:
waterTableArray[row,col]=-1e+20

#plot water table
plt.imshow(waterTableArray)

<matplotlib.image.AxesImage at 0x28f1f819190>``````
``````# H1 determination

H1[lay] = np.where((mf.modelgrid.top>waterTableArray)&(waterTableArray>mf.modelgrid.botm[lay]),
mf.modelgrid.top-waterTableArray,0)

# H2 determination

H2[lay] = np.where((mf.modelgrid.top>waterTableArray)&(waterTableArray>mf.modelgrid.zcellcenters[lay]),
waterTableArray-mf.modelgrid.zcellcenters[lay],0)

# HCell determination

rhoUnsat = 15 # in kN/m3 as reference value
rhoSat = 19 # in kN/m3 as reference value
rhoWater = 9.8 # in kN/m3

effectiveStress = rhoUnsat*H1 + rhoSat*H2 - rhoWater*HCell

#plot effective stresses for lay 1
plt.imshow(effectiveStress)

<matplotlib.image.AxesImage at 0x28f1f867ac0>``````
``````#plot effective stresses for row 25
plt.imshow(effectiveStress[:,25,:])``````

### Working on the VTK

``````# get information about the drain cells
drnCells = mf.drn.stress_period_data

# get information about the ghb cells
ghbCells = mf.ghb.stress_period_data``````
``````# add the arrays to the vtkObject
vtkObject = hataripy.export.vtk.Vtk3D(mf,'../vtuFiles/',verbose=True)

# add the effective stress array
``````# Create the VTKs for model output, boundary conditions and water table
vtkObject.modelMesh('modelDrn.vtu',smooth=True,cellvalues=['drn'],boundary='drn',avoidpoint=True)
vtkObject.modelMesh('modelGhb.vtu',smooth=True,cellvalues=['ghb'],boundary='ghb',avoidpoint=True)
vtkObject.waterTable('waterTable.vtu',smooth=True)``````
``````Writing vtk file: modelMesh.vtu
Number of point is 419520, Number of cells is 52440

Writing vtk file: modelDrn.vtu
Number of point is 2072, Number of cells is 259

Writing vtk file: modelGhb.vtu
Number of point is 11064, Number of cells is 1383

Writing vtk file: waterTable.vtu
Number of point is 10488, Number of cells is 2622``````

## Related Posts

#### Geospatial representation drone camera coordinates with Python and Folium – Tutorial

If you have a set of drone imagery and you want to know the location and direction of the camera, this tutorial might be interesting for you. We have done an applied example that retrieves the geospatial metadata from the

#### Evaporation and Transpiration as you have never seen before

A consistent part of root zone and surface water evaporates and returns to the atmosphere to eventually form clouds and precipitation again. The process follows

#### Are hydrogeologists / numerical modelers forbidden from cloud services?

In the semantic space the forbidden word is related to refuse to allow, and thus can create a controversy with the freedom of hydrogeologist / numerical modelers to choose whichever tool that can make the work done or the simulation

#### River-Well-Aquifer Geospatial Groundwater Flow Model with Voronoi Mesh – Tutorial

This is an applied example of a fully geospatial groundwater flow model with the MODFLOW 6 Disv discretized by vertices option. Model was constructed with Python and Flopy from a series of shapefiles and it has a progressive refinement for

Boost Inflight Internet- Business Hosting- Secure Email Account- Dropcatch Domain Names- Antisnoop Email- Free Secure Email- Cheap VOIP Calls- Free Hosting and Email- Aero Connectivity- Premium Domains for Sale- Transcom Telecom- Satphone Airtime- Mobile Plans- Free Domain Names- Organic Products- Aviation News
Transcom ISP - Transcom VOIP - Free Secure Email - Dropcatch Software - FastApn Inflight - Aero Connect - Premium Domains