Navigation

dynwavestate,dynwaveflux

dynwave

Dynamic wave equation as state and flux.

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 chezy: [m1/2sec-1]
if manning: [m1/3sec-1]
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 --chezy (default) or --manning, 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

Table Of Contents

Previous topic

downstreamdist

Next topic

dynamicwave

Navigation

© Copyright 1990-2013, PCRaster R&D Team. Created using Sphinx 1.1.3.