dynwavestate,dynwaveflux,lookuppotential,lookupstate

state = dynwavestate(DynWaveProfileTable,
                     sectionID,
                     ldd,
                     prevState,
                     inflow,
                     bottomlevel,
                     roughness,
                     segmentlength,
                     nrtimeslices,
                     timestepinseconds,
                     constantState)

flux = dynwaveflux(DynWaveProfileTable,
                   sectionID,
                   ldd,
                   prevState,
                   inflow,
                   bottomlevel,
                   roughness,
                   segmentlength,
                   nrtimeslices,
                   timestepinseconds,
                   constantState);

potential = lookuppotential(DynWaveProfileTable,
                            sectionID,
                            bottomlevel,
                            segmentlength,
                            state)

state = lookupstate(DynWaveProfileTable,
                    sectionID,
                    bottomlevel,
                    segmentlength,
                    potential)

Argument

Type

[Units] Remarks

DynWaveProfileTable

ascii table

chanel profile table, described in text below

sectionID

spatial, nominal,

section IDs of the channels

ldd

spatial, ldd

prevState

scalar, spatial

initial state (state of the previous timestep) [m3]

inflow

scalar, spatial

lateral inflow into the channel segment [m3/channellength/sec]

bottomLevel

scalar, spatial

bottom levels of channel segments [m]

roughnes

scalar, spatial
[-] > 0, if manning: Manning coefficient

if chezy: Chezy coefficient

segmentlength

scalar, spatial

actual length of the river segment [m]

nrtimeslices

scalar, spatial non-spatial

number of timeslices in dynwave iterator. The value at the pit defines the nrtimeslices for its whole catchment.

timestepinseconds

non-spatial scalar

length of the pcraster timestep [sec]

constantState

boolean,spatial non-spatial

keep the state constant in the locations which are ‘True’. This allows for defining constant heads (translated in states) at outlet points.

state

scalar, spatial

[m3]

flux

scalar, spatial

[m3]

Global options:

The algortithm is influenced by the options --manning (default) or --chezy, to select the dynamic flow equation (see below for these equations).

Operation

The dynamic wave is a hydraulic routing function that calculates flows through the channel network as a function of water level slope, channel profile and channel roughness characteristics (as opposed to kinematic routing, which takes channel bottom slope as an approximation of water level slope). For a detailed discussion of the dynamic wave is referred to the hydrology and hydraulics handbooks (e.g. Chow), here it is assumed that the concepts of the dynamic wave are known.

In hydrology, we commonly do not deal with the state of water (the amount of water) in any given cell in a given timestep, but instead we deal with water levels or water depths. Water levels are the driving forces of the dynamic wave, but for a correct water massbalance, water levels have to be converted to water levels (potentials) and vice-versa. The relation between water levels and water amounts is determined by the channel profile and section length.

This description of the channel profile is stored in a DynWaveProfileTable. The profile table is an ASCII-formatted table, with each line listing (in this order):

  1. cross section ID

  2. water depth [m]

  3. cross section [m2]

  4. wetted perimeter [m]

Note that the DynWaveProfileTable lists water depth above bottomlevel, and bottomlevel itself is no part of the DynWaveTable.

Example DynWaveProfileTable:

1 0 0   0.0
1 1 1   2.82842712475
1 2 4   5.65685424949
1 3 9   8.48528137424
1 4 16 11.313708499
1 5 25 14.1421356237
1 6 36 16.9705627485
….
2 0 0   0.0
2 1 1   2.82842712475
2 2 4   5.65685424949
2 3 9   8.48528137424
2 4 16 11.313708499
2 5 25 14.1421356237
2 6 36 16.9705627485
…

The conversion between water level (a potential) and water state (an amount) is performed with the functions lookuppotential and lookupstate.

The dynamic wave is in PCRaster solved as an iterative process. This process is defined by the number of iterations (nrtimeslices) and the length of one PCRaster timestep. The number of iterations is a positive integer, and defines the internal iterations for the dynamic wave. Because PCRaster does not explicitly defines the length of the timestep and the units of the time parameters. For the solving the dynamic wave, this information is necessary, so the length of the PCRaster timestep has to be defined in the dynamic wave functions.

The current implementation of the dynamic wave function in PCRaster is not a complete dynamic wave formulation according to the summary of the Saint Venant equations as discussed by Ven Te Chow. The implementation currently consists of the friction force term, the gravity force term and the pressure force term and should therefore be correctly characterised as a diffusion wave formulation. The equations are solved as an explicit, finite forward-difference scheme. A straightforward iteration using an Euler solution scheme is used to solve these equations. These characteristics can not guarantee a stable solution (if the number of iterations is choosen too low, results may show oscillations), but are relatively fast and can be used to maintain a closed water balance.

Pits are considered by the dynamic wave as a closed boundary. Pits will thus act as a ‘storage of water’, meaning water will enter the pit cell through the fluxes from upstream cells, but no outflow from the pit is defined and the water remains in the pit. Over time, the pit thus accumulates water, and this water will eventually build a high enough potential to prevent other water from flowing into the cell.

In reality, two types of boundary conditions can be defined. The one boundary condition is a constant flow boundary condition, the other is a constant head boundary. Both boundary conditions can be implemented in the dynwave as follows.

Constant flow boundary condition in pits

A straightforward implementation of a constant flow boundary in the dynamic wave is to explicitly extract this flow from the pit. Model code implementing this behavior is:

binding
  constantOutflow = ….;
dynamic
  state,flux = dynwavestate,dynwaveflux(DynWaveProfileTable,
                                        sectionID,
                                        lddMap,
                                        prevState,
                                        inflow,
                                        bottomlevel,
                                        roughness,
                                        segmentlength,
                                        10,
                                        1000,
                                        0);

  # Find the pit and subtract a constant value from state
  prevState = if(pit(lddmap) ne 0 then state - constantOutflow else state);

Constant head boundary condition in pits

A straightforward implementation of a constant head boundary in the dynamic wave is to explicitly set the pit to this head, using the function lookupstate to go from a head to a state. Model code implementing this behavior is:

binding
  constantHead = ….;
dynamic
  state,flux = dynwavestate,dynwaveflux(DynWaveProfileTable,
                sectionID,
                lddMap,
                prevState,
                inflow,
                bottomlevel,
                roughness,
                segmentlength,
                10,
                timeinsecs,
                constantHead);

# Find the pit and subtract a constant value from state
prevState = if(pit(lddmap) ne 0 then lookupstate(DynWaveProfileTable,
                                                 sectionID,
                                                 bottomlevel,
                                                 segmentlength,
                                                 potential)
                                                 else state);

Examples

  1. • pcrcalc
    #! –manning
    binding
    State2Manning = State2Manning.map;
    Flux2Manning = Flux2Manning.map;
    Table = Table.txt;
    ProfileId2 = ProfileId2.map;
    Ldd = Ldd.map;
    OldState2 = OldState2.map;
    initial
    report State2Manning, Flux2Manning = dynwavestate,dynwaveflux(
    Table,
    ProfileId2, # with a MV
    Ldd,
    OldState2,
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    );

    • python
    setglobaloption(“manning”)
    Table = “Table.txt”
    ProfileId2 = readmap(“ProfileId2.map”)
    Ldd = readmap(“Ldd.map”)
    OldState2 = readmap(“OldState2.map”)

    State2Manning = dynwavestate(
    Table,
    ProfileId2, # with a MV
    Ldd,
    OldState2,
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )
    Flux2Manning = dynwaveflux(
    Table,
    ProfileId2, # with a MV
    Ldd,
    OldState2,
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )

    State2Manning.map

    Flux2Manning.map

    Table.txt

    ProfileId2.map

    Ldd.map

    OldState2.map
    ../../_images/dynamicwave_State2Manning.png ../../_images/dynamicwave_Flux2Manning.png 
    14 0 0 0.0
14 1 1 2.82842712475
14 2 4 5.65685424949
14 3 9 8.48528137424
14 4 16 11.313708499
14 5 25 14.1421356237
14 6 36 16.9705627485
.....................
14 996 992016 2817.11341625
14 997 994009 2819.94184337
14 998 996004 2822.7702705
14 999 998001 2825.59869762
    		../../_images/dynamicwave_ProfileId2.png ../../_images/accu_Ldd.png ../../_images/dynamicwave_OldState2.png

  2. • pcrcalc
    #! –chezy
    binding
    State2Chezy = State2Chezy.map;
    Flux2Chezy = Flux2Chezy.map;
    Table = Table.txt;
    ProfileId2 = ProfileId2.map;
    Ldd = Ldd.map;
    OldState2 = OldState2.map;
    initial
    report State2Chezy, Flux2Chezy = dynwavestate,dynwaveflux(
    Table,
    ProfileId2, # with a MV
    Ldd,
    OldState2,
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )
    ;

    • python
    setglobaloption(“chezy”)
    Table = “Table.txt”
    ProfileId2 = readmap(“ProfileId2.map”)
    Ldd = readmap(“Ldd.map”)
    OldState2 = readmap(“OldState2.map”)

    State2Chezy=dynwavestate(
    Table,
    ProfileId2, # with a MV
    Ldd,
    OldState2,
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )
    Flux2Chezy=dynwaveflux(
    Table,
    ProfileId2, # with a MV
    Ldd,
    OldState2,
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )

    State2Chezy.map

    Flux2Chezy.map

    Table.txt

    ProfileId2.map

    Ldd.map

    OldState2.map
    ../../_images/dynamicwave_State2Chezy.png ../../_images/dynamicwave_Flux2Chezy.png 
    14 0 0 0.0
14 1 1 2.82842712475
14 2 4 5.65685424949
14 3 9 8.48528137424
14 4 16 11.313708499
14 5 25 14.1421356237
14 6 36 16.9705627485
.....................
14 996 992016 2817.11341625
14 997 994009 2819.94184337
14 998 996004 2822.7702705
14 999 998001 2825.59869762
    		../../_images/dynamicwave_ProfileId2.png ../../_images/accu_Ldd.png ../../_images/dynamicwave_OldState2.png

  3. • pcrcalc
    #! –manning
    binding
    StateManning = StateManning.map;
    FluxManning = FluxManning.map;
    Table = Table.txt;
    Ldd = Ldd.map;
    Material = Material.map;
    initial
    report StateManning, FluxManning = dynwavestate,dynwaveflux(
    Table,
    14,# profileId
    Ldd,
    Material, # oldState
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )
    ;

    • python
    setglobaloption(“manning”)
    Table = “Table.txt”
    Ldd = readmap(“Ldd.map”)
    Material = readmap(“Material.map”)

    StateManning = dynwavestate(
    Table,
    14,# profileId
    Ldd,
    Material, # oldState
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )
    FluxManning = dynwaveflux(
    Table,
    14,# profileId
    Ldd,
    Material, # oldState
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )

    StateManning.map

    FluxManning.map

    Table.txt

    Ldd.map

    Material.map
    ../../_images/dynamicwave_StateManning.png ../../_images/dynamicwave_FluxManning.png 
    14 0 0 0.0
14 1 1 2.82842712475
14 2 4 5.65685424949
14 3 9 8.48528137424
14 4 16 11.313708499
14 5 25 14.1421356237
14 6 36 16.9705627485
.....................
14 996 992016 2817.11341625
14 997 994009 2819.94184337
14 998 996004 2822.7702705
14 999 998001 2825.59869762
    		../../_images/accu_Ldd.png ../../_images/accufraction_Material.png

  4. • pcrcalc
    #! –chezy
    binding
    State3Chezy = State3Chezy.map;
    Flux3Chezy = Flux3Chezy.map;
    Table = Table.txt;
    Ldd = Ldd.map;
    OldState3 = OldState3.map;
    initial
    report State3Chezy, Flux3Chezy = dynwavestate,dynwaveflux(
    Table,
    14,# profileId
    Ldd,
    OldState3,
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )
    ;

    • python
    setglobaloption(“chezy”)
    Table = “Table.txt”
    Ldd = readmap(“Ldd.map”)
    OldState3 = readmap(“OldState3.map”)

    State3Chezy = dynwavestate(
    Table,
    14,# profileId
    Ldd,
    OldState3,
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )
    Flux3Chezy = dynwaveflux(
    Table,
    14,# profileId
    Ldd,
    OldState3,
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )

    State3Chezy.map

    Flux3Chezy.map

    Table.txt

    Ldd.map

    OldState3.map
    ../../_images/dynamicwave_State3Chezy.png ../../_images/dynamicwave_Flux3Chezy.png 
    14 0 0 0.0
14 1 1 2.82842712475
14 2 4 5.65685424949
14 3 9 8.48528137424
14 4 16 11.313708499
14 5 25 14.1421356237
14 6 36 16.9705627485
.....................
14 996 992016 2817.11341625
14 997 994009 2819.94184337
14 998 996004 2822.7702705
14 999 998001 2825.59869762
    		../../_images/accu_Ldd.png ../../_images/dynamicwave_OldState3.png

  5. • pcrcalc
    #! –chezy
    binding
    StateChezy = StateChezy.map;
    FluxChezy = FluxChezy.map;
    Table = Table.txt;
    Ldd = Ldd.map;
    Material = Material.map;
    initial
    report StateChezy, FluxChezy = dynwavestate,dynwaveflux(
    Table,
    14, # profileId
    Ldd,
    Material, # oldState
    0, # inflow
    0, # bottomLevel
    0.1, # roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )
    ;

    • python
    setglobaloption(“chezy”)
    Table = “Table.txt”
    Ldd = readmap(“Ldd.map”)
    Material = readmap(“Material.map”)

    StateChezy = dynwavestate(
    Table,
    14, # profileId
    Ldd,
    Material, # oldState
    0, # inflow
    0, # bottomLevel
    0.1, # roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )
    FluxChezy = dynwaveflux(
    Table,
    14, # profileId
    Ldd,
    Material, # oldState
    0, # inflow
    0, # bottomLevel
    0.1, # roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )

    StateChezy.map

    FluxChezy.map

    Table.txt

    Ldd.map

    Material.map
    ../../_images/dynamicwave_StateChezy.png ../../_images/dynamicwave_FluxChezy.png 
    14 0 0 0.0
14 1 1 2.82842712475
14 2 4 5.65685424949
14 3 9 8.48528137424
14 4 16 11.313708499
14 5 25 14.1421356237
14 6 36 16.9705627485
.....................
14 996 992016 2817.11341625
14 997 994009 2819.94184337
14 998 996004 2822.7702705
14 999 998001 2825.59869762
    		../../_images/accu_Ldd.png ../../_images/accufraction_Material.png

  6. • pcrcalc
    #! –manning
    binding
    State3Manning = State3Manning.map;
    Flux3Manning = Flux3Manning.map;
    Table = Table.txt;
    Ldd = Ldd.map;
    OldState3 = OldState3.map;
    initial
    report State3Manning, Flux3Manning = dynwavestate,dynwaveflux(
    Table,
    14,# profileId
    Ldd,
    OldState3,
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )
    ;

    • python
    setglobaloption(“manning”)
    Table = “Table.txt”
    Ldd = readmap(“Ldd.map”)
    OldState3 = readmap(“OldState3.map”)

    State3Manning = dynwavestate(
    Table,
    14,# profileId
    Ldd,
    OldState3,
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )
    Flux3Manning = dynwaveflux(
    Table,
    14,# profileId
    Ldd,
    OldState3,
    0,# inflow
    0,# bottomLevel
    0.1,# roughness
    1,# segmentLength
    10,# nrTimeSlices
    1, # timestepInSecs
    0 # constantState
    )

    State3Manning.map

    Flux3Manning.map

    Table.txt

    Ldd.map

    OldState3.map
    ../../_images/dynamicwave_State3Manning.png ../../_images/dynamicwave_Flux3Manning.png 
    14 0 0 0.0
14 1 1 2.82842712475
14 2 4 5.65685424949
14 3 9 8.48528137424
14 4 16 11.313708499
14 5 25 14.1421356237
14 6 36 16.9705627485
.....................
14 996 992016 2817.11341625
14 997 994009 2819.94184337
14 998 996004 2822.7702705
14 999 998001 2825.59869762
    		../../_images/accu_Ldd.png ../../_images/dynamicwave_OldState3.png