NCSU GIS 714:
Geospatial Computing and Simulations

Surface water flow and erosion simulation

In this assignment we explore simulation of spatial pattern of overland flow depth by running the SIMWE model which employs path sampling method to solve the shallow water flow equations. We also compute sediment transport and net erosion and deposition using the path sampling model.

Resources:

For animating in GRASS GIS 7, see the video instructions or instructions for the Spatio-temporal data handling and visualization in GRASS GIS workshop

Start GRASS GIS

Start GRASS - click on GRASS icon or type
grass

In startup pannel set GRASS GIS Database Directory to path to datasets, for example on MS Windows, C:\Users\myname\grassdata.
For GRASS Location select nc_spm_08_grass7 (North Carolina, State Plane, meters) and for GRASS Mapset create a new mapset (called e.g. HW_interpolation_2).
Click Start GRASS session.

Change working directory:
Settings > GRASS working environment > Change working directory > select/create any directory
or type cd (stands for change directory) into the GUI Console and hit Enter: 

cd

Overland flow depth and discharge

Set the computational region to a small agricultural watershed and define the 2m resolution: 

g.region rural_1m res=2 -p

Calculate partial derivatives to define the gradient vector of elevation surface: 

v.surf.rst -d input=elev_lid792_bepts elevation=elev_lid792_2m slope=dx_2m aspect=dy_2m tension=15 smooth=1.5 npmin=150
Note that partial derivatives can be also computed using r.slope.aspect. Why do we need to compute elevation surface gradient vector for flow simulation? What does it control?

Compute spatial pattern of overland flow depth and discharge by running the SIMWE model implemented in the r.sim.water module - read the manual page to understand the parameters. We first use uniform rainfall excess, infiltration and land cover. 

r.sim.water -t elevation=elev_lid792_2m dx=dx_2m dy=dy_2m rain_value=50 infil_value=0 man_value=0.05 depth=wdp_2m discharge=disch_2m nwalkers=100000 niterations=30 output_step=2
Display the results: 
d.rast wdp_2m.30
d.rast disch_2m.30
Add legend to the final raster maps and save the images. Comment on the resulting pattern - where is the water ponding and what is causing it at this location?

Optionally, animate the time series: File > Animation tool

Peak runoff with preferential flow direction in the stream channel

Compute peak runoff with predefined flow direction along the stream inluding a culvert under the road.

First, compute direction (aspect) of the given streams: 

v.to.rast streams output=streams_dir_2m use=dir

Compute stream dx and dy using the stream direction angle and a slope angle equal to 2 degrees: 

r.mapcalc "dxstr_2m = tan(2.)*cos(streams_dir_2m)"
r.mapcalc "dystr_2m = tan(2.)*sin(streams_dir_2m)"

Compute flow gradient vector by combining dx and dy derived from the DEM and stream: 

r.mapcalc "dxdemstr_2m = if(isnull(dxstr_2m), dx_2m, dxstr_2m)"
r.mapcalc "dydemstr_2m = if(isnull(dystr_2m), dy_2m, dystr_2m)"

Run the overland flow model: 

r.sim.water -t elevation=elev_lid792_2m dx=dxdemstr_2m dy=dydemstr_2m rain_value=50 infil_value=0 man_value=0.05 depth=wdpstr_2m discharge=dischstr_2m nwalkers=100000 niterations=30 output_step=2
d.rast dischstr_2m.30
Add legend to the final raster maps and save the images. Comment on the resulting pattern - is the water still ponding? If not, why is it now flowing under the road?

Optionally, use Animation tool to animate the discharge.

Runoff for spatially variable landcover and rainfall excess

Compute runoff for spatially variable landcover represented by spatially variable Mannings coefficient and rainfall excess (rainfall intensity - infiltration rate).

We will use variable Mannings coefficient defined by reclassifying land cover class. Find out the land use categories in the map landcover_1m: 

r.category landcover_1m
d.rast landcover_1m
d.legend landcover_1m
We assign each category a Mannings roughness coefficient (published in literature) by recoding the landcover map using the r.recode module and reclassification rules stored in this file land_to_mannings.txt and also shown here: 
1:1:0.9:0.9
2:2:0.5:0.5
3:3:0.01:0.01
4:4:0.03:0.03
5:5:0.01:0.01
6:6:0.05:0.05
7:7:0.1:0.1
8:8:0.1:0.1
9:9:0.9:0.9
10:10:0.03:0.03
11:11:0.5:0.5

r.recode input=landcover_1m output=mancover_2m rules=land_to_mannings.txt
Similarly, we will create raster with spatially variable rainfall excess rates based on the land cover classes. We use file named land_to_rain.txt to specify the rates for individual classes. 
1:1:50:50
2:2:5:5
3:3:40:40
4:4:35:35
5:5:50:50
6:6:40:40
7:7:25:25
8:8:15:15
9:9:50.:50.
10:10:40:40
11:11:10:10
Again, we use the file as rules for the r.recode module. 
r.recode input=landcover_1m output=raincover_2m rules=land_to_rain.txt

Run the model: 

r.sim.water -t elevation=elev_lid792_2m dx=dxdemstr_2m dy=dydemstr_2m rain=raincover_2m infil_value=0 man=mancover_2m depth=wdpstrcov_2m discharge=distrcov_2m nwalkers=100000 niterations=30 output_step=2
Display the static results and save the images for your paper. 
d.rast wdpstrcov_2m.30
d.legend wdpstrcov_2m.30
d.rast distrcov_2m.30
d.legend distrcov_2m.30
Comment on the resulting pattern - how is the flow depth related to the landcover? Which landcover holds most water?

Optionaly animate the time series using the Animation tool.

Sediment flow rate, erosion and deposition

Compute sediment flow rate and net erosion/deposition using sediment transport part of the SIMWE model implemented in the the r.sim.sediment module. To make the computations faster, set region just to the upper part of the watershed: 
g.region s=s+290

Compute input transport capacity and detachment coefficient maps: 

r.mapcalc "tranin = 0.001"
r.mapcalc "detin = 0.001"

Compute input critical shear stress: 

r.mapcalc "tauin = 0.01"

Run the model using the last depth from previous run: 

g.copy rast=wdp_2m.30,wdp_2m
r.sim.sediment elevation=elev_lid792_2m dx=dx_2m dy=dy_2m water_depth=wdp_2m detachment_coeff=detin transport_coeff=tranin shear_stress=tauin man_value=0.05 nwalkers=1000000 niterations=30 transport_capacity=tcapacity tlimit_erosion_deposition=erdepmax sediment_flux=sedflow erosion_deposition=erdepsimwe

Display these results after few seconds: 

d.rast tcapacity
d.rast erdepmax

Display the final results: 

d.rast sedflow
d.legend sedflow
d.rast erdepsimwe
d.legend erdepsimwe