Integrated Hydrology

This is a series of problems to test and demonstrate the integrated hydrology components of ATS. Integrated hydrology consists of Richards equation in the subsurface coupled to a diffusion wave approximation on the surface system. Coupling is done through assuming continuous pressure and flux across the surface-subsurface interface.

Tests demonstrate a 1D column (showing dry-to-wet and wet-to-dry transitions), a 2D “open book” problem (showing that water flows down hill) and a typical 2D “hillslope” problem comparable to the Sulis et al (WRR) benchmark problems.

Author: Ethan Coon

%matplotlib inline

import sys,os
sys.path.append(os.path.join(os.environ['ATS_SRC_DIR'],'tools', 'utils'))
sys.path.append(os.path.join(os.environ['ATS_SRC_DIR'],'docs', 'documentation', 'source', 'ats_demos'))

import ats_xdmf
import plot_column_data

import colors
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import colorbar
import pandas
import run_demos

Example 1: 1D infiltration

In the first example, we rain on a column for half a day, then turn off the rain. A continuous, constant baseflow is applied as a sink at the bottom of the column. Initially the column is partially saturated. As rainfall > baseflow, the column saturates and ponded water expresses on the surface. After rainfall stops, baseflow takes the water table below the surface again.

This demonstrates the dry-to-wet and wet-to-dry transitions.

# run the simulation
run_demos.run_demo_local('column')

  Test log file : LOGS/ats-tests-2024-04-22_12-15-32.testlog
Running ATS demo problems :
04_integrated_hydro  | .

----------------------------------------------------------------------
Regression test summary:
    Total run time: 5.01331 [s]
    Total tests : 1
    Tests run : 1
    All tests passed.

# load the ATS solution
directory = "./column.demo"
vis = ats_xdmf.VisFile(directory, time_unit='d')
vis.loadMesh(columnar=True)

vis_surf = ats_xdmf.VisFile(directory, domain='surface', time_unit='d')
vis.loadMesh(columnar=True)

# plot the solution
fig, axs = plt.subplots(1,2, figsize=(10,4))

plot_column_data.plot_subsurface(vis, ['saturation_liquid',], axs[0], None, cmap='jet')
axs[0].set_xlabel('saturation [-]')
axs[0].set_ylabel('z [m]')
axs[0].set_xlim(0.7,1.05)
axs[0].set_ylim(0,10)

plot_column_data.plot_surface(vis_surf, ['ponded_depth',], axs[1], 'ponded depth [m]','k')
axs[1].set_xlabel("time [days]")
axs[1].set_ylabel("ponded depth [m]")

plt.show()

Example 2: IH-MIP Runoff Generation

These are lowest-resolution runs (80m in x, 0.2m in z) from the Integrated Hydrology MIP #1 paper, Maxwell et al WRR, 2013, and demonstrate saturation-limited and infiltration-limited runs of integrated hydrology. Note there is an error in that paper – units of alpha in the van Genuchten curve are [1/m], not [1/cm].

# run the simulation
run_demos.run_demo_local(['saturation_limited1', 'saturation_limited2', 'infiltration_limited1', 'infiltration_limited2'])

  Test log file : LOGS/ats-tests-2024-04-22_12-15-38.testlog
Running ATS demo problems :
04_integrated_hydro  | ....

----------------------------------------------------------------------
Regression test summary:
    Total run time: 18.0182 [s]
    Total tests : 4
    Tests run : 4
    All tests passed.

# plot the hydrograph
d1 = 'saturation_limited1.demo'
d2 = 'saturation_limited2.demo'
d3 = 'infiltration_limited1.demo'
d4 = 'infiltration_limited2.demo'
fig, axs2 = plt.subplots(1,2, figsize=(10,3))

dat = pandas.read_csv(os.path.join(d1, 'surface_outlet_flux.dat'), comment='#')
axs2[0].plot(dat['time [s]'][:]/60.0, dat['surface outlet flux'][:]/55389.0*60., 'b', label='0.5m IC')
dat2 = pandas.read_csv(os.path.join(d2, 'surface_outlet_flux.dat'), comment='#')
axs2[0].plot(dat2['time [s]'][:]/60.0, dat2['surface outlet flux'][:]/55389.0*60., 'r', label='1m IC')
axs2[0].legend()
axs2[0].set_xlabel('time [min]')
axs2[0].set_ylabel('flow rate [m^3 / min]')

dat3 = pandas.read_csv(os.path.join(d3, 'surface_outlet_flux.dat'), comment='#')
axs2[1].plot(dat3['time [s]'][:]/60.0, dat3['surface outlet flux'][:]/55389.0*60., 'goldenrod', label='k_sat = 6.94e-5')
dat4 = pandas.read_csv(os.path.join(d4, 'surface_outlet_flux.dat'), comment='#')
axs2[1].plot(dat4['time [s]'][:]/60.0, dat4['surface outlet flux'][:]/55389.0*60., 'forestgreen', label='k_sat = 6.94e-6')
axs2[1].legend()
axs2[1].set_xlabel('time [min]')
axs2[1].set_ylabel('flow rate [m^3 / min]')

plt.tight_layout()
plt.show()

Example 3: Sulis et al comparison

These are lowest-resolution runs (80m in x, 0.2m in z) from the predecessor to that paper, Sulis et al 2010, which investigated the comparison across a variety of slopes as well. These demonstrate saturation-limited and infiltration-limited runs of integrated hydrology.

# run the demos.  Note this should only run two runs if you have previously run the last set of demos.
run_demos.run_demo_local(['saturation_limited1', 'saturation_limited3', 'infiltration_limited1', 'infiltration_limited3'])

  Test log file : LOGS/ats-tests-2024-04-22_12-15-56.testlog
Running ATS demo problems :
04_integrated_hydro  | ..

----------------------------------------------------------------------
Regression test summary:
    Total run time: 8.01383 [s]
    Total tests : 2
    Tests run : 2
    All tests passed.

# plot the hydrograph
d1 = 'saturation_limited1.demo'
d2 = 'saturation_limited3.demo'
d3 = 'infiltration_limited1.demo'
d4 = 'infiltration_limited3.demo'
fig, axs3 = plt.subplots(1,2, figsize=(10,3))

n_liq = 55347.3783

dat = pandas.read_csv(os.path.join(d1, 'surface_outlet_flux.dat'), comment='#')
axs3[0].plot(dat['time [s]'][:]/60.0, dat['surface outlet flux'][:]/n_liq*60., 'b', label='0.05% slope')
dat2 = pandas.read_csv(os.path.join(d2, 'surface_outlet_flux.dat'), comment='#')
axs3[0].plot(dat2['time [s]'][:]/60.0, dat2['surface outlet flux'][:]/n_liq*60., 'r', label='5% slope')
axs3[0].set_xlabel('time [min]')
axs3[0].set_ylabel('flow rate [m^3 / min]')
axs3[0].legend()

dat3 = pandas.read_csv(os.path.join(d3, 'surface_outlet_flux.dat'), comment='#')
axs3[1].plot(dat3['time [s]'][:]/60.0, dat3['surface outlet flux'][:]/n_liq*60., 'goldenrod', label='0.05% slope')
dat4 = pandas.read_csv(os.path.join(d4, 'surface_outlet_flux.dat'), comment='#')
axs3[1].plot(dat4['time [s]'][:]/60.0, dat4['surface outlet flux'][:]/n_liq*60., 'forestgreen', label='5% slope')
axs3[1].legend()
axs3[1].set_xlabel('time [min]')
axs3[1].set_ylabel('flow rate [m^3 / min]')

plt.tight_layout()
plt.show()

Example 4: V-catchment transect

In the next example, we rain on a transect of a V-catchment in 2D, allowing the water to pond. This simply demonstrates that water runs downhill (in a coupled environment) and pools on the surface. The plots show saturation (in colors) in the subsurface and the ponded depth (a line above the catchment).

At time zero, a water table is established. Infiltration occurs through the 30 hour plot. By 60 hours, the domain is saturated and water is ponding on the surface.

# run the demo
run_demos.run_demo_local('transect')

  Test log file : LOGS/ats-tests-2024-04-22_12-16-04.testlog
Running ATS demo problems :
04_integrated_hydro  | .

----------------------------------------------------------------------
Regression test summary:
    Total run time: 6.01323 [s]
    Total tests : 1
    Tests run : 1
    All tests passed.

# plot the saturation at a few times
directory = "transect.demo"

vis = ats_xdmf.VisFile(directory, time_unit='d')
vis.filterIndices(range(0,21,4))
vis.loadMeshPolygons()

vis_surf = ats_xdmf.VisFile(directory, domain='surface', time_unit='d')
vis_surf.filterIndices(range(0,21,4))
vis_surf.loadMesh(order=['x','z'])

fig, axs = plt.subplots(2,3,figsize=(12,8))

def plot(i, ax):
    cax = colorbar.make_axes(ax, 'right')[0]
    cax.set_ylabel("saturation [-]")

    sat = vis.get("saturation_liquid", vis.cycles[i])
    poly = vis.getMeshPolygons(cmap='jet_r', linewidth=1)
    poly.set_array(sat)
    poly.set_clim(0.4,1)
    ax.add_collection(poly)

    elev = vis_surf.get('surface-elevation', vis.cycles[i])
    pd = vis_surf.get('surface-ponded_depth', vis.cycles[i])

    ax.plot(vis_surf.centroids[:,0], elev+pd, 'b', linewidth=2)
    ax.set_title(f"{3*i} hours")

plot(0, axs[0,0])
plot(1, axs[0,1])
plot(2, axs[0,2])
plot(3, axs[1,0])
plot(4, axs[1,1])
plot(5, axs[1,2])

Example 5: Superslab problem

The last demo is based off of the superslab problem from Kollet et all 2018 IH-MIP Benchmark 2 paper. The goal of this demo is to show how to extend the above integrated hydrology problems into cases with multiple materials or soil types.

Note this run is a bit bigger (10k grid cells) and is natively run on 2 cores. It takes about 20 minutes to run on a laptop with debug code.

# run the demo
run_demos.run_demo_local('superslab')

  Test log file : LOGS/ats-tests-2024-04-22_12-16-11.testlog
Running ATS demo problems :
04_integrated_hydro  | .

----------------------------------------------------------------------
Regression test summary:
    Total run time: 103.022 [s]
    Total tests : 1
    Tests run : 1
    All tests passed.

# plot the saturation at a few times
directory = "superslab.demo"

vis = ats_xdmf.VisFile(directory, time_unit='d')
vis.filterIndices(list(range(0,11,2)))
vis.loadMeshPolygons()

vis_surf = ats_xdmf.VisFile(directory, domain='surface', time_unit='d')
vis_surf.filterIndices(list(range(0,11,2)))
vis_surf.loadMesh(order=['x','z'])

fig, axs = plt.subplots(2,3,figsize=(12,8))

def plot(i, ax):
    cax = colorbar.make_axes(ax, 'right')[0]
    cax.set_ylabel("saturation [-]")

    sat = vis.get("saturation_liquid", vis.cycles[i])
    poly = vis.getMeshPolygons(cmap='jet_r', linewidth=0)
    poly.set_array(sat)
    poly.set_clim(0.4,1)
    ax.add_collection(poly)

    elev = vis_surf.get('surface-elevation', vis.cycles[i])
    pd = vis_surf.get('surface-ponded_depth', vis.cycles[i])

    ax.plot(vis_surf.centroids[:,0], elev+pd, 'b', linewidth=2)
    ax.set_title(f"{3*i} hours")

plot(0, axs[0,0])
plot(1, axs[0,1])
plot(2, axs[0,2])
plot(3, axs[1,0])
plot(4, axs[1,1])
plot(5, axs[1,2])