Lesson 5 - Land surface experiments

Overview

In Lesson 4 we experimented with different initial and boundary conditions (precipitation).

In this lesson, we experiment with different terrain physics options as well as manipulate a few parameters 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.

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.

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

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.

Now run the baseline model simulation as-is.

Check that the simulation 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.

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

Check that the simulation 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.

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

Check that the simulation 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

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.

Plot the hydrographs

We see definite event differences between the streamflow responses with and without overland routing active (black dashed and blue solid lines). When overland routing is turned off, excess surface runoff from the land model is scraped off the cells and added to the stream directly, so we tend to see earlier and flashier responses to precipitation. The overland routing processes slowdown the runoff and allow water to re-infiltrate where soil storage space is available.

When subsurface routing is also turned off (red solid line), the earlier/flashier peaks are also present, but we lose the secondary pulse and have lower flow after the event. Where is that water going?

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.

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.

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

Load the GWOUT groundwater bucket model output files

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 the simplified physics configuration without lateral routing processes, infiltrated water is not permitted to flow laterally. It has to make its way vertically through the 2-meter soil column, then through the deep groundwater "buckets" before it becomes streamflow. These are often much slower flowpaths than lateral flowpaths. Therefore, this water is first wetting up the soil column and the bucket stores, then will eventually make its way to the channel as baseflow.

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:

In this lesson, we will manipulate parameters in the soil_properties.nc (refkdt, dksat, bexp, slope) and Fulldom_hires.nc (LKSATFAC, RETDEPRTFAC) files.

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.

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.

Launch the baseline run and make sure it completes 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.

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.

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

SLOPE

There are a number of Noah/NoahMP model parameters that affect water partitioning. One important 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.

Here we will use the NCO command ncap2 to increase the slope value.

First, we check the current parameter values using ncdump.

The default value for slope is a global 0.1. We modify the parameter values using ncap2 to open up the bottom boundary (slope=0.2) and increase the deeper baseflow component. We will discuss the baseflow model more in Lesson 5b.

We check to make sure the changes are as expected.

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

REFKDT

An 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 can become 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.

The default value for refkdt is a global 3.0. We modify the parameter values using ncap2 to a lower value (1.8).

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

DKSAT and BEXP

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. Initial values are estimated based on soil texture class, but reported ranges have large (many orders of magnitude) variability. This is a common calibration parameter, along with the related bexp parameter that controls how actual conductivity is scaled from saturated conductivity based on soil water content.

Here we will use the NCO command ncap2 to increase the dksat value by a factor of 3 and bexp by a factor of 2.

First, we check the current parameter values using ncdump.

Note that this parameter's initial values vary in space based on soil texture class. We modify the parameter values using ncap2 to increase by a factor of 3. This maintains the spatial distribution, but encourage faster soil water movement overall.

Use ncdump to check the starting values of bexp. Note that these also vary spatially by soil type. We will double these values.

Step 3: Modify terrain routing parameters using NCO tools

Now we will turn our attention to the parameters that affect lateral routing behavior.

Step 3a: Use the NCO command ncap2 to modify the LKSATFAC parameter values in the Fulldom_hires.nc file

LKSATFAC

The Fulldom_hires.nc file contains two parameters that are also important for lateral flow processes. The LKSATFAC parameter is a multiplier on the prescribed lateral saturated hydraulic conductivity values specified in hydro2dtbl.nc. By default, lateral conductivity in hydro2dtbl.nc matches vertical conductivity specified in soil_properties.nc. However, in the real world we frequently see many orders of magnitude higher conductivities in the lateral direction vs. the vertical direction (due to soil stratigraphy, preferential flowpaths caused by roots and animals, etc.). LKSATFAC is an easy way to adjust this anisotropy, and by default it is set to 1,000.

We will use the NCO command ncap2 to increase the LKSATFAC value to 2000.

First, we check the current parameter values.

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

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 increase the RETDEPRTFAC value to 10000 (max of 10mm).

First, we check the current parameter values.

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

Step 3b: Run the simulation

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

Check to make sure your run completed 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

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.

Plot the hydrographs

Discussion

Our baseline simulation with default parameters had an initial peak that was a bit too early and flashy and a muted secondary peak. We increased slope to better blend the hydrograph peaks and modified refkdt to better match our activated terrain physics. We increased dksat and bexp (controlling the vertical drainage rate) and LKSATFAC (controlling the lateral drainage rate) to add more spread to the flow peak. We increased RETDEPRTFAC to hold more ponded water at the surface. Our modified parameter experiment results do better match the gage observations, though we have a delayed peak which might be addressed through additional parameter manipulation (e.g., groundwater parameters).

IMPORTANT NOTE: We abbreviated this parameter calibration exercise to fit within a short lesson. In practice, you would want to give the model time to adjust to a new parameter set by running an extended "spin-up" period before the time period you are evaluating. You would also want to expose the model to a wider range of conditions than the short event demostrated here. Model parameters calibrated to a short event may not transfer well to other time periods. Good practice is to calibrate your model to a wide range of conditions to minimize the impacts of forcing, gage, or model physics/parameter uncertainties.

Next up - Do it yourself!

This concludes Lesson 5. 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 Kernel -> Shut Down Kernel in JupyterLab.

© UCAR 2020