The wflow_routing Model

Introduction

The wflow routing module uses the pcraster kinematic wave to route water over a DEM. By adding a bankfull level and a floodplainwidth to the configuration the model can also include estimated flow over a floodplain. In addition, simple reservoirs can be configured.

Method

A simple river profile is defined using a river width a bankfull heigth and a floodplainwidth. A schematic drawing is shown in the figure below.

_images/wflow_routing.png

First the maximum bankfull flow for each cell is determined using:

\[Q_b = (\frac{H_{b}}{\alpha{_{ch}}} * Bw)^{1/\beta}\]

Next the channel flow is determined by taking the mimumm value of the total flow and the maximum banfull flow and the floodplain flow is determined by subtracting the bankfull flow from the total flow.

In normal conditions (below bankfull), the waterlevel in the river is determined as follows:

\[H_{ch} = \alpha{_{ch}} {Q_{ch}}^{\beta}/ Bw\]

Where \(H_{ch}\) is the water level for the channel, \(\alpha{_{ch}}\) is the kinematic wave coefficient for the channel, \(Q_{ch}\) is the discharge in the channel and \(Bw\) is the width of the river.

If the water level is above bankfull the water level on the floodplains is calculated as follows:

\[H_{fp} = \alpha{_{fp}} {Q_{fp}}^{\beta}/ (Bw + P_{fp})\]

where \(H_{ch}\) is the water level on the floodplain, \(Q_{fp}\) is the discharge on the floodplain and \(P_{fp}\) is the wetted perimiter of the floodplain and \(\alpha{_{fp}}\) is the kinematic wave coefficient for the floodplain,

The wetted perimiter of the channel, \(P_{ch}\), is determined by:

\[P_{ch} = 2.0 H_{ch} + Bw\]

The wetted perimiter of the floodplain is defined as follows:

\[ \begin{align}\begin{aligned}N = max(0.0001,1.0/(1.0 + exp(-c * H_{fp}) - 0.5) * 2.0\\P_{fp} = N W_{fp}\end{aligned}\end{align} \]

This first equation defines the upper half of an S or sigmoid curve and will return values between 0.001 and 1.0. The c parameter defines the sharpness of the function, a high value of c will turn this into a step-wise function while a low value will make the function more smooth. The default value for c = 0.5. For example, with this default value a floodplain level of 1m will result in an N value of 0.25 and 2m will return 0.46. In the second equation this fraction is multiplied with the maximum floodplain width \(W_{fp}\).

The \(\alpha\) for the channel and floodplain are calculated as follows:

\[ \begin{align}\begin{aligned}\alpha_{ch} = (n_{ch}/\sqrt{slope})^\beta P_{ch}^{(2.0 / 3.0)\beta}\\\alpha_{fp} = (n_{fp}/\sqrt{slope})^\beta P_{fp}^{(2.0 / 3.0)\beta}\end{aligned}\end{align} \]

In which slope is the slope of the river bed and floodplain and \(n_{ch}\) and \(n_{fp}\) represent the manning’s n for the channel and floodplain respectively.

A compound \(\alpha_{total}\) is estimated by first calculating a compound n value \(n_{total}\):

\[ \begin{align}\begin{aligned}n_{total} = (P_{ch}/P_{total} n_{ch}^{3/2} + P_{fp}/P_{total} n_{fp}^{3/2})^{2/3}\\\alpha_{total} = (n_{total}/\sqrt{slope})^\beta (P_{fp} + P_{ch})^{(2.0 / 3.0)\beta}\end{aligned}\end{align} \]

The \(\alpha_{total}\) is used in the pcraster kinematic function to get the discharge for the next timestep.

The code is implemented in the updateRunoff attribute of the model class as follows:

self.Qbankfull = pow(self.bankFull/self.AlphaCh * self.Bw,1.0/self.Beta)
self.Qchannel = min(self.SurfaceRunoff,self.Qbankfull)
self.floodcells  = boolean(ifthenelse(self.WaterLevelCH > self.bankFull, boolean(1), boolean(0)))
self.Qfloodplain = max(0.0,self.SurfaceRunoff - self.Qbankfull)

self.WaterLevelCH = self.AlphaCh * pow(self.Qchannel, self.Beta) / (self.Bw)
self.WaterLevelFP = ifthenelse(self.River,self.AlphaFP * pow(self.Qfloodplain, self.Beta) / (self.Bw + self.Pfp),0.0)
self.WaterLevel = self.WaterLevelCH + self.WaterLevelFP

# Determine Qtot as a check
self.Qtot = pow(self.WaterLevelCH/self.AlphaCh * self.Bw,1.0/self.Beta) + pow(self.WaterLevelFP/self.AlphaFP * (self.Pfp + self.Bw),1.0/self.Beta)
# wetted perimeter (m)
self.Pch = self.wetPerimiterCH(self.WaterLevelCH,self.Bw)
self.Pfp = ifthenelse(self.River,self.wetPerimiterFP(self.WaterLevelFP,self.floodPlainWidth,sharpness=self.floodPlainDist),0.0)

# Alpha
self.WetPComb = self.Pch + self.Pfp
self.Ncombined = (self.Pch/self.WetPComb*self.N**1.5 + self.Pfp/self.WetPComb*self.NFloodPlain**1.5)**(2./3.)
self.AlpTermFP = pow((self.NFloodPlain / (sqrt(self.SlopeDCL))), self.Beta)
self.AlpTermComb = pow((self.Ncombined / (sqrt(self.SlopeDCL))), self.Beta)
self.AlphaFP = self.AlpTermFP * pow(self.Pfp, self.AlpPow)
self.AlphaCh = self.AlpTerm * pow(self.Pch, self.AlpPow)
self.Alpha = ifthenelse(self.River,self.AlpTermComb * pow(self.Pch + self.Pfp, self.AlpPow),self.AlphaCh)
self.OldKinWaveVolume = self.KinWaveVolume
self.KinWaveVolume = (self.WaterLevelCH * self.Bw * self.DCL) + (self.WaterLevelFP * (self.Pfp + self.Bw) * self.DCL)

Reservoirs

Simple reservoirs can be included within the kinematic wave routing by supplying a map with the locations of the reservoirs in which each reservoir has a unique id. Furthermore a set of lookuptables must be defined linking the reservoir-id’s to reservoir characteristics:

  • ResTargetFullFrac.tbl - Target fraction full (of max storage) for the reservoir: number between 0 and 1

  • ResTargetMinFrac.tbl - Target minimum full fraction (of max storage). Number between 01 and 1 < ResTargetFullFrac

  • ResMaxVolume.tbl - Maximum reservoir storage (above which water is spilled) [m^3]

  • ResDemand.tbl - Water demand on the reservoir (all combined) m^3/s

  • ResMaxRelease.tbl - Maximum Q that can be released if below spillway [m^3/s]

By default the reservoirs are not included in the model. To include them put the following lines in the .ini file of the model.

[modelparameters]
# Add this if you want to model reservoirs
ReserVoirLocs=staticmaps/wflow_reservoirlocs.map,staticmap,0.0,0
ResTargetFullFrac=intbl/ResTargetFullFrac.tbl,tbl,0.8,0,staticmaps/wflow_reservoirlocs.map
ResTargetMinFrac=intbl/ResTargetMinFrac.tbl,tbl,0.4,0,staticmaps/wflow_reservoirlocs.map
ResMaxVolume=intbl/ResMaxVolume.tbl,tbl,0.0,0,staticmaps/wflow_reservoirlocs.map
ResMaxRelease=intbl/ResMaxRelease.tbl,tbl,1.0,0,staticmaps/wflow_reservoirlocs.map
ResDemand=intbl/ResDemand.tbl,tblmonthlyclim,1.0,0,staticmaps/wflow_reservoirlocs.map

In the above example most values are fixed thought the year, only the demand is given per month of the year.

Configuration

The default name for the file is wflow_routing.ini.

Subcatchment flow

Normally the the kinematic wave is continuous throughout the model. By using the the SubCatchFlowOnly entry in the model section of the ini file all flow is at the subcatchment only and no flow is transferred from one subcatchment to another. This can be handy when connecting the result of the model to a water allocation model such as Ribasim.

Example:

[model]
SubCatchFlowOnly = 1

Forcing data

The model needs one set of forcing data: IW (water entering the model for each cell in mm). The name of the mapstack is can be defined in the ini file. By default it is inmaps/IW

See below for an example:

[inputmapstacks]
Inwater = /run_default/outmaps/IW
Inflow = /inmaps/IF

[run]
starttime = 1995-01-31 00:00:00
endtime = 1995-02-28 00:00:00
timestepsecs = 86400
reinit = 0

[model]
modeltype = routing
AnnualDischarge = 2290
Alpha = 120
WIMaxScale = 0.8
Tslice = 1
UpdMaxDist = 300000.0
reinit = 1
fewsrun = 0
OverWriteInit = 0
updating = 0
updateFile = no_set
sCatch = 0
intbl = intbl
timestepsecs = 86400
MaxUpdMult = 1.3
MinUpdMult = 0.7
UpFrac = 0.8
SubCatchFlowOnly = 0
wflow_subcatch = staticmaps/wflow_subcatch.map
wflow_dem = staticmaps/wflow_dem.map
wflow_ldd = staticmaps/wflow_ldd.map
wflow_river = staticmaps/wflow_river.map
wflow_riverlength = staticmaps/wflow_riverlength.map
wflow_riverlength_fact = staticmaps/wflow_riverlength_fact.map
wflow_gauges = staticmaps/wflow_gauges.map
wflow_inflow = staticmaps/wflow_inflow.map
wflow_riverwidth = staticmaps/wflow_riverwidth.map
wflow_floodplainwidth = staticmaps/wflow_floodplainwidth.map
wflow_bankfulldepth = staticmaps/wflow_bankfulldepth.map
wflow_landuse = staticmaps/wflow_landuse.map
wflow_soil = staticmaps/wflow_soil.map

[framework]
outputformat = 1
debug = 0
netcdfinput = None
netcdfoutput = None
netcdfstaticoutput = None
netcdfstaticinput = None
EPSG = EPSG:4326

[layout]
sizeinmetres = 0

[outputmaps]
self.SurfaceRunoff = _run
self.Qfloodplain = _qfp
self.Qchannel = _qch
self.Qbankfull = _qbnk
self.WaterLevelFP = _levfp
self.WaterLevelCH = _levch
self.InwaterMM = _IW
self.floodcells = fcel
self.Qtot = QQQ
self.Pch = ch
self.Pfp = fp
self.Alpha = al
self.AlphaCh = alch
self.AlphaFP = alfp
self.Ncombined = nc
self.MassBalKinWave = wat

[outputcsv_0]
samplemap = None

[outputtss_0]
samplemap = None

A description of the implementation of the kinematic wave is given on the pcraster website at http://pcraster.geo.uu.nl/pcraster/4.0.2/doc/manual/op_kinematic.html

In addition to the settings in the ini file you need to give the model additional maps or lookuptables in the staticmaps or intbl directories:

Lookup tables

N.tbl

Manning’s N for all no-river cells. Defaults to 0.072

N_River.tbl

Manning’s N for the river cells. Defaults to 0.036

N_FloodPlain.tbl

Manning’s N for the floodplain. A floodplain is always linked to a river cell. Defaults to 2* N of the river

ResDemand.tbl

Lookup table of demand (\(m^3/s\)) for reservoir locations.

ResMaxRelease.tbl

Lookup table of maximum release capcity (\(m^3/s\)) for reservoir locations.

ResMaxVolume.tbl

Lookup table of reservoi maximum volume (\(m^3\)) for reservoir locations.

ResTargetFullFrac

Lookup table of the target maximum full fraction (0-1) (\(m^3/s\)) for reservoir locations.

ResTargetMinFrac

Lookup table of the target minimum full fraction (0-1) (\(m^3/s\)) for reservoir locations.

As with all models the lookup tables can be replaced by a map with the same name (but with the .map extension) in the staticmaps directory.

staticmaps

wflow_subcatch.map

Map of the subcatchment in the area. Usually shared with the hydrological model

wflow_dem.map

The digital elevation model. Usually shared with the hydrological model

wflow_ldd.map

The D8 local drainage network.

wflow_river.map

Definition of the river cells.

wflow_riverlength.map

Optional map that defines the actual legth of the river in each cell.

wflow_riverlength_fact.map

Optional map that defines a multiplication factor for the river length in each cell.

wflow_gauges.map

Map of river gauges that can be used in outputting timeseries

wflow_inflow.map

Optional map of inflow points into the surface water. Limited testing.

wflow_riverwidth.map

Optional map of the width of the river for each river cell.

wflow_floodplainwidth.map

Optional map of the width of the floodplain for each river cell.

wflow_bankfulldepth.map

Optional map of the level at which the river starts to flood and water will also be conducted over the floodplain.

wflow_floodplaindist.map

Optional map that defines the the relation between bankfulldepth and the floodplaindepth. Default = 0.5

wflow_landuse.map

Required map of landuse/land cover. This map is used in the lookup tables to relate parameters to landuse/landcover. Usually shared with the hydrological model

wflow_soil.map

Required map of soil type. Usually shared with the hydrological model

initial conditions

The model needs the following files with initial conditions:

WaterLevelCH

Water level in the channel or the grid-cell water level for non-river cells.

WaterLevelFP

Water level in the floodplain

SurfaceRunoff

Discharge in each grid-cell

ReservoirVolume

Volume of each reservoir in \(m^3\)

The total waterlevel is obtained by adding the two water levels.

wflow_routing module documentation

Run the wflow_routing model..

usage

wflow_routing [-h][-v level][-F runinfofile][-L logfile][-C casename][-R runId]
      [-c configfile][-T last_step][-S first_step][-s seconds][-l loglevel]


-X: save state at the end of the run over the initial conditions at the start

-T: Set end time of the run: yyyy-mm-dd hh:mm:ss

-S: Set start time of the run: yyyy-mm-dd hh:mm:ss

-s: Set the model timesteps in seconds

-I: re-initialize the initial model conditions with default

-C: set the name  of the case (directory) to run

-R: set the name runId within the current case

-L: set the logfile

-c: name of wflow the configuration file (default: Casename/wflow_routing.ini).

-h: print usage information

-P: set parameter change string (e.g: -P 'self.FC = self.FC * 1.6') for non-dynamic variables

-p: set parameter change string (e.g: -P 'self.Precipitation = self.Precipitation * 1.11') for
    dynamic variables

-l: loglevel (must be one of DEBUG, WARNING, ERROR)
wflow_routing.main(argv=None)

Perform command line execution of the model.

wflow_routing.usage(*args)