Lesson 5 - Land Surface Experiments

Overview

In Lessons 4 we experimented with different initial and boundary conditions (precipitation). In this lesson, we will experiment with different terrain physics options to evaluate impacts on streamflow. Click here for reference material on terrain routing and groundwater in WRF-Hydro.

NOTE: If you have not completed Lessons 1 through 4, please stop and do so now.

Software and conventions

The easiest way to run these lessons is via the wrfhydro/training Docker container, which has all software dependencies and data pre-installed.

For a complete description of the software environment used for this training please see Getting started.

You may either execute commands by running each cell of this notebook. Alternatively, you may open a terminal in Jupyter Lab by selecting New -> Terminal in your Home tab of Jupyter Lab and input the commands manually if you prefer. You can also use your own terminal by logging into the container with the command docker exec -it wrf-hydro-training bash

All paths used in this lesson assume that the lesson materials are located under your home directory in a folder named wrf-hydro-training. If your materials are located in another directory, you will not be able to run the commands in this notebook inside Jupyter and will need to type them manually in your terminal session.

Compiling WRF-Hydro

Below are the commands to compile WRF-Hydro. We are doing a quick short-cut to edit the setEnvar.sh to make sure HYDRO_D and SPATIAL_SOIL are both active. See Lesson 1 if you have any questions on these settings.

NOTE: You only need to do these steps if you do not already have a compiled binary. If you do, you can skip the compile steps.

In [1]:
%%bash
# Change to the trunk/NDHMS directory and configure for gfort
cd ~/wrf-hydro-training/wrf_hydro_nwm_public/trunk/NDHMS/; ./configure 2

# Make a copy of the template environment variable file, setEnvars.sh
cp ~/wrf-hydro-training/wrf_hydro_nwm_public/trunk/NDHMS/template/setEnvar.sh \
~/wrf-hydro-training/wrf_hydro_nwm_public/trunk/NDHMS/setEnvar.sh

# Edit setEnvar.sh to turn on debug (HYDRO_D=1) and spatial soils (SPATIAL_SOIL=1)
sed -i -e 's#export HYDRO_D=0#export HYDRO_D=1#g' ~/wrf-hydro-training/wrf_hydro_nwm_public/trunk/NDHMS/setEnvar.sh
sed -i -e 's#export SPATIAL_SOIL=0#export SPATIAL_SOIL=1#g' ~/wrf-hydro-training/wrf_hydro_nwm_public/trunk/NDHMS/setEnvar.sh

# Compile the code
./compile_offline_NoahMP.sh setEnvar.sh >> compile.log 2>&1

# Check to make sure it completed successfully and confirm the correct compile-time options were set.
tail -13 ~/wrf-hydro-training/wrf_hydro_nwm_public/trunk/NDHMS/compile.log
Configured: gfort
*****************************************************************
Make was successful

*****************************************************************
The environment variables used in the compile:
HYDRO_D=1
NCEP_WCOSS=0
NETCDF=/usr/local
SPATIAL_SOIL=1
WRFIO_NCD_LARGE_FILE_SUPPORT=1
WRF_HYDRO=1
WRF_HYDRO_NUDGING=0
WRF_HYDRO_RAPID=0

Experiment with Terrain Routing Options

Background

In most global and regional land surface models, model cells are simulated as independent columns and "runoff" is interpreted as excess surface and subsurface water that is "scraped" from each cell after vertical infiltration and drainage processes complete. WRF-Hydro adds functionality to allow this excess water to instead move laterally around the land surface based on topographic and head gradients. WRF-Hydro can explicitly represent overland flow processes, where infiltration and saturation excess water propagates over the surface through a diffusive wave formulation. In addition, WRF-Hydro includes shallow subsurface flow, allowing saturated soil water to move from cell to cell through a Boussinesq approximation. For more detailed information refer to WRF-Hydro V5 Technical Description.

Objective

Experiment with lateral flow physics options, then compare hydrographs with and without terrain routing processes active.

Step 1: Create and run first baseline simulation

As in Lesson 4, we will create a template directory to use for the terrain physics experiments.

Step 1a: Create a template run directory

In [2]:
%%bash
# Make a new directory for our baseline simulation
mkdir -p ~/wrf-hydro-training/output/lesson5/run_gridded_template

# Copy our model files to the simulation directory
cp ~/wrf-hydro-training/wrf_hydro_nwm_public/trunk/NDHMS/Run/*.TBL \
~/wrf-hydro-training/output/lesson5/run_gridded_template
cp ~/wrf-hydro-training/wrf_hydro_nwm_public/trunk/NDHMS/Run/wrf_hydro.exe \
~/wrf-hydro-training/output/lesson5/run_gridded_template
cp ~/wrf-hydro-training/example_case/Gridded/namelist.hrldas \
~/wrf-hydro-training/output/lesson5/run_gridded_template
cp ~/wrf-hydro-training/example_case/Gridded/hydro.namelist \
~/wrf-hydro-training/output/lesson5/run_gridded_template

# Create symbolic links to large domain files
cp -as $HOME/wrf-hydro-training/example_case/FORCING \
~/wrf-hydro-training/output/lesson5/run_gridded_template
cp -as $HOME/wrf-hydro-training/example_case/Gridded/RESTART \
~/wrf-hydro-training/output/lesson5/run_gridded_template
cp -as $HOME/wrf-hydro-training/example_case/Gridded/DOMAIN \
~/wrf-hydro-training/output/lesson5/run_gridded_template

Step 1b: Now make a copy of the template directory and run the baseline simulation

First, we remind ourselves what we are running as a baseline. Note in particular that OVRTSWCRT=1 and SUBRTSWCRT=1, meaning that both overland and subsurface flow routing options are active.

In [3]:
%%bash
# Make a copy of the template directory for the baseline run
cp -r ~/wrf-hydro-training/output/lesson5/run_gridded_template \
~/wrf-hydro-training/output/lesson5/run_gridded_baseline
# View the hydro.namelist settings
cat ~/wrf-hydro-training/output/lesson5/run_gridded_baseline/hydro.namelist
&HYDRO_nlist
!!!! ---------------------- SYSTEM COUPLING ----------------------- !!!!

! Specify what is being coupled:  1=HRLDAS (offline Noah-LSM), 2=WRF, 3=NASA/LIS, 4=CLM
sys_cpl = 1

!!!! ------------------- MODEL INPUT DATA FILES ------------------- !!!!

! Specify land surface model gridded input data file (e.g.: "geo_em.d01.nc")
GEO_STATIC_FLNM = "./DOMAIN/geo_em.d01.nc"

! Specify the high-resolution routing terrain input data file (e.g.: "Fulldom_hires.nc")
GEO_FINEGRID_FLNM = "./DOMAIN/Fulldom_hires.nc"

! Specify the spatial hydro parameters file (e.g.: "hydro2dtbl.nc")
! If you specify a filename and the file does not exist, it will be created for you.
HYDROTBL_F = "./DOMAIN/hydro2dtbl.nc"

! Specify spatial metadata file for land surface grid. (e.g.: "GEOGRID_LDASOUT_Spatial_Metadata.nc")
LAND_SPATIAL_META_FLNM = "./DOMAIN/GEOGRID_LDASOUT_Spatial_Metadata.nc"

! Specify the name of the restart file if starting from restart...comment out with '!' if not...
RESTART_FILE  = 'RESTART/HYDRO_RST.2017-11-11_00:00_DOMAIN1'

!!!! --------------------- MODEL SETUP OPTIONS -------------------- !!!!

! Specify the domain or nest number identifier...(integer)
IGRID = 1

! Specify the restart file write frequency...(minutes)
! A value of -99999 will output restarts on the first day of the month only.
rst_dt = 14400

! Reset the LSM soil states from the high-res routing restart file (1=overwrite, 0=no overwrite)
! NOTE: Only turn this option on if overland or subsurface rotuing is active!
rst_typ = 1

! Restart file format control
rst_bi_in = 0       !0: use netcdf input restart file (default)
                    !1: use parallel io for reading multiple restart files, 1 per core
rst_bi_out = 0      !0: use netcdf output restart file (default)
                    !1: use parallel io for outputting multiple restart files, 1 per core

! Restart switch to set restart accumulation variables to 0 (0=no reset, 1=yes reset to 0.0)
RSTRT_SWC = 0

! Specify baseflow/bucket model initialization...(0=cold start from table, 1=restart file)
GW_RESTART = 1

!!!! -------------------- MODEL OUTPUT CONTROL -------------------- !!!!

! Specify the output file write frequency...(minutes)
out_dt = 60

! Specify the number of output times to be contained within each output history file...(integer)
!   SET = 1 WHEN RUNNING CHANNEL ROUTING ONLY/CALIBRATION SIMS!!!
!   SET = 1 WHEN RUNNING COUPLED TO WRF!!!
SPLIT_OUTPUT_COUNT = 1

! Specify the minimum stream order to output to netcdf point file...(integer)
! Note: lower value of stream order produces more output.
order_to_write = 1

! Flag to turn on/off new I/O routines: 0 = deprecated output routines (use when running with Noah LSM),
! 1 = with scale/offset/compression, ! 2 = with scale/offset/NO compression,
! 3 = compression only, 4 = no scale/offset/compression (default)
io_form_outputs = 4

! Realtime run configuration option:
! 0=all (default), 1=analysis, 2=short-range, 3=medium-range, 4=long-range, 5=retrospective,
! 6=diagnostic (includes all of 1-4 outputs combined)
io_config_outputs = 0

! Option to write output files at time 0 (restart cold start time): 0=no, 1=yes (default)
t0OutputFlag = 1

! Options to output channel & bucket influxes. Only active for UDMP_OPT=1.
! Nonzero choice requires that out_dt above matches NOAH_TIMESTEP in namelist.hrldas.
! 0=None (default), 1=channel influxes (qSfcLatRunoff, qBucket)
! 2=channel+bucket fluxes    (qSfcLatRunoff, qBucket, qBtmVertRunoff_toBucket)
! 3=channel accumulations    (accSfcLatRunoff, accBucket) *** NOT TESTED ***
output_channelBucket_influx = 0

! Output netcdf file control
CHRTOUT_DOMAIN = 0           ! Netcdf point timeseries output at all channel points (1d)
                             !      0 = no output, 1 = output
CHANOBS_DOMAIN = 1           ! Netcdf point timeseries at forecast points or gage points (defined in Routelink)
                             !      0 = no output, 1 = output at forecast points or gage points.
CHRTOUT_GRID = 0             ! Netcdf grid of channel streamflow values (2d)
                             !      0 = no output, 1 = output
                             !      NOTE: Not available with reach-based routing
LSMOUT_DOMAIN = 0            ! Netcdf grid of variables passed between LSM and routing components (2d)
                             !      0 = no output, 1 = output
                             !      NOTE: No scale_factor/add_offset available
RTOUT_DOMAIN = 1             ! Netcdf grid of terrain routing variables on routing grid (2d)
                             !      0 = no output, 1 = output
output_gw = 1                ! Netcdf GW output
                             !      0 = no output, 1 = output
outlake  = 0                 ! Netcdf grid of lake values (1d)
                             !      0 = no output, 1 = output
frxst_pts_out = 0            ! ASCII text file of forecast points or gage points (defined in Routelink)
                             !      0 = no output, 1 = output

!!!! ------------ PHYSICS OPTIONS AND RELATED SETTINGS ------------ !!!!

! Specify the number of soil layers (integer) and the depth of the bottom of each layer... (meters)
! Notes: In Version 1 of WRF-Hydro these must be the same as in the namelist.input file.
!      Future versions will permit this to be different.
NSOIL=4
ZSOIL8(1) = -0.10
ZSOIL8(2) = -0.40
ZSOIL8(3) = -1.00
ZSOIL8(4) = -2.00

! Specify the grid spacing of the terrain routing grid...(meters)
DXRT = 100.0

! Specify the integer multiple between the land model grid and the terrain routing grid...(integer)
AGGFACTRT = 10

! Specify the channel routing model timestep...(seconds)
DTRT_CH = 10

! Specify the terrain routing model timestep...(seconds)
DTRT_TER = 10

! Switch to activate subsurface routing...(0=no, 1=yes)
SUBRTSWCRT = 1

! Switch to activate surface overland flow routing...(0=no, 1=yes)
OVRTSWCRT = 1

! Specify overland flow routing option: 1=Seepest Descent (D8) 2=CASC2D (not active)
! NOTE: Currently subsurface flow is only steepest descent
rt_option = 1

! Switch to activate channel routing...(0=no, 1=yes)
CHANRTSWCRT = 1

! Specify channel routing option: 1=Muskingam-reach, 2=Musk.-Cunge-reach, 3=Diff.Wave-gridded
channel_option = 3

! Specify the reach file for reach-based routing options (e.g.: "Route_Link.nc")
!route_link_f = "./DOMAIN/Route_Link.nc"

! If using channel_option=2, activate the compound channel formulation? (Default=.FALSE.)
! This option is currently only supported if using reach-based routing with UDMP=1.
compound_channel = .FALSE.

! Specify the lake parameter file (e.g.: "LAKEPARM.nc").
! Note REQUIRED if lakes are on.
!route_lake_f = "./DOMAIN/LAKEPARM.nc"

! Switch to activate baseflow bucket model...(0=none, 1=exp. bucket, 2=pass-through)
GWBASESWCRT = 1

! Groundwater/baseflow 2d mask specified on land surface model grid (e.g.: "GWBASINS.nc")
!Note: Only required if baseflow  model is active (1 or 2) and UDMP_OPT=0.
gwbasmskfil = "./DOMAIN/GWBASINS.nc"

! Groundwater bucket parameter file (e.g.: "GWBUCKPARM.nc")
GWBUCKPARM_file = "./DOMAIN/GWBUCKPARM.nc"

! User defined mapping, such NHDPlus: 0=no (default), 1=yes
UDMP_OPT = 0

! If on, specify the user-defined mapping file (e.g.: "spatialweights.nc")
!udmap_file = "./DOMAIN/spatialweights.nc"

/

&NUDGING_nlist

! Path to the "timeslice" observation files.
timeSlicePath = "./nudgingTimeSliceObs/"

nudgingParamFile = "DOMAIN/nudgingParams.nc"

! Nudging restart file = "nudgingLastObsFile"
! nudgingLastObsFile defaults to '', which will look for nudgingLastObs.YYYY-mm-dd_HH:MM:SS.nc
!   **AT THE INITALIZATION TIME OF THE RUN**. Set to a missing file to use no restart.
!nudgingLastObsFile = '/a/nonexistent/file/gives/nudging/cold/start'

!! Parallel input of nudging timeslice observation files?
readTimesliceParallel = .TRUE.

! temporalPersistence defaults to true, only runs if necessary params present.
temporalPersistence = .FALSE.

! The total number of last (obs, modeled) pairs to save in nudgingLastObs for
! removal of bias. This is the maximum array length. (This option is active when persistBias=FALSE)
! (Default=960=10days @15min obs resolution, if all the obs are present and longer if not.)
nLastObs = 960

! If using temporalPersistence the last observation persists by default.
! This option instead persists the bias after the last observation.
persistBias = .FALSE.

! AnA (FALSE)  vs Forecast (TRUE) bias persistence.
! If persistBias: Does the window for calculating the bias end at
! model init time (=t0)?
! FALSE = window ends at model time (moving),
! TRUE = window ends at init=t0(fcst) time.
! (If commented out, Default=FALSE)
! Note: Perfect restart tests require this option to be .FALSE.
biasWindowBeforeT0 = .FALSE.

! If persistBias: Only use this many last (obs, modeled) pairs. (If Commented out, Default=-1*nLastObs)
! > 0: apply an age-based filter, units=hours.
! = 0: apply no additional filter, use all available/usable obs.
! < 0: apply an count-based filter, units=count
maxAgePairsBiasPersist = -960

! If persistBias: The minimum number of last (obs, modeled) pairs, with age less than
! maxAgePairsBiasPersist, required to apply a bias correction. (default=8)
minNumPairsBiasPersist = 8

! If persistBias: give more weight to observations closer in time? (default=FALSE)
invDistTimeWeightBias = .TRUE.

! If persistBias: "No constructive interference in bias correction?", Reduce the bias adjustment
! when the model and the bias adjustment have the same sign relative to the modeled flow at t0?
! (default=FALSE)
! Note: Perfect restart tests require this option to be .FALSE.
noConstInterfBias = .FALSE.

/

Now run the baseline model simulation as-is.

In [4]:
%%bash
# Run the simulation
cd ~/wrf-hydro-training/output/lesson5/run_gridded_baseline
mpirun -np 2 ./wrf_hydro.exe >> run.log 2>&1

Check that the simulation finished successfully

In [5]:
%%bash
tail -1 ~/wrf-hydro-training/output/lesson5/run_gridded_baseline/diag_hydro.00000
 The model finished successfully.......

Step 2: Create and run the "overland routing off" simulation

We will run the same experiment as the above baseline run, but with the overland routing module turned off.

Step 2a: Create a new run directory

We first create a run_overland_routing_off directory from the template.

In [6]:
%%bash
cp -r ~/wrf-hydro-training/output/lesson5/run_gridded_template \
~/wrf-hydro-training/output/lesson5/run_overland_routing_off

We will make one modification to hydro.namelist to turn off the overland routing modules.

Step 2b: Edit the hydro.namelist file

For this experiment, we will set the OVRTSWCRT physics option to 0, which deactivates overland terrain routing modules.

! Switch to activate surface overland flow routing...(0=no, 1=yes)
OVRTSWCRT = 0

Step 2c: Run the simulation and make sure it finishes successfully

In [7]:
%%bash
cd ~/wrf-hydro-training/output/lesson5/run_overland_routing_off
mpirun -np 2 ./wrf_hydro.exe >> run.log 2>&1

Check that the simulation finished successfully

In [8]:
%%bash
tail -1 ~/wrf-hydro-training/output/lesson5/run_overland_routing_off/diag_hydro.00000
 The model finished successfully.......

Step 3: Create and run the "no terrain routing" simulation

Now we run an experiment with all terrain routing options turned off. This will approximate a standard column physics configuration.

Step 3a: Create a new run directory

We first create a run_no_terrain_routing directory from the template.

In [9]:
%%bash
cp -r ~/wrf-hydro-training/output/lesson5/run_gridded_template \
~/wrf-hydro-training/output/lesson5/run_no_terrain_routing

Step3b: Edit the hydro.namelist file We will make three modifications to hydro.namelist to turn off all of the terrain routing modules.

For this experiment, we will set the OVRTSWCRT and SUBRTSWCRT physics options to 0, which deactivates both overland and subsurface terrain routing modules.

! Switch to activate subsurface routing...(0=no, 1=yes)
SUBRTSWCRT = 0
! Switch to activate surface overland flow routing...(0=no, 1=yes)
OVRTSWCRT = 0

Since we are deactivating all terrain routing physics, we also need to tell the model not to overwrite initial conditions with the routing grids using the rst_typ option. Deactivate this function by setting rst_typ to 0.

! Reset the LSM soil states from the high-res routing restart file (1=overwrite, 0=no overwrite)
! NOTE: Only turn this option on if overland or subsurface rotuing is active!
rst_typ = 0

Step 3c: Run the simulation

In [10]:
%%bash
cd ~/wrf-hydro-training/output/lesson5/run_no_terrain_routing
mpirun -np 2 ./wrf_hydro.exe >> run.log 2>&1

Check that the simulation finished successfully

In [11]:
%%bash
tail -1 ~/wrf-hydro-training/output/lesson5/run_no_terrain_routing/diag_hydro.00000
 The model finished successfully.......

Results

We will now look at the differences in streamflow between our baseline run with all terrain routing options on and our experiments with different terrain routing options turned off.

We will use Python and the xarray library to load the data and plot hydrographs. For an intro to these tools, please see Lesson 3.

Load the xarray python package

In [12]:
# Load the xarray package
%matplotlib inline
import xarray as xr
import matplotlib.pyplot as plt
import pandas as pd

Load the CHANOBS streamflow datasets

We are going to use the CHANOBS files (vs. CHRTOUT files) because it will limit the number of reaches to only those which have a gage.

In [13]:
chanobs_control = xr.open_mfdataset('/home/docker/wrf-hydro-training/output/lesson5/run_gridded_baseline/*CHANOBS*',
                            combine='by_coords')
chanobs_ov_off = xr.open_mfdataset('/home/docker/wrf-hydro-training/output/lesson5/run_overland_routing_off/*CHANOBS*',
                            combine='by_coords')
chanobs_terrain_off = xr.open_mfdataset('/home/docker/wrf-hydro-training/output/lesson5/run_no_terrain_routing/*CHANOBS*',
                            combine='by_coords')
obs = pd.read_csv('/home/docker/wrf-hydro-training/example_case/obs_event.csv',dtype=str)
obs['dateTime'] = pd.to_datetime(obs['DATE'])
obs['streamflow_cms'] = pd.to_numeric(obs['streamflow_cms'])

Plot the hydrographs

In [14]:
fig, axes = plt.subplots(ncols=1,figsize=(12, 6))
plt.suptitle('Hydrographs for terrain routing physics options',fontsize=24)
chanobs_control.sel(feature_id = 3).streamflow.plot(label='All Terrain Routing On',
                                                        color='black',
                                                        linestyle='--')
chanobs_ov_off.sel(feature_id = 3).streamflow.plot(label='Overland Routing Off',
                                                        color='blue',
                                                        linestyle='-')
chanobs_terrain_off.sel(feature_id = 3).streamflow.plot(label='No Terrain Routing',
                                                        color='red',
                                                        linestyle='-')

obs[obs['Location'] == 'Desague'].plot(x='dateTime',
                                       y='streamflow_cms',
                                       ax=axes,
                                       label='Observed',
                                       color='grey')
plt.ylim(0,16)
plt.legend()
plt.show()
/home/docker/miniconda3/lib/python3.7/site-packages/pandas/plotting/_matplotlib/converter.py:102: FutureWarning: Using an implicitly registered datetime converter for a matplotlib plotting method. The converter was registered by pandas on import. Future versions of pandas will require you to explicitly register matplotlib converters.

To register the converters:
	>>> from pandas.plotting import register_matplotlib_converters
	>>> register_matplotlib_converters()
  warnings.warn(msg, FutureWarning)

Load the LDASOUT land model output files

In our simulations, we specified land surface model output (LDASOUT) files every hour. We will read those in using xarray.

In [15]:
ldasout_terrain_on = xr.open_mfdataset('/home/docker/wrf-hydro-training/output/lesson5/run_gridded_baseline/*.LDASOUT*',
                            combine='by_coords')
ldasout_terrain_off = xr.open_mfdataset('/home/docker/wrf-hydro-training/output/lesson5/run_no_terrain_routing/*.LDASOUT*',
                            combine='by_coords')

Plot the soil moisture states

Among other variables, the LDASOUT files include soil moisture states (SOIL_M) in each of the 4 soil layers. We will select and plot soil moisture in the top soil layer from each simulation at 2011-08-29 00Z. For reference, we also plot the soil types from the model input geogrid file.

In [16]:
# Select data for 2017-11-16 12Z at surface layer
ldasout_terrain_on_time = ldasout_terrain_on.sel(time = '2017-11-16T12:00:00').sel(soil_layers_stag = 0)
ldasout_terrain_off_time = ldasout_terrain_off.sel(time = '2017-11-16T12:00:00').sel(soil_layers_stag = 0)

# Plot the top layer soil moisture states
fig, axes = plt.subplots(ncols=2,figsize=(12, 6))
plt.suptitle('Surface soil moisture states for terrain routing on and off',fontsize=24)
ldasout_terrain_on_time.SOIL_M.plot(ax=axes[0],vmin=0.30,vmax=0.45,cmap='viridis_r')
axes[0].set_title('Terrain routing on')
ldasout_terrain_off_time.SOIL_M.plot(ax=axes[1],vmin=0.30,vmax=0.45,cmap='viridis_r')
axes[1].set_title('Terrain routing off')
plt.show()
In [17]:
# Load the geogrid dataset for reference
geogrid = xr.open_dataset('~/wrf-hydro-training/output/lesson5/run_gridded_baseline/DOMAIN/geo_em.d01.nc')
In [18]:
# Plot the dominant soil type
fig, axes = plt.subplots(figsize=(6, 6))
geogrid.SCT_DOM.plot(cmap="tab20")
axes.set_title('Soil Type')
plt.show()

The top-level soil moisture state shows significant differences with and without terrain routing active. Overall, we see much more heterogeneity in the "terrain routing on" run (above, left) due to lateral redistribution of water. In the "terrain routing off" run (above, right), soil moisture distribution is largely controlled by vertical soil properties (e.g., saturated hydraulic conductivity, porosity), and you distinctly see the soil type pattern in the soil moisture states.

Plot the mean soil moisture time series

In [19]:
# Calculate the mean bottom-layer soil moisture across the domain
smois_terron_avg = ldasout_terrain_on.SOIL_M.sel(soil_layers_stag = 3).mean(dim=('y','x'))
smois_terroff_avg = ldasout_terrain_off.SOIL_M.sel(soil_layers_stag = 3).mean(dim=('y','x')) 

# Plot the soil moisture time series
fig, axes = plt.subplots(ncols=1,figsize=(12, 6))
plt.suptitle('Average Soil Moisture: Bottom Layer',fontsize=24)
smois_terron_avg.plot(label='Terrain Routing On', color='black', linestyle='-')
smois_terroff_avg.plot(label='Terrain Routing Off', color='red', linestyle='-')
plt.legend()
plt.show()

Load the GWOUT groundwater bucket model output files

In [20]:
gwbucket_terrain_on = xr.open_mfdataset('/home/docker/wrf-hydro-training/output/lesson5/run_gridded_baseline/*GWOUT*',
                            combine='by_coords')
gwbucket_terrain_off = xr.open_mfdataset('/home/docker/wrf-hydro-training/output/lesson5/run_no_terrain_routing/*GWOUT*',
                            combine='by_coords')
gwbucket_ov_off = xr.open_mfdataset('/home/docker/wrf-hydro-training/output/lesson5/run_overland_routing_off/*GWOUT*',
                            combine='by_coords')

Plot the mean groundwater bucket level time series

Note that groundwater bucket "basins" vary in size, so a simple mean is not a good way to track mass. We do it here for simplicity and to get a general sense of the bucket behavior across the domain.

In [21]:
# Calculate the mean bucket level across the domain
gwlevel_terron_avg = gwbucket_terrain_on.depth.mean(dim=('feature_id'))
gwlevel_terroff_avg = gwbucket_terrain_off.depth.mean(dim=('feature_id')) 

# Plot the bucket level time series
fig, axes = plt.subplots(ncols=1,figsize=(12, 6))
plt.suptitle('Average Groundwater Bucket Level',fontsize=24)
gwlevel_terron_avg.plot(label='Terrain Routing On', color='black', linestyle='-')
gwlevel_terroff_avg.plot(label='Terrain Routing Off', color='red', linestyle='-')
plt.legend()
plt.show()

IMPORTANT NOTE: This particular test case has been calibrated to the full terrain routing physics options, so will obviously be sub-optimal for a different physics configuration. You should adjust your parameters to your particular model configuration. You should also re-spinup the model to give it time to adjust to the new physics. We have done neither of those steps in this short example.

Experiment with Modified Parameters

There are a number of key parameters that impact water partitioning, storage, and movement through the model system. We have pulled many of the most important model parameters into NetCDF files to ease parameter display and manipulation, as well as to allow the parameters to vary independently in space. Key terrain routing parameter files include:

  • soil_properties.nc - NoahMP soil and vegetation properties (LSM grid)
  • hydro2dtbl.nc - Lateral routing model soil and surface parameters (LSM grid)
  • Fulldom_hires.nc - Lateral routing model high-res parameters (routing grid)
  • GWBUCKPARM.nc - Groundwater baseflow bucket model parameters (groundwater basin objects)

In this lesson, we will manipulate parameters in the soil_properties.nc (refkdt, dksat, bexp, slope) file.

Step 1: Create a new template directory and and run default parameter test case

As in the first section, we will make a new simulation directory and use this as a template for creating multiple new simulation directories.

Step 1a: Create a new template directory for the parameter experiments.

In [22]:
%%bash
# Make a new directory for our default baseline simulation
mkdir -p ~/wrf-hydro-training/output/lesson5/run_default_template

# Copy our model files to the simulation directory
cp ~/wrf-hydro-training/wrf_hydro_nwm_public/trunk/NDHMS/Run/*.TBL \
~/wrf-hydro-training/output/lesson5/run_default_template
cp ~/wrf-hydro-training/wrf_hydro_nwm_public/trunk/NDHMS/Run/wrf_hydro.exe \
~/wrf-hydro-training/output/lesson5/run_default_template
cp ~/wrf-hydro-training/output/lesson5/run_gridded_template/namelist.hrldas \
~/wrf-hydro-training/output/lesson5/run_default_template
cp ~/wrf-hydro-training/output/lesson5/run_gridded_template/hydro.namelist \
~/wrf-hydro-training/output/lesson5/run_default_template

# Create symbolic links to large domain files
cp -as $HOME/wrf-hydro-training/example_case/FORCING \
~/wrf-hydro-training/output/lesson5/run_default_template
cp -as $HOME/wrf-hydro-training/example_case/Gridded/RESTART \
~/wrf-hydro-training/output/lesson5/run_default_template
cp -as $HOME/wrf-hydro-training/example_case/Gridded/DOMAIN \
~/wrf-hydro-training/output/lesson5/run_default_template

Step 1b: Setup and run the default parameter baseline experiment

We now make a new run directory to run the baseline simulation with default (uncalibrated) parameters.

In [23]:
%%bash
cp -r ~/wrf-hydro-training/output/lesson5/run_default_template \
~/wrf-hydro-training/output/lesson5/run_parameter_baseline

Launch the baseline run and make sure it completes successfully.

In [24]:
%%bash
cd ~/wrf-hydro-training/output/lesson5/run_parameter_baseline
mpirun -np 2 ./wrf_hydro.exe >> run.log 2>&1
In [25]:
%%bash
tail -1 ~/wrf-hydro-training/output/lesson5/run_parameter_baseline/diag_hydro.00000
 The model finished successfully.......

Step 1c: Plot the hydrograph for the baseline (default) run

We want to take a quick look at the hydrograph for the default parameter run so we see what behavior we might want to adjust.

In [26]:
# Load the xarray package
%matplotlib inline
import xarray as xr
import matplotlib.pyplot as plt
import pandas as pd

# Pull the channel output files into xarray objects 
chanobs_control = xr.open_mfdataset('/home/docker/wrf-hydro-training/output/lesson5/run_parameter_baseline/*CHANOBS*',
                            combine='by_coords')
obs = pd.read_csv('/home/docker/wrf-hydro-training/example_case/obs_event.csv',dtype=str)
obs['dateTime'] = pd.to_datetime(obs['DATE'])
obs['streamflow_cms'] = pd.to_numeric(obs['streamflow_cms'])

# Plot the baseline hydrograph
fig, axes = plt.subplots(ncols=1,figsize=(12, 6))
plt.suptitle('Hydrographs for parameter experiment',fontsize=24)
chanobs_control.sel(feature_id = 3).streamflow.plot(label='Control',
                                                        color='black',
                                                        linestyle='--')
obs[obs['Location'] == 'Desague'].plot(x='dateTime',
                                       y='streamflow_cms',
                                       ax=axes,
                                       label='Observed',
                                       color='grey')
plt.ylim(0,30)
plt.legend()
plt.show()

Step 2: Modify NoahMP parameters using NCO tools

We will create a new simulation directory for our parameter manipulation experiment. All of the parameter file edits will be done using NCO. NCO (NetCDF Operators, http://nco.sourceforge.net/) is a set of useful utilities to manipulate NetCDF files.

Step 2a: Setup the parameter experiment run directory

First make a new run directory where we can start modifying parameter files.

In [27]:
%%bash
cp -r ~/wrf-hydro-training/output/lesson5/run_default_template \
~/wrf-hydro-training/output/lesson5/run_parameter_mods

Step 2b: Use the NCO command ncap2 to modify the refkdt parameter values in the soil_properties.nc file.

REFKDT

There are a number of model parameters that affect lateral terrain flow. One important Noah/NoahMP parameter that we commonly adjust when activating terrain routing is refkdt. Refkdt controls how easily precipitation reaching the surface infiltrates into the soil column vs. stays on the surface where it will be "scraped" off as surface runoff. Higher values of refkdt lead to more infiltration and less surface (fast) runoff. This tunable parameter can be set to a relatively high value (e.g., 3.0) suitable for running the column land surface model only. When activating terrain routing to explicitly model these processes, we often reduce this parameter. In addition, if you are calling the land surface model on a small timestep (e.g., seconds to minutes), you may want to reduce this parameter to compensate for the more frequent calls to the vertical infiltration scheme.

Here we will use the NCO command ncap2 to decrease the refkdt value.

First, we check the current parameter values using ncdump.

In [28]:
%%bash
ncdump -v refkdt ~/wrf-hydro-training/output/lesson5/run_parameter_mods/DOMAIN/soil_properties.nc \
| tail -n 10
  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3 ;
}

The default value for refkdt is a global 3.0. We modify the parameter values using ncap2 to a lower value (0.1) to encourage less infiltration and more flow available for exfiltration to the channel.

In [29]:
%%bash
ncap2 -O -s "refkdt=refkdt*0.0+0.1" \
~/wrf-hydro-training/output/lesson5/run_parameter_mods/DOMAIN/soil_properties.nc \
~/wrf-hydro-training/output/lesson5/run_parameter_mods/DOMAIN/soil_properties.nc

We check to make sure the changes are as expected.

In [30]:
%%bash
ncdump -v refkdt ~/wrf-hydro-training/output/lesson5/run_parameter_mods/DOMAIN/soil_properties.nc \
| tail -n 10
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1 ;
}

Step 3b: Run the simulation

Now we are ready to run a simulation with our new modified parameters.

In [31]:
%%bash
cd ~/wrf-hydro-training/output/lesson5/run_parameter_mods
mpirun -np 2 ./wrf_hydro.exe >> run.log 2>&1

Check to make sure your run completed successfully.

In [32]:
%%bash
tail -1 ~/wrf-hydro-training/output/lesson5/run_parameter_mods/diag_hydro.00000
 The model finished successfully.......

Results

We will now look at the differences in streamflow between our baseline run with default parmeters and the simulation using our new modified parameters.

We will use Python and the xarray library to load the data and plot hydrographs. For an intro to these tools, please see Lesson 3.

In [33]:
chanobs_control = xr.open_mfdataset('/home/docker/wrf-hydro-training/output/lesson5/run_parameter_baseline/*CHANOBS*',
                            combine='by_coords')
chanobs_param_mods = xr.open_mfdataset('/home/docker/wrf-hydro-training/output/lesson5/run_parameter_mods/*CHANOBS*',
                            combine='by_coords')
In [34]:
fig, axes = plt.subplots(ncols=1,figsize=(12, 6))
plt.suptitle('Hydrographs for parameter experiment',fontsize=24)
chanobs_control.sel(feature_id = 3).streamflow.plot(label='Control, refkdt=3',
                                                        color='black',
                                                        linestyle='--')
chanobs_param_mods.sel(feature_id = 3).streamflow.plot(label='refkdt=0.1',
                                                        color='blue',
                                                        linestyle='-')
obs[obs['Location'] == 'Desague'].plot(x='dateTime',
                                       y='streamflow_cms',
                                       ax=axes,
                                       label='Observed',
                                       color='grey')
plt.ylim(0,15)
plt.legend()
plt.show()

Step 2c: Use the NCO command ncap2 to modify the bexp parameter values in the soil_properties.nc file.

BEXP

The bexp parameter controls how actual conductivity is scaled from saturated conductivity based on soil water content. Initial values are estimated based on soil texture class, but reported ranges have large (many orders of magnitude) variability.

The equation for the application of the BEXP parameter in the calculation of unsaturated hydaulic conductivity (WCND) the NoahMP code is as follows:

 WCND  = parameters%DKSAT(ISOIL) * FACTR2 ** EXPON
 where: 
   parameters%DKSAT(ISOIL) is the saturated hydraulic conductivity (DKSAT) value
   FACTR2 is fraction of actual to saturated soil water content [MAX(0.01, SMC/parameters%SMCMAX(ISOIL))]
   EXPON is 2.0*parameters%BEXP(ISOIL) + 3.0

First we make a new directory and repeat our previous change to refkdt so that we are building on the previous experiment.

In [85]:
%%bash
cp -r ~/wrf-hydro-training/output/lesson5/run_default_template \
~/wrf-hydro-training/output/lesson5/run_parameter_mods2

ncap2 -O -s "refkdt=refkdt*0.0+0.1" \
~/wrf-hydro-training/output/lesson5/run_parameter_mods2/DOMAIN/soil_properties.nc \
~/wrf-hydro-training/output/lesson5/run_parameter_mods2/DOMAIN/soil_properties.nc

First, we check the current parameter values using ncdump. Note that these also vary spatially by soil type.

In [90]:
%%bash
ncdump -v bexp ~/wrf-hydro-training/output/lesson5/run_parameter_mods2/DOMAIN/soil_properties.nc \
| tail -n 10
    5.25,
  11.55, 11.55, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 
    5.25,
  11.55, 11.55, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 
    5.25,
  11.55, 11.55, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 
    5.25,
  5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25,
  5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25 ;
}

Here we will use the NCO command ncap2 to decrease by a factor of bexp by a factor of 0.5 to promote more runoff.

In [38]:
%%bash
ncdump -v bexp ~/wrf-hydro-training/output/lesson5/run_parameter_mods2/DOMAIN/soil_properties.nc \
| tail -n 10
    5.25,
  11.55, 11.55, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 
    5.25,
  11.55, 11.55, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 
    5.25,
  11.55, 11.55, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 
    5.25,
  5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25,
  5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25, 5.25 ;
}
In [91]:
%%bash
ncap2 -O -s "bexp=bexp*0.5" \
~/wrf-hydro-training/output/lesson5/run_parameter_mods2/DOMAIN/soil_properties.nc \
~/wrf-hydro-training/output/lesson5/run_parameter_mods2/DOMAIN/soil_properties.nc
ncdump -v bexp ~/wrf-hydro-training/output/lesson5/run_parameter_mods2/DOMAIN/soil_properties.nc \
| tail -n 10
    2.625, 2.625, 2.625,
  5.775, 5.775, 2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 
    2.625, 2.625, 2.625,
  5.775, 5.775, 2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 
    2.625, 2.625, 2.625,
  2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 
    2.625, 2.625, 2.625,
  2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 2.625, 
    2.625, 2.625, 2.625 ;
}
In [92]:
%%bash
cd ~/wrf-hydro-training/output/lesson5/run_parameter_mods2
mpirun -np 2 ./wrf_hydro.exe >> run.log 2>&1
In [93]:
chanobs_param_mods2 = xr.open_mfdataset('/home/docker/wrf-hydro-training/output/lesson5/run_parameter_mods2/*CHANOBS*',
                            combine='by_coords')
In [94]:
fig, axes = plt.subplots(ncols=1,figsize=(12, 6))
plt.suptitle('Hydrographs for parameter experiment',fontsize=24)
chanobs_control.sel(feature_id = 3).streamflow.plot(label='Control, w/refkdt=3',
                                                        color='black',
                                                        linestyle='--')
chanobs_param_mods.sel(feature_id = 3).streamflow.plot(label='refkdt=0.1',
                                                        color='blue',
                                                        linestyle='-')
chanobs_param_mods2.sel(feature_id = 3).streamflow.plot(label='refkdt=0.1, bexp*0.5',
                                                        color='green',
                                                        linestyle='-')
obs[obs['Location'] == 'Desague'].plot(x='dateTime',
                                       y='streamflow_cms',
                                       ax=axes,
                                       label='Observed',
                                       color='grey')
plt.ylim(0,20)
plt.legend()
plt.show()

Step 2d: Use the NCO command ncap2 to modify the dksat parameter values in the soil_properties.nc file.

dksat

As with most physically based hydrological models, the soil saturated hydraulic conductivity, dksat, controls the speed at which water moves through the subsurface. This is a sensitive parameter in the model, and while easy to measure at the point scale, dksat is tricky to estimate at the scale of kilometers. This is a common calibration parameter.
First we make a new directory and modify refkdt and bexp to the same values as the previous experiments so that we are building on our previous simulations.

In [116]:
%%bash
cp -r ~/wrf-hydro-training/output/lesson5/run_default_template \
~/wrf-hydro-training/output/lesson5/run_parameter_mods3

ncap2 -O -s "refkdt=refkdt*0.0+0.1" \
~/wrf-hydro-training/output/lesson5/run_parameter_mods3/DOMAIN/soil_properties.nc \
~/wrf-hydro-training/output/lesson5/run_parameter_mods3/DOMAIN/soil_properties.nc

ncap2 -O -s "bexp=bexp*0.5" \
~/wrf-hydro-training/output/lesson5/run_parameter_mods3/DOMAIN/soil_properties.nc \
~/wrf-hydro-training/output/lesson5/run_parameter_mods3/DOMAIN/soil_properties.nc

Here we will view the default dksat and then reduce dksat by a factor of 0.5 using the ncap2 NCO command.

In [117]:
%%bash
ncdump -v dksat ~/wrf-hydro-training/output/lesson5/run_parameter_mods3/DOMAIN/soil_properties.nc \
| tail -n 10
    3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06,
  9.74e-07, 9.74e-07, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 
    3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06,
  9.74e-07, 9.74e-07, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 
    3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06,
  3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 
    3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06,
  3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 
    3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06, 3.38e-06 ;
}
In [118]:
%%bash
ncap2 -O -s "dksat=dksat*0.5" \
~/wrf-hydro-training/output/lesson5/run_parameter_mods3/DOMAIN/soil_properties.nc \
~/wrf-hydro-training/output/lesson5/run_parameter_mods3/DOMAIN/soil_properties.nc
ncdump -v dksat ~/wrf-hydro-training/output/lesson5/run_parameter_mods3/DOMAIN/soil_properties.nc \
| tail -n 10
    1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06,
  4.87e-07, 4.87e-07, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 
    1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06,
  4.87e-07, 4.87e-07, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 
    1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06,
  1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 
    1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06,
  1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 
    1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06, 1.69e-06 ;
}

Step 2e: Use the NCO command ncap2 to modify the slope parameter values in the soil_properties.nc file.

SLOPE

Another important Noah/NoahMP parameter that we commonly adjust is slope. Originally estimated based on land surface topography (hence the name slope), the slope parameter actually controls how open or closed the bottom boundary of the soil column is. Values range from 0 to 1, where 0 is a completely closed bottom boundary and 1 is completely open. Lower slope values will keep more water in the soil column, while higher values will allow more water to drain to the channel or to deeper baseflow stores, depending on the selected baseflow physics options.

As with refkdt and dksat above, we will use the NCO command ncap2 to modify the slope value.

First, we check the current parameter values using ncdump.

In [119]:
%%bash
ncdump -v slope ~/wrf-hydro-training/output/lesson5/run_parameter_mods3/DOMAIN/soil_properties.nc \
| tail -n 10
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
  0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1 ;
}

The default value for slope is a global 0.1. We modify the parameter values using ncap2 to open the bottom boundary (slope*0+0.5) and increase the deeper baseflow component.

In [120]:
%%bash
ncap2 -O -s "slope=slope*0+0.5" \
~/wrf-hydro-training/output/lesson5/run_parameter_mods3/DOMAIN/soil_properties.nc \
~/wrf-hydro-training/output/lesson5/run_parameter_mods3/DOMAIN/soil_properties.nc
ncdump -v slope ~/wrf-hydro-training/output/lesson5/run_parameter_mods3/DOMAIN/soil_properties.nc \
| tail -n 10
  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5 ;
}

Note: we are running this in the same dir as before, but will be able to plot each of the experiments because the datasets are loaded.

In [121]:
%%bash
cd ~/wrf-hydro-training/output/lesson5/run_parameter_mods3
mpirun -np 2 ./wrf_hydro.exe >> run.log 2>&1
In [122]:
chanobs_param_mods3 = xr.open_mfdataset('/home/docker/wrf-hydro-training/output/lesson5/run_parameter_mods3/*CHANOBS*',
                            combine='by_coords')
In [123]:
fig, axes = plt.subplots(ncols=1,figsize=(12, 6))
plt.suptitle('Hydrographs for parameter experiment',fontsize=24)
chanobs_control.sel(feature_id = 3).streamflow.plot(label='Control, w/refkdt=3, slope=0.1',
                                                        color='black',
                                                        linestyle='--')
chanobs_param_mods.sel(feature_id = 3).streamflow.plot(label='refkdt=0.1, slope=0.1',
                                                        color='blue',
                                                        linestyle='-')
chanobs_param_mods2.sel(feature_id = 3).streamflow.plot(label='refkdt=0.1, slope=0.1, bexp*0.5',
                                                        color='green',
                                                        linestyle='-')
chanobs_param_mods3.sel(feature_id = 3).streamflow.plot(label='refkdt=0.1, bexp*0.5, slope=0.5, dksat*0.5',
                                                        color='orange',
                                                        linestyle='-')

obs[obs['Location'] == 'Desague'].plot(x='dateTime',
                                       y='streamflow_cms',
                                       ax=axes,
                                       label='Observed',
                                       color='grey')
plt.ylim(0,15)
plt.legend()
plt.show()

OPTIONAL::: RETDEPRTFAC

The RETDEPRTFAC parameter is a multiplier on the maximum retention depth. Ponded water on the surface above this retention depth threshold can be moved around the landscape via overland flow. The default value in the code is quite small (~0.001mm, though variable by terrain slope) to allow almost all ponded water to be available for routing. However, in many regions landscape features like wetlands, small detention ponds, and heavy vegetation litter/debris can trap water on the land surface. Increasing the RETDEPRTFAC multiplier will hold more ponded water on the surface before it becomes runoff.

We will use the NCO command ncap2 to reduce RETDEPRTFAC value to 0.

First, we check the current parameter values.

In [141]:
%%bash
cp -r ~/wrf-hydro-training/output/lesson5/run_default_template \
~/wrf-hydro-training/output/lesson5/run_parameter_mods4

ncap2 -O -s "refkdt=refkdt*0.0+0.1" \
~/wrf-hydro-training/output/lesson5/run_parameter_mods4/DOMAIN/soil_properties.nc \
~/wrf-hydro-training/output/lesson5/run_parameter_mods4/DOMAIN/soil_properties.nc

ncap2 -O -s "bexp=bexp*0.5" \
~/wrf-hydro-training/output/lesson5/run_parameter_mods4/DOMAIN/soil_properties.nc \
~/wrf-hydro-training/output/lesson5/run_parameter_mods4/DOMAIN/soil_properties.nc

ncap2 -O -s "dksat=dksat*0.5" \
~/wrf-hydro-training/output/lesson5/run_parameter_mods4/DOMAIN/soil_properties.nc \
~/wrf-hydro-training/output/lesson5/run_parameter_mods4/DOMAIN/soil_properties.nc

ncap2 -O -s "slope=slope*0+0.5" \
~/wrf-hydro-training/output/lesson5/run_parameter_mods4/DOMAIN/soil_properties.nc \
~/wrf-hydro-training/output/lesson5/run_parameter_mods4/DOMAIN/soil_properties.nc
In [142]:
%%bash
ncdump -v RETDEPRTFAC ~/wrf-hydro-training/output/lesson5/run_parameter_mods4/DOMAIN/Fulldom_hires.nc \
| tail -n 10
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1,
  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
    1, 1, 1, 1, 1, 1, 1, 1, 1 ;
}

Then, we modify the parameter values using ncap2 and confirm our changes.

In [143]:
%%bash
ncap2 -O -s "RETDEPRTFAC=RETDEPRTFAC*10000.0" \
~/wrf-hydro-training/output/lesson5/run_parameter_mods4/DOMAIN/Fulldom_hires.nc \
~/wrf-hydro-training/output/lesson5/run_parameter_mods4/DOMAIN/Fulldom_hires.nc
ncdump -v RETDEPRTFAC ~/wrf-hydro-training/output/lesson5/run_parameter_mods4/DOMAIN/Fulldom_hires.nc \
| tail -n 10
    10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 
    10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 
    10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 
    10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 
    10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 
    10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 
    10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 
    10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 
    10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000 ;
}

Step 3b: Run the simulation

Now we are ready to run a simulation with our new modified parameters.

In [144]:
%%bash
cd ~/wrf-hydro-training/output/lesson5/run_parameter_mods4
mpirun -np 2 ./wrf_hydro.exe >> run.log 2>&1

Check to make sure your run completed successfully.

In [145]:
%%bash
tail -1 ~/wrf-hydro-training/output/lesson5/run_parameter_mods4/diag_hydro.00000
 The model finished successfully.......

Results

We will now look at the differences in streamflow between our baseline run with default parmeters and the simulation using our new modified parameters.

We will use Python and the xarray library to load the data and plot hydrographs. For an intro to these tools, please see Lesson 3.

Load the CHANOBS streamflow datasets and plot the results.

We are going to use the CHANOBS files because it will limit the number of grid cells to only those which we have specified have a gage.

In [146]:
chanobs_param_mods4 = xr.open_mfdataset('/home/docker/wrf-hydro-training/output/lesson5/run_parameter_mods4/*CHANOBS*',
                            combine='by_coords')
In [147]:
fig, axes = plt.subplots(ncols=1,figsize=(12, 6))
plt.suptitle('Hydrographs for parameter experiment',fontsize=24)
chanobs_control.sel(feature_id = 3).streamflow.plot(label='Control, w/refkdt=3, slope=0.1',
                                                        color='black',
                                                        linestyle='--')
chanobs_param_mods.sel(feature_id = 3).streamflow.plot(label='refkdt=0.1, slope=0.1',
                                                        color='blue',
                                                        linestyle='-')
chanobs_param_mods2.sel(feature_id = 3).streamflow.plot(label='refkdt=0.1, slope=0.1, bexp*2',
                                                        color='green',
                                                        linestyle='-')
chanobs_param_mods3.sel(feature_id = 3).streamflow.plot(label='refkdt=0.1, slope=0.1, bexp*0.5, dksat*0.5',
                                                        color='orange',
                                                        linestyle='-')
chanobs_param_mods4.sel(feature_id = 3).streamflow.plot(label='refkdt=0.1, bexp*0.5, dksat*0.5, slope=0.5, retdeprtfac*10000',
                                                        color='red',
                                                        linestyle='-')


obs[obs['Location'] == 'Desague'].plot(x='dateTime',
                                       y='streamflow_cms',
                                       ax=axes,
                                       label='Observed',
                                       color='grey')
plt.ylim(0,20)
plt.legend()
plt.show()

Next up - Do it yourself!

This concludes Lesson 5a. Spend some time creating your own parameter and physics experiments.

IT IS BEST TO EITHER SHUTDOWN THIS LESSON OR CLOSE IT BEFORE PROCEEDING TO THE NEXT LESSON TO AVOID POSSIBLY EXCEEDING ALLOCATED MEMORY. Shutdown the lesson be either closing the browser tab for the lesson or selecting KERNAL->SHUTDOWN in the jupyter notebook toolbar.