3.8. Dynamic stochastic modelling: infiltration model¶

3.8.1. Introduction¶

Important note: questions in this section are NOT available in Blackboard so please answer in a text document that you hand in afterwards. This section is also not part of the Land Surface Process Modelling course.

3.8.2. Model structure and field data¶

In the exercises, you will use a rainfall-runoff model simulating the Hortonian overland flow during a single rainstorm event. This section provides a short introduction to the model. Note that you need at least 2 GB diskspace for the exercises.

You will use data from an approximately rectangular 7500 m2 sandy loam hillslope with a vineyard, located in the Ouveze river basin, S. France. The model simulates rainfall, infiltration and overland flow. Interception is ignored here for reasons of simplicity (note that it is rather small anyway on a vineyard), as is surface storage.

Open the model runoff.py and have a look at the code. Run the model.

Display the input maps

aguila dem.map + ldd.map out.map + ldd.map &

The map dem.map is the digital elevation model of the modelling area, ldd.map is the local drain direction map with the drainage directions. out.map gives the ouflow location of the catchment.

Rainfall is modelled as a block rainstorm (constant rainfall at the start of the modelling period abruptly finishing half way):

aguila --timesteps=[1,360,1] --scenarios={1} p

Right-click on the legend for options and to open a timeseries (click on map to get the values for a location).

It is assumed that for each cell and each time step, infiltration capacity (m/h) equals the saturated conductivity ( $Z(\bf{s})$ , m/h). For each cell, the actual infiltration ( $I(\bf{s})$ , m/h) is the minimum value of the saturated conductivity and the water available for infiltration (m/h). The water available for infiltration is rainfall ( $P(\bf{s})$ , m/h) plus runon from upstream cells ( $R(\bf{s})$ , m/h). In an equation:

$I(\bf{s}) = min(K(\bf{s}),P(\bf{s}) + R(\bf{s}))$

where min(a,b) assigns the minimum value of a and b.

Runoff is routed using the Manning equation and the kinematic wave in downstream direction over the local drain direction map (ldd.map).

Saturated conductivity was measured at 10 locations on the hillslope using ring infiltrometer experiments. The average saturated conductivity was 0.05 m/h (n=10), which is quite normal for this type of soil. The standard deviation of the saturated conductivity data set was very high. In addition, discharge was measured at the outflow point of the catchment. With rainstorms having similar characteristics as the one used here, peak discharge was mostly about 0.02 m3/s. In the following exercises, you will try to simulate the discharge from the catchment, using an average K value of 0.02 m/h. The question is whether simulated discharge is in the same order of magnitude as measured discharge at the outflow point.

3.8.3. Deterministic modelling using mean K¶

The run with the model you did in the previous section was done using a saturated conductivity value of 0.05 m/h for all cells. This value is the same as the average value derived from the ring infiltrometre measurements. The question is how much runoff was generated by this run.

Let’s first have a look at the change through time of actual infiltration (m/h, the files i0000000.001, representing the first time step, up to i0000000.360, representing the last time step), runoff (m3/h, q0000000.001 up to q0000000.360), and actual infiltration as a percentage of the saturated conductivity (%, ip000000.001, ..). Note that these are written to the folder 1 but you do not need to care about this as it is indicated by scenarios in the aguila commane. Type from the folder where you were running the script:

aguila --timesteps=[1,360,1] --scenarios={1} i ip q + ldd.map &

and animate the maps. By right clicking on the legend you can open a time series view for a location on the map (just click on the map to change the location for which you want the time series). Note that the time step duration is 10 seconds. Click on the location of the outflow point to show the time series for that location.

Question

What is the peak discharge at the outflow point?

Question

Explain why such a small amount of discharge is generated (compared to the measured peak discharges for these kind of rainstorms).

3.8.4. The stochastic model¶

The deterministic model used in the previous exercise does not generate runoff since it ignores spatial variation of saturated conductivity on the field. To give a more realistic representation of runoff, it is needed to include this spatial variation in our model. The problem is that there is insufficient data to represent this in a deterministic way, using a map with the actually occurring pattern of saturated conductivity on the field. The only way to solve the problem is to represent spatial variation of saturated conductivity as a stochastic variable, having a defined spatial probability distribution.

We assume a lognormal distribution of saturated conductivity by defining:

$K(\bf{s}) = e^{Z(\bf{s})}$

It is assumed that $Z(\bf{s})$ is a multivariable normal and stationary random spatial function. Thus, $Z(\bf{s})$ is defined by its expectation $m_{Z(\bf{s})}$ , its variance $\sigma_{Z(\bf{s})}$ and its spatial correlation structure. The spatial correlation structure could be defined by a variogram, but we will use a simplified approach here to create random fields very similar to thos created with a variogram.

The expectation (mean) of $K(\bf{s})$ is:

$m_{K(\bf{s})} = e^{m_{Z(\bf{s})} + 0.5 \sigma_{Z(\bf{s})}}$

Increasing $\sigma_{Z(\bf{s})}$ results in a probability distribution of $K(\bf{s})$ with both a higher variance and a higher skewness.

The model can be run in stochastic mode, using user defined spatial probability distributions of $K(\bf{s})$ . By modifying the following components of the script runoff.py, you can create scenario’s with different spatial probability distributions of $K(\bf{s})$ . In runoff.py you can define:

The number of Monte Carlo loops. The number of Monte Carlo loops is set in the following line at the upper part of the script:

nrSamples = 2

For instance, to run 50 loops, change the line to:

nrSamples = 50

Somewhat down, the value of $m_{K(\bf{s})}$ :

# mean of Z
mean = -2.9957

The value of $\sigma_{Z(\bf{s})}$ :

# variance of Z
var = 0.000001

The spatial correlation structure of $Z(\bf{s})$ . The random field is generated by drawing independent realizations for each pixel. This results in a random field corresponding to a variogram without range (only nugget). Spatially correlated versions are created by applying a smoothing window. This creates realizations with a spatial correlation structure. The larger the window, the larger the scale of variation (larger variogram range). This is set by the value of approximateRangeInPixels. A value of one corresponds to no spatial correlation (nugget only variogram), other possible values are 3, 5, 7, etc. It is now set to 1 (no spatial correlation):

# approximate range of random field (pixels), use 1, 3, 5, 7, ..
# note that 1 implies 'white noise', i.e. no spatial correlation
approximateRangeInPixels = 1

Using the information given above, make a copy of runoff.py and save it as runoff_stoch.py to run the model with:

50 Monte Carlo loops
white noise, no spatial correlation
a value of $\sigma_{Z(\bf{s})}$ of 1.0
an expectation (mean) of $K(\bf{s})$ of $m_{K(\bf{s})} = 0.05$ m/h

Question:

With the input values given above, what value needs to be used for $m_{Z(\bf{s})}$ ?

Write down all changes you have made to create runoff_stoch.py. You will need these in the next exercises!

Run the script runoff_stoch.py.

After running the script, display realizations of $K(\bf{s})$ :

aguila --scenarios='{1,2,3,4,5,6,7,8,9,10,11,12}' --multi=3x4 ks

Right-click on the legend for options.

You can also display the cumulative probability distributions as they are calculated at the bottom of the script from the realisations:

aguila --quantiles=[0.05,0.95,0.05] ks

And the mean value:

aguila ks-ave.map

Question

Based on visual interpretation of the values and patterns on the maps of $K(\bf{s})$ , do you think these agree with the parameters of the spatial probability distribution of $K(\bf{s})$ which you defined in the script?

3.8.5. Evaluating individual realizations of the model¶

For each Monte Carlo loop, the rainfall-runoff model is run with the realization of $K(\bf{s})$ for that loop, resulting in a realization of all outputs of the model. To minimize the runtime of the model, only limited number of output variables is stored on the harddisk.

Animate the actual infiltration i, runoff q, actual infiltration as a percentage of K ip, and the saturated conductivity ks, for 12 Monte Carlo loops, type:

aguila --timesteps=[1,360,1] --scenarios='{1,2,3,4,5,6,7,8,9,10,11,12}' --multi=3x4 i + ldd.map  q + ldd.map ip + ldd.map ks + ldd.map

Right click on the legend for options, e.g. to open a timeseries plot.

Question

While the deterministic model did not generate runoff, it is clear that most realizations of the stochastic model do generate runoff (using the same mean saturated conductivity). Explain this difference.

Question

In an animation, stop at time step 100. Explain in detail the relation between infiltration, runoff, and saturated conductivity for this timestep. Include the spatial patterns in your answer.

3.8.6. Evaluating the stochastic output¶

Since interpretation of individual realizations is difficult, it is better to calculate sampling statistics of the realizations. In a dynamic (temporal) model, this means that sample statistics are calculated for each time step (over all Monte Carlo loops). Files with these sample statistics are stored in the main directory (where your Python script is stored). The script calculates these statistics only for eacht 10-th time step to reduce runtime and storage.

Plot cumulative probability distributions of the discharge q (m3/s) and infiltration i (m/h):

aguila --timesteps=[10,360,10] --quantiles=[0.05,0.95,0.05] q + ldd.map i + ldd.map

Use the options explained previously in stochastic modelling to plot the timeseries of the 25 percentile at the outflow point. You will need to get the probability plot (right click on the menu).

Question

Explain how the 25 percentile of the hydrograph is calculated.

Question

Somebody has bought quite a small flume for discharge measurements to be installed at the outflow point of the catchment. The manual of the flume says it cannot measure discharge above 0.009 m3/s. During which period(s) of the rainstorm is it certain (using a 50% percent confidence interval) that the discharge at the outflow point is above 0.009 m3/s?

Question

Plot the median discharge curve (50% percentile) at the outflow point. Write down the modelled median peak flow. Compare it with the measured peak flow of 0.05 m3/s. Do you think our model is ‘good’?

Question

Write down the boundaries (discharge, m3/h) of the 50% confidence interval of the peak discharge and calculate the width of the confidence interval.

3.8.7. The effect of the distribution of saturated conductivity on discharge¶

It is interesting to study the relation between 1) the parameters describing the spatial probability distribution of $K(\bf{s})$ and 2) the hydrograph. Both the variance of $K(\bf{s})$ as well as the scale of the variation in $K(\bf{s})$ has an effect on the hydrograph.

Let’s first look at the effect of the variance in saturated conductivity. Copy runoff_stoch.py and save as runoff_stoch_low_var.py. In runoff_stoch_low_var.py, make changes such that the value of $\sigma_{Z(\bf{s})}$ becomes 0.5 (it was 1.0). Keep the expectation (mean) of $K(\bf{s})$ at $m_{K(\bf{s})} = 0.05$ m/h. Note that you need to change both the values of var and mean in the script. Run the script.

Question

What is the effect of decreasing the variance of saturated conductivity on 1) the median and 2) the 50% confidence interval of the peak discharge?

Now let’s look at the effect of the spatial pattern in saturated conductivity. Copy runoff_stoch.py and save it as runoff_stoch_large_range.py. Change the value of approximateRangeInPixels to a value of 3.

After running the script, display realizations of $K(\bf{s})$ :

aguila --scenarios='{1,2,3,4,5,6,7,8,9,10,11,12}' --multi=3x4 ks

You will notice that the spatial pattern has changed to a spatially correlated pattern, with a larger patch size.

Question

What is the effect of increasing the spatial range of saturated conductivity on 1) the median and 2) the 50% confidence interval of the peak discharge?

For a detailed explanation of the simulated relations between saturated conductivity and infiltration (and thus discharge), refer to https://doi.org/10.1016/j.advwatres.2005.06.012