Index: doc/wflow_routing.rst =================================================================== diff -u -r0fc7439329d0ea7ec14ea9360d97a18ebbcf82c0 -r3f28b3d315ef8871a89ef8575ad2594dd7410b8b --- doc/wflow_routing.rst (.../wflow_routing.rst) (revision 0fc7439329d0ea7ec14ea9360d97a18ebbcf82c0) +++ doc/wflow_routing.rst (.../wflow_routing.rst) (revision 3f28b3d315ef8871a89ef8575ad2594dd7410b8b) @@ -1,4 +1,4 @@ -THe wflow_routing Model +The wflow_routing Model ======================= @@ -8,38 +8,112 @@ and a floodplainwidth to the configuration the model can also include estimated flow over a floodplain. Method -====== +------ - """ - Updates the kinematic wave reservoir. Should be run after updates to Q +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. - WL = A * pow(Q,B)/Bw - WL/A * Bw = pow(Q,B) - Q = pow(WL/A * Bw,1/B) - """ - # TODO: This is not correct, need to subtract Channel runoff - 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) +.. figure:: _images/wflow_routing.png - self.WaterLevelCH = self.AlphaCh * pow(self.Qchannel, self.Beta) / (self.Bw) - self.WaterLevelFP = self.AlphaFP * pow(self.Qfloodplain, self.Beta) / (self.Bw + self.Pfp) - #self.WaterLevelFP = max(0.0,self.WaterLevelCH - self.bankFull) - #self.Alpha * pow(self.SurfaceRunoff-self.Qchannel, self.Beta) / (self.Bw + self.Pch) + self.bankFull) - self.Qtot = pow(self.WaterLevelCH/self.AlphaCh * self.Bw,1.0/self.Beta) + pow(self.WaterLevelFP/self.AlphaFP * self.Pfp,1.0/self.Beta) +First the maximum bankfull flow for each cell is determined using: - # wetted perimeter (m) - self.Pch = self.wetPerimiterCH(self.WaterLevelCH,self.Bw) - self.Pfp = self.wetPerimiterFP(self.WaterLevelFP,self.floodPlainWidth) - # Alpha - self.AlphaFP = self.AlpTerm * pow(self.Pfp, self.AlpPow) - self.AlphaCh = self.AlpTerm * pow(self.Pch, self.AlpPow) - self.Alpha = self.AlpTerm * pow(self.Pch + self.Pfp, self.AlpPow) - self.OldKinWaveVolume = self.KinWaveVolume - self.KinWaveVolume = self.WaterLevelCH * self.Bw * self.DCL +.. math:: 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: + +.. math:: H_{ch} = \alpha{_{ch}} {Q_{ch}}^{\beta}/ Bw + +Where :math:`H_{ch}` is the water level for the channel, :math:`\alpha{_{ch}}` is the kinematic wave coefficient for the channel, +:math:`Q_{ch}` is the discharge in the channel and :math:`Bw` is the width of the river. + +If the water level is above bankfull the water level on the floodplains is calculated as follows: + +.. math:: H_{fp} = \alpha{_{fp}} {Q_{fp}}^{\beta}/ (Bw + P_{fp}) + +where :math:`H_{ch}` is the water level on the floodplain, :math:`Q_{fp}` is the discharge on the floodplain +and :math:`P_{fp}` is the wetted perimiter of the floodplain and :math:`\alpha{_{fp}}` is the kinematic wave +coefficient for the floodplain, + +The wetted perimiter of the channel, :math:`P_{ch}`, is determined by: + +.. math:: P_{ch} = 2.0 H_{ch} + Bw + +The wetted perimiter of the floodplain is defined as follows: + +.. math:: + + N = max(0.0001,1.0/(1.0 + exp(-c * H_{fp}) - 0.5) * 2.0 + + P_{fp} = N W_{fp} + +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 :math:`W_{fp}`. + +The :math:`\alpha` for the channel and floodplain are calculated as follows: + + +pow((self.NFloodPlain / (sqrt(self.SlopeDCL))), self.Beta) + +.. math:: + + \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} + +In which slope is the flope of the river bed and floodplain and :math:`n_{ch}` and :math:`n_{fp}` represent +the manning's n for the channel and floodplain respectively. + +A compound :math:`\alpha` is estimated by first calculating a compound n value :math:`n_{total}`: + +.. math:: + + n_{total} = (P_{ch}/P_{total} n_{ch}^{3/2} + P_{fp}/P_{total} n_{fp}^{3/2})^{2/3} + + + + + + + + + +""" +Updates the kinematic wave reservoir. Should be run after updates to Q + +WL = A * pow(Q,B)/Bw +WL/A * Bw = pow(Q,B) +Q = pow(WL/A * Bw,1/B) +""" +# TODO: This is not correct, need to subtract Channel runoff +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 = self.AlphaFP * pow(self.Qfloodplain, self.Beta) / (self.Bw + self.Pfp) + +#self.WaterLevelFP = max(0.0,self.WaterLevelCH - self.bankFull) +#self.Alpha * pow(self.SurfaceRunoff-self.Qchannel, self.Beta) / (self.Bw + self.Pch) + self.bankFull) +self.Qtot = pow(self.WaterLevelCH/self.AlphaCh * self.Bw,1.0/self.Beta) + pow(self.WaterLevelFP/self.AlphaFP * self.Pfp,1.0/self.Beta) + +# wetted perimeter (m) +self.Pch = self.wetPerimiterCH(self.WaterLevelCH,self.Bw) +self.Pfp = self.wetPerimiterFP(self.WaterLevelFP,self.floodPlainWidth) +# Alpha +self.AlphaFP = self.AlpTerm * pow(self.Pfp, self.AlpPow) +self.AlphaCh = self.AlpTerm * pow(self.Pch, self.AlpPow) +self.Alpha = self.AlpTerm * pow(self.Pch + self.Pfp, self.AlpPow) +self.OldKinWaveVolume = self.KinWaveVolume +self.KinWaveVolume = self.WaterLevelCH * self.Bw * self.DCL + + Dependencies ------------ T @@ -133,8 +207,23 @@ 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 + +The total waterlevel is obtained by adding the two water levels. + + wflow_routing module documentation ----------------------------------