departement
Mobiliteit en Openbare Werken
Development of conceptual models for an integrated catchment management Subreport 1: LITERATURE REVIEW OF CONCEPTUAL MODELSTRUCTURES
00_131
WL Rapporten
Vlaamse overheid
Development of conceptual models for an integrated catchment management Subreport 1 – Literature review of conceptual modelstructures Meert, P.; Nossent, J.; Vanderkimpen, P.; Pereira, F.; Delgado, R.; Mostaert, F.
November 2014
WL2014R00_131_1
F-WL-PP10-1 Versie 04 GELDIG VANAF: 12/11/2012
This publication must be cited as follows:
Meert, P.; Nossent, J.; Vanderkimpen, P.; Pereira, F.; Delgado, R.; Mostaert, F. (2014). Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures. Version 4.0. WL Rapporten, 00_131. Flanders Hydraulics Research: Antwerp, Belgium.
Waterbouwkundig Laboratorium Flanders Hydraulics Research
Berchemlei 115 B-2140 Antwerp Tel. +32 (0)3 224 60 35 Fax +32 (0)3 224 60 36 E-mail:
[email protected] www.waterbouwkundiglaboratorium.be
Nothing from this publication may be duplicated and/or published by means of print, photocopy, microfilm or otherwise, without the written consent of the publisher.
F-WL-PP10-1 Versie 04 GELDIG VANAF: 12/11/2012
Document identification Title:
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Customer:
Waterbouwkundig Laboratorium
Keywords (3-5):
Conceptuele modellering, waterbeheer, waterkwantiteit
Text (p.):
83
Confidentiality:
☐ Yes
Ref.:
WL2014R00_131_1
Appendices (p.):
/ ☐ Customer
Exceptions:
☐ Internal ☐ Flemish government Released as from: / ☒ No
☒ Available online
Approval Author
Reviser
Project Leader
Meert, P.
Nossent, J. Vanderkimpen, P.
Pereira, F.
Research & Consulting Manager
Head of Division Mostaert, F.
Verwaest, T.
Revisions Nr.
Date
Definition
Author(s)
1.0
24/01/2014
Concept version
Meert, P.
2.0
05/08/2014
Substantive revision
Nossent, J.; Vanderkimpen, P.
3.0
16/10/2014
Revision customer
Pereira, F.
4.0
18/11/2014
Final version
Pereira, F.
Abstract The evolution towards integral and integrated catchment modelling is impeded by a large number of factors, whereof the two most important are the incompatibility of the existing models and the computational effort needed to run simulations. Incompatibilities arise because of differences between software suppliers, time and space scales, the level of detail, the intentions of the different models and many more. Calculation times of full hydrodynamic models limit the use of integrated catchment models for flood prediction systems, real-time control, uncertainty analysis and other applications. The intention of the present study is to develop and apply a methodology for integrated catchment modelling based on the use of conceptual models for each subcomponent. These conceptual models will replace the existing (hydrodynamic) models for simulating river courses, sewer systems, estuaries, rainfall-runoff, groundwater hydrology and water quality components. Conceptual models can be situated somewhere between empirical and physically-based models: they have a limited physical basis but need to be calibrated with measurements or model results. This report summarizes the results of the first work package of this study. A literature review was set up to identify the possible conceptual model structures that can be used for modelling every subcomponent of the integrated model. The overview contains a description of the derivation, algorithms and application of a large number of conceptual models. One of the key elements here was the possibility of these models to conduct longterm simulations in a very short period. A generalization of the selected model structures was then carried out to allow an application of one technique to multiple (sub)components of the system. This has led to a proposal of a limited number of conceptual model structures that will be used further on in the study for constructing the integrated catchment model F-WL-PP10-1 Versie 04 GELDIG VANAF: 12/11/2012
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Contents 1.
2.
3.
4.
5.
6.
7.
Literatuuroverzicht conceptuele modelstructuren ..................................................................................... 1 1.1.
Inleiding ............................................................................................................................................. 1
1.2.
Neerslag afvoer modellen ................................................................................................................. 2
1.3.
Waterlopen ........................................................................................................................................ 2
1.4.
Afwatering in verstedelijkte gebieden ................................................................................................ 4
1.5.
Estuaria en tijgebonden rivieren ........................................................................................................ 5
Introduction ............................................................................................................................................... 6 2.1.
General framework ............................................................................................................................ 6
2.2.
Model types ....................................................................................................................................... 7
Rainfall – runoff ...................................................................................................................................... 10 3.1.
NAM ................................................................................................................................................ 10
3.2.
PDM ................................................................................................................................................ 13
3.3.
SRM ................................................................................................................................................ 16
3.4.
VHM ................................................................................................................................................ 17
3.5.
Hydromax ........................................................................................................................................ 18
3.6.
Conclusion ...................................................................................................................................... 20
River Courses ......................................................................................................................................... 21 4.1.
Channel flow routing ....................................................................................................................... 21
4.2.
Reservoir routing ............................................................................................................................. 36
4.3.
Application ....................................................................................................................................... 43
Urban drainage systems......................................................................................................................... 45 5.1.
Components of the urban water system ......................................................................................... 45
5.2.
Software packages for integrated urban catchment modelling ....................................................... 60
5.3.
Discussion ....................................................................................................................................... 64
Estuaries................................................................................................................................................. 66 6.1.
Muskingum water-stage methods ................................................................................................... 66
6.2.
CALAM – model .............................................................................................................................. 68
6.3.
Idealized models ............................................................................................................................. 71
6.4.
Response surface models ............................................................................................................... 73
6.5.
Conclusion ...................................................................................................................................... 74
References ............................................................................................................................................. 76
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
I
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
List of tables Table 2-1. Schematic overview of different model types. After (Willems, 2000).............................................. 8
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
II
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
List of figures Figure 2-1. Scheme used for building a conceptual model based on the results of a process-based model. . 9 Figure 3-1. Structure of the NAM – model. .................................................................................................... 11 Figure 3-2. Structure of the PDM rainfall-runoff model (a) and storage representation of any point in the catchment (b). ................................................................................................................................................ 13 Figure 3-3. Catchment representation by storage elements of different depth and their associated probability density function (a); Direct runoff production from a population of stores (b). ............................. 14 Figure 3-4. VHM model structure. (Willems, 2014)........................................................................................ 17 Figure 3-5. Structure of the Hydromax forecasting model ............................................................................. 19 Figure 4-1. Derivation of relationships for convex routing method. ............................................................... 22 Figure 4-2. Effects of attenuation by a linear reservoir (red line) and translation by a linear channel (green line) on an inflow hydrograph (blue line). ...................................................................................................... 24 Figure 4-3. Storage-outflow slope used in the modified Att-Kin routing method (SCS, 1986). ...................... 26 Figure 4-4. Hysteresis on rating curve ........................................................................................................... 28 Figure 4-5. Characteristic length for Kalinin – Milyukov method.................................................................... 28 Figure 4-6. Muskingum prism and wedge storage concept. .......................................................................... 30 Figure 4-7. Definition sketch of molecule in a space-time grid (Smith, 1980). .............................................. 32 Figure 4-8. Arrangements of linear reservoirs that represent the transfer function structures: ..................... 33 Figure 4-9. Computational procedure of a multilinear model. After (Kim, 2011). .......................................... 34 Figure 4-10. Storage-discharge relations: linear versus non-linear reservoirs. ............................................. 35 Figure 4-11. Schematic overview of Euler’s method ..................................................................................... 41 Figure 4-12. Results obtained by the selected computation schemes (left) and comparison of the accuracy of various schemes and the dependency on the time step used (right) ........................................................ 41 Figure 4-13. Comparison between Laurenson-Pilgrim method (LP) and fourth order Runge-Kutta method (RK) in terms of CPU time (left) and relative error on peak water level (right), for different values of the time step. ............................................................................................................................................................... 42 Figure 4-14. Transfer function methodology for conceptual models of river courses, ................................... 43 Figure 4-15. Hydraulic structures methodology for conceptual models of river courses, consisting of a river reach (in green) and a floodplain component (in red). ................................................................................... 44 Figure 5-1. Loss model with initial (Pi) and permanent (Pp) losses and net rainfall (Pnet). After (Achleitner, 2006).............................................................................................................................................................. 47 Figure 5-2. Runoff from pervious surface resulting from a block rain event (Solvi, 2007). ............................ 48 Figure 5-3. Unit hydrograph of the rational method ....................................................................................... 50 Figure 5-4. Time-area diagram and areas contributing to the total flow at the outlet point. (Berlamont, 2004) ....................................................................................................................................................................... 51 Figure 5-5. Longitudinal profile of a sewer system on which the static and dynamic storage is indicated. ... 55 Figure 5-6. Piecewise linear ‘dynamic’ approach (with slopes kstat and kdyn) for the storage/throughflow relationship of a small gravitary sewer system .............................................................................................. 55 Figure 5-7. Partially linearized throughflow - storage relationship. After (Kamradt, 2008). ........................... 56
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
III
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 5-8. Comparison of hydrodynamic and conceptual routing methods for an event with extreme backwater conditions. (Sartor, 1999) ............................................................................................................. 57 Figure 5-9. Backflow model withcombiner-splitter approach (Solvi, 2007) .................................................... 58 Figure 5-10. Longitudinal profile of a sewer system on which the static, dynamic and extra storage during the overflow event are indicated (Vaes, 1999). ............................................................................................. 58 Figure 5-11. Simulating the flow in a Combined Sewer Overflow:................................................................. 59 Figure 5-12. Schematic representation of the reservoir modelling system Remuli (Vaes, 1999). ................. 61 Figure 5-13. Elements and processes within the KOSIM-WEST model. (Solvi, 2007). ................................ 62 Figure 5-14. Screenshot of the CITY DRAIN 2.0 block library (Achleitner, 2007b). ...................................... 64 Figure 6-1. Accuracy of water-stage simulating with bidirectional Muskingum water-stage method, during validation phase (Si-min et al., 2009). ........................................................................................................... 68 Figure 6-2. Schematization of the Rhine-Meuse Delta in the CALAM - model, with nodes and segments ... 69 Figure 6-3. Top view of an idealized estuary or embayment. ........................................................................ 72 Figure 6-4. Diagram of a cross section with tidal flats as used in idealized models. ..................................... 72
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
IV
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
1. Literatuuroverzicht conceptuele modelstructuren Nederlandstalige samenvatting 1.1. Inleiding 1.1.1.
Doelstellingen
Dankzij het werk van de verschillende waterbeheerders in Vlaanderen in de voorbije jaren, beschikt de regio over een uitgebreid model instrumentarium ter ondersteuning van de verschillende waterbeheerders en het waterbeleid. Deze modellen beperken zich echter meestal tot individuele deelcomponenten van het watersysteem en hebben nauwelijks interactie met elkaar. Bovendien zijn de studiegebieden van deze modellen meestal overlappend wat resulteert in gebieden die 2 à 3 keer gemodelleerd zijn met verschillende soorten modellen en door verschillende administraties. Zo zijn in Vlaanderen gedetailleerde hydrodynamische modellen opgebouwd voor de bevaarbare en de onbevaarbare waterlopen, met behulp van verschillende software en in handen van verschillende administraties. De meeste stroomgebieden hebben ook meer dan één hydrologisch model, en het rioolstelsel is gemodelleerd met nog andere software. Rekening houdende met de evolutie van kennis en modellen verwacht men eerder dat het gebruik van verschillende soorten modellen nog zal toenemen (waterkwaliteit, sediment, ecologische modellen, waterbalans, enz.) niet alleen in Vlaanderen maar ook in de naburige regio’s en landen (Wallonië, Brussel, Frankrijk en Nederland). Om de huidige evolutie naar een integraal en geïntegreerd waterbeheer op stroomgebiedsniveau voldoende te kunnen ondersteunen is er een integratie nodig van de hogervermelde modellen, rekening houdend met de volgende aspecten: • • • • •
Verschillende doelstelling van de modellen. Verschillende tijds- en ruimteschaal van de processen. Accumulatie van onzekerheden bij elk systeemmodel. Interactie van de modellen, waardoor een verbinding tussen modellen noodzakelijk is. Verschillend detailniveau van de modellen.
Deze “integrale modellering” wordt vandaag aangepakt door de verschillende leveranciers van software met de beperking dat ze compatibel moeten blijven met de huidige en toekomstige producten van die specifieke leverancier, wat in veel gevallen niet compatibel is met de software van andere leveranciers. Vertrekkend vanuit de bovenvermelde noden heeft de huidige studie als hoofdoel: •
De ontwikkeling en toepassing van een methodologie voor de integraalmodellering van het watersysteem gebaseerd op het gebruik van conceptuele modellen.
De methodologie en toepassing moet toelaten om maximaal gebruik te maken van de bestaande modellen van de verschillende deelsystemen, en moet voldoende ‘open’ zijn om elke uitbreiding met andere deelmodellen mogelijk te maken. De systeemontwikkeling beperkt zich in eerste instantie tot de waterbeheersaspecten: waterkwantiteit en fysico-chemische waterkwaliteit. Het systeem wordt voldoende ‘open’ opgebouwd om een uitbreiding naar andere aspecten zoals ecologie mogelijk te maken. Tijdens de ontwikkeling kunnen rivier-oppervlaktewater, grondwater, rioolstelsels en RWZI’s worden beschouwd. Het voorliggende rapport kadert in de eerste deelopdracht van de studie, nl. een identificatie en veralgemening van mogelijke conceptuele modelstructuren. Hierbij is er gekeken naar modellen en technieken die op wereldschaal dikwijls toegepast worden voor de simulatie van de verschillende onderdelen van het geïntegreerde watersysteem, en die bovendien toelaten om langdurige simulaties uit te voeren in een zeer korte termijn. De beschouwde componenten van het watersysteem zijn: waterwegen en rivieren (met overstromingsgebieden), estuaria en tijgebonden rivieren, rioleringsstelsels, oppervlakte hydrologie en waterkwaliteitscomponenten (opgeloste zuurstof, biochemische zuurstofvraag, ammonium, nitraat, fosforcomponenten en chloride).
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
1
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
1.1.2.
Soorten modellen
Hydrologische en hydraulische modellen die de waterhuishouding beschrijven kunnen onderverdeeld worden in drie categorieën: empirische modellen, fysisch-gebaseerde modellen en conceptuele modellen, waarbij de grootste verschillen te vinden zijn in het detaileringsniveau en de fysische basis van de modelstructuur. Empirische modellen beschrijven de modeluitvoer y in functie van de modelinvoer x en enkele fysische grootheden, via eenvoudige wiskundige vergelijkingen. De fysische processen die aan de basis liggen van het verband F tussen x en y worden niet beschreven. Dit type van modellen wordt daarom ook vaak ‘blackbox’ modellen genoemd. Fysisch-gebaseerde modellen vormen de tegenpool van empirische modellen: de fysische processen worden nu wel zo goed mogelijk beschreven door wiskundige vergelijkingen. Deze fysische basis heeft als voordeel dat een extrapolatie naar extreme gebeurtenissen buiten het kalibratie gebied betrouwbare resultaten zal opleveren, wat mogelijks niet het geval is voor empirische modellen. Omwille van de duidelijke beschrijving van het systeem en de tegenstelling met ‘black-box’ modellen, duidt men dit soort modellen ook aan met ‘white-box’ modellen. Conceptuele modellen kunnen gesitueerd worden tussen de ‘black-box’ en de ‘white-box’ modellen en worden daarom ook wel ‘grey-box’ modellen genoemd. Dit type modellen stelt de fysische werkelijkheid voor met behulp van een beperkt aantal processen, die de werkelijkheid aggregeren in tijd en ruimte. De parameters van deze modellen zijn meestal niet te bepalen op basis van directe metingen, maar moeten gekalibreerd worden op basis van optimalisatie technieken. Ze hebben echter wel een min of meer fysische betekenis, wat enkele extra voorwaarden oplegt bij de kalibratie.
1.2. Neerslag afvoer modellen Neerslag afvoer modellen simuleren de hoeveelheid water die de waterlopen bereikt als gevolg van neerslag afstroming over de oppervlakte, de onverzadigde en verzadigde zone van de bodem en het grondwater. De meeste conceptuele neerslag afvoer modellen hebben een structuur die gebaseerd is op het bestaan van een aantal reservoirs, die de berging van water in de verschillende onderdelen van de hydrologische cyclus beschrijven. De parameters van deze modellen beschrijven de reservoirs en hun onderlinge interacties. Neerslag afvoer modellen die op dit moment gebruikt worden door waterbeheerders in Vlaanderen, en bij uitbreiding België, als onderdeel van voorspellings- en waarschuwingssystemen zijn reeds van het conceptuele type. In Vlaanderen maakt Waterbouwkundig Laboratorium (WL) gebruik van het Deense NAM-model, en de Vlaamse Milieu Maatschappij (VMM) gebruikt het Britse PDM-model en het vereenvoudigd SRM-model dat hiervan afgeleid is. Voorspellingen over debieten en waterstanden in het Waalse gedeelte van het Maasbekken gebeuren met HYDROMAX. Al deze modellen zijn in detail verder besproken in hoofdstuk 3. Aangezien deze modellen reeds uitvoerig gebruikt worden en bijgevolg gekalibreerd zijn voor de meeste gebieden is geopteerd om het literatuuroverzicht te beperken tot deze modellen. Dit biedt ook het voordeel dat in de uiteindelijke toepassing van een geïntegreerd stroomgebiedsmodel de bestaande conceptuele neerslagafvoer modellen en hun parameters overgenomen kunnen worden.
1.3. Waterlopen Conceptuele modellering van rivieren en waterlopen gaat terug tot de begindagen van rivier modellering. Door het gebrek aan rekenkracht werden eenvoudige modellen, met een beperkte fysische basis, ontwikkeld om op korte termijn voorspellingen te kunnen maken van waterpeilen en debieten in de waterloop. Deze technieken (en verdere ontwikkelingen ervan) kunnen in dit onderzoek eveneens gebruikt worden. Ze zijn onder te verdelen in twee categorieën, ‘channel routing’ en ‘reservoir routing’, die hieronder verder besproken worden.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
2
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
1.3.1.
Channel routing
Channel routing technieken voor een rivierpand kunnen beschouwd worden als een eenvoudige berekening van het afwaartse debiet, op basis van het opwaartse en afwaartse debiet op de huidige en vorige tijdstap. Paragraaf 4.1 beschrijft een aantal van deze technieken en modellen die meestal gebaseerd zijn op een lineair regressie model. Deze modellen kunnen zo goed als allemaal samengevat worden in het ‘unified coefficient routing model’ dat uit drie termen bestaat, zie vergelijking ( 4.54 ). De modellen verschillen onderling echter wel in hun afleiding, de bepaling van de parameter waarden en hun eventuele fysische betekenis. Het ‘unified coefficient routing model’ kan nog verder veralgemeend worden door niet enkel de debieten op de huidige en vorige tijdstap te beschouwen, maar ook rekening te houden met debieten op eerdere tijdstippen. Dit leidt tot het transfer functie model van vergelijking ( 4.59 ), wat in principe ook een lineair regressie model is. Het kan echter bewezen worden dat de structuur en parameterwaarden van een transferfunctie herleid kunnen worden tot een serie- en/of parallel schakeling van lineaire reservoirs, wat deze methode een fysische basis verleent. Het nadeel van de transferfunctie methode (en de andere lineaire regressie modellen) is dat een in wezen (sterk) niet-lineair systeem gesimuleerd wordt met een lineair model, wat in sommige situaties (bijvoorbeeld bij overstromingen of niet-laminaire stromingen) tot grove fouten kan leiden. Dit kan verholpen worden door gebruik te maken van niet-lineaire transfer functie modellen, waarbij de parameter waarden kunnen variëren in de tijd als functie van de input. In het meest eenvoudige geval hoeft dit slechts voor één parameter doorgevoerd te worden, wat dan overeenkomt met een transformatie van de input. De niet-lineaire transferfunctie methodologie wordt daarom beschouwd als de meest veralgemeende techniek voor de simulatie van de waterbeweging in natuurlijke rivierpanden. De methode biedt een brede waaier aan modellen van zeer eenvoudige met slechts één parameter tot meer complexe waarbij de parameters kunnen variëren in de tijd. Het gebruik van de methode blijft echter beperkt tot rivierpanden met weinig tot geen artificiële invloeden, zoals bv. regelbare stuwen. Deze zullen namelijk opstuwing veroorzaken, waardoor het lineaire reservoir model dat aan de basis ligt van de methode niet langer van toepassing is. 1.3.2.
Reservoir routing
Wanneer er sprake is van belangrijke opstuwingseffecten kan de methode op basis van transferfuncties niet langer gebruikt worden en dient overgeschakeld te worden naar ‘reservoir routing’ technieken. Deze technieken zijn gebaseerd op numerieke oplossingen van de continuïteitsvergelijking ( 4.63 ). Dit is een differentiaalvergelijking die de verandering van het volume, en hieraan gekoppeld de verandering van het waterpeil, in een reservoir beschrijft op basis van de in- en uitgaande debieten. Het uitgaande debiet van dit reservoir kan dan gekoppeld worden aan het volume in het reservoir, meestal via het verband van beide met het waterpeil. Grof geschetst zijn er twee methoden beschikbaar om deze continuïteitsvergelijking op te lossen. Een eerste methode maakt gebruik van de trapeziumregel, waarbij het gemiddelde genomen wordt van zowel het inkomend en uitgaand debiet als het volume, over de vorige en huidige tijdstap. Dit levert twee onbekenden op (volume en uitgaand debiet), die gerelateerd zijn via een (mogelijks sterk) niet-lineair verband, wat het oplossen tijdrovend kan maken. Verschillende iteratieve en niet-iteratie methodes kunnen gebruikt worden om dit niet-lineair, impliciet stelsel op te lossen. Een tweede mogelijkheid is om gebruik te maken van expliciete methodes, waarbij het volume op de huidige tijdstap berekend wordt op basis van gegevens van één of meer eerdere tijdstappen. Ook hier zijn er verschillende oplossingsmethodes voor handen, elk met hun eigen complexiteit en nauwkeurigheid. De meest ideale methode voor ‘reservoir routing’ moet een optimale balans vertonen tussen nauwkeurigheid (en bijgevolg ook het vermijden van instabiliteiten) en efficiëntie, wat bepaald wordt door de complexiteit van de methode en de gebruikte tijdsstap. Een zeer eenvoudige methode kan aanleiding geven tot instabiliteiten of te grote afwijkingen en bijgevolg een kleine tijdstap nodig hebben om dit te vermijden, terwijl een complexere methode nauwkeurigere resultaten kan opleveren en daarom met een grotere tijdstap kan rekenen. Dit is reeds onderzocht in een paar eerdere studies, al wordt hier telkens een beperkt aantal methodes vergeleken en zijn deze niet toegespitst op de toepassing van conceptuele modellering zoals die in deze studie beoogd wordt.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
3
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
1.3.3.
Toepassing
In twee eerdere onderzoeken werd reeds een strategie ontwikkeld waarin de technieken van ‘channel routing’ en ‘reservoir routing’ toegepast worden voor het opstellen van conceptuele modellen voor waterlopen en de aangrenzende waterlichamen zoals overstromingsgebieden en wachtbekkens. Wanneer er geen sprake is van opstuwing veroorzaakt door belangrijke artificiële invloeden, zoals bv. regelbare stuwen, is de transferfunctie methode de meest aangewezen. Structuren met een permanente invloed op de stroming, zoals bruggen en duikers, kunnen eventueel ook (impliciet) opgenomen worden in de methode. In geval van belangrijke opstuwingseffecten is ‘reservoir routing’, waarbij het riviernetwerk opgedeeld wordt in reservoirs, de aangewezen techniek. Deze techniek zal ook gebruikt worden voor wachtbekkens en overstromingsgebieden, aangezien er ook hier meestal een eenduidig verband bestaat tussen volume en waterpeil.
1.4. Afwatering in verstedelijkte gebieden De afvoer van water in verstedelijkte gebieden kan onderverdeeld worden in drie of vier componenten: een opwaarts afstromingsgebied, het rioleringsnetwerk op zichzelf, de rioolwaterzuiveringsinstallatie (RWZI) en de waterlichamen waarin het (gezuiverde) afvalwater uiteindelijk geloosd wordt. De conceptuele modellen voor deze componenten zijn gelijkaardig aan de eerder besproken modellen voor neerslagafvoer en waterlopen. 1.4.1.
Stedelijk afwateringsgebied
Het opwaarts verstedelijkte afstromingsgebied bepaalt de hoeveelheid water die het rioleringsnetwerk binnen stroomt. Dit kan opgedeeld worden in twee delen: droog weer afvoer en neerslag afvoer, die afhankelijk van het type riolering gemengd of gescheiden afgevoerd worden. Droog weer afvoer omvat het afvalwater dat afkomstig is van huishoudelijk, industrieel en commercieel gebruik en is daarom sterk afhankelijk van het aantal inwoners, hun verbruik en het landgebruik in het afwateringsgebied. De droog weer afvoer verschilt sterk doorheen de dag, maar vertoont per dag wel een gelijkaardig profiel, met als gevolg dat de twee meest gebruikte methodes om de droog weer afvoer te simuleren hier dan ook gebruik van maken. Een eerste methode laat de gebruiker toe om de dagelijkse variatie zelf vast te leggen, terwijl de tweede methode deze benadert met behulp van Fourier-reeksen. In beide gevallen kan een extra variatie ten opzichte van het standaardprofiel voorzien worden, via een random number generator. In het geval van rioleringen wordt in neerslag afvoer modellen enkel de oppervlakkige afstroming beschouwd. De andere afvoer componenten worden verondersteld de naburige waterloop te bereiken. De bepaling van de oppervlakkige afstroming gebeurt in twee stappen: een transformatie van de bruto neerslag in een netto neerslag die bijdraagt tot de oppervlakkige afstroming, gevolgd door een routing naar het punt waar de neerslag afstroming het rioleringsnetwerk bereikt. De berekening van de netto neerslag verloopt via sterk vereenvoudigde empirische modellen. Voor de routing van de oppervlakkige afstroming wordt meestal gebruik gemaakt van één van de eenvoudige ‘channel routing’ technieken, zoals het lineaire reservoir model of het Muskingum routing model. 1.4.2.
Rioleringsnetwerk
Het rioleringsnetwerk verzamelt het afval- en regenwater en transporteert dit naar een stroomafwaarts punt, wat een RWZI en/of een ontvangende waterloop is. Drie soorten stromingen in rioleringen kunnen onderscheiden worden: stroming met vrij oppervlak, stroming onder druk wanneer de leidingen vol zijn en opstuwingen wanneer het inkomend debiet groter is dan de doorvoer capaciteit. De methodes voor de conceptuele modellering van de stroming met vrij oppervlak zijn vrij gelijkaardig aan de methodes die eerder besproken werden voor waterlopen. In eerste instantie zijn er de eenvoudige lineaire modellen, met als meest gekende voorbeeld het lineaire reservoir model of een cascade van dergelijke modellen. Deze modellen kunnen opnieuw veralgemeend worden tot de transferfunctie methode, om zoveel mogelijk situaties te kunnen beschrijven. Het rioleringssysteem kan echter een sterk niet-lineair gedrag vertonen, als gevolg van turbulente stromingen, leidingen die onder druk stromen, e.d.. Om dit nietlineaire karakter in rekening te brengen kan gebruik gemaakt worden van multi-lineaire of niet-lineaire modellen. Het gebruik van niet-lineaire transferfuncties, met variërende parameter waarden, kan mogelijks Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
4
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
een oplossing bieden om het niet-lineaire karakter van het systeem te verwerken in de bestaande methode. Dit is tot op heden echter nog niet getest. Gelijkaardig aan wat bij waterlopen werd vastgesteld zijn deze modellen niet of amper in staat om de effecten van opstuwing in rekening te brengen. Een vaak voorkomende techniek is het opleggen van een maximaal doorvoer debiet op een bepaalde locatie in het netwerk. Indien het inkomende debiet hoger is dan dit doorvoer debiet wordt het overtollige water opgeslagen in een extra reservoir. Het reservoir kan pas geledigd worden, wanneer het debiet terug onder de maximale waarde zakt. Een andere methode bestaat erin om de totale berging in het netwerk te berekenen en hieraan dan het doorvoerdebiet te koppelen, cfr. de ‘reservoir routing’ technieken beschreven bij waterlopen. Beide technieken zijn ook toepasbaar voor de berekening van doorvoer- en overloop debiet ter plaatse van overstortlocaties. 1.4.3.
Rioolwaterzuiveringsinstallatie
Conceptuele modellen voor de water huishouding in RWZI’s lijken momenteel niet te bestaan. De meeste modellen zijn dan ook veel meer gericht op de fysische en chemische processen die optreden tijdens de waterzuivering. Op basis van de opbouw van een RWZI (vb. bezinkingsbekkens) wordt verwacht dat reservoir modellen zullen volstaan voor de modellering van de hoeveelheid water in een RWZI. Indien de berging in de riolering beduidend groter is dan de berging in de RWZI, is een expliciete modellering mogelijks niet nodig. 1.4.4.
Toepassing
Een vergelijking van de conceptuele modellen gebruikt voor waterlopen en rioleringssystemen toont aan dat voor beide systemen quasi dezelfde methodes en modellen gebruikt (kunnen) worden. Dit biedt grote voordelen om tot een veralgemeende methode te komen, waarbij beide systemen op dezelfde manier gemodelleerd worden. De voorgestelde methode verloopt als volgt: op locaties die waarschijnlijk onderhevig zijn aan opstuwingen (bv. knijpleidingen en overstorten) wordt gebruik gemaakt van reservoir modellen, terwijl voor de andere delen de transferfunctie methode toegepast wordt. Afhankelijk van de situatie en de modelresultaten dient hier dan gekozen te worden tussen lineaire of niet-lineaire transferfuncties.
1.5. Estuaria en tijgebonden rivieren Estuaria en tijgebonden rivieren verschillen van de eerder besproken waterlopen doordat de waterpeilen en debieten niet alleen bepaald worden door de stroomopwaartse neerslagafstromingsdebieten, maar ook sterk beïnvloed worden door het afwaartse waterpeil. Beide invloeden dienen dus in rekening gebracht te worden in het conceptuele model, wat niet mogelijk is met behulp van de klassieke methodes die toegepast worden op niet-tijgebonden rivieren. Hoofdstuk 6 bespreekt enkele conceptuele en vereenvoudigde modellen die mogelijks gebruikt kunnen worden voor de simulatie van waterpeilen in tijgebonden rivieren en estuaria. Het ‘Muskingum water-stage’ model en het CALAM-model zijn gebaseerd op reservoir modellen gebruikt voor gewone waterlopen, waarbij enkele aanpassingen gemaakt zijn om de invloed van het afwaartse waterpeil in rekening te brengen. Deze kunnen dus beschouwd worden als conceptuele modellen. De resultaten van het eerste model geven aan dat de nauwkeurigheid te wensen overlaat, waardoor dit model waarschijnlijk weinig bruikbaar is. Het CALAM-model daarentegen vertoont wel voldoende nauwkeurige resultaten en kan bovendien gekoppeld worden aan de ‘reservoir routing’ technieken van delen 1.3.2 en 4.2. De twee andere modellen die besproken worden (geïdealiseerde modellen en de responsoppervlak methode) zijn geen conceptuele modellen, maar zijn toch opgenomen in dit rapport omwille van hun eenvoudige structuur en hun eerdere toepassing in Vlaanderen. Geïdealiseerde modellen leveren een analytische oplossing voor de voortplanting van het getij in een estuarium, wanneer het getij voorgesteld wordt door een eenvoudig oscillerende sinus of cosinus functie. Dit soort modellen is bijgevolg eerder geschikt voor lange termijn voorspellingen, eerder dan de korte termijn waarbij stormopzet en de stand van de maan het getij verstoren. De responsoppervlak methode is een empirisch model waarmee getracht wordt om, via een eenvoudige wiskundige vergelijking, het waterpeil op een bepaalde locatie te voorspellen, op basis van een beperkt aantal controlerende variabelen. Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
5
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
2. Introduction 2.1. General framework Over the last decades the different water managers in the Flanders region have developed and constructed a large number of models to support their tasks in water management and policy. Despite the undeniable advantages of these models, problems arise when one tries to combine them: the models are mostly limited to individual components of the water system, can hardly interact with each other and are all computationally demanding. Furthermore, there exists an overlap in the study areas of the models, resulting in areas that are modelled twice or more with different model types and by different water managers or public services. In the Flanders region, detailed hydrodynamic models have been constructed of the navigable and the unnavigable river courses, with different software packages and by different water managers. Most catchments have more than one rainfall-runoff model and sewer systems are modelled with yet another software package. It is expected that the number of models will not decrease over the next years, but quite the opposite: thanks to the increase in computer power and extended knowledge of the different processes (e.g. water quality, sediment transport, ecological models, …) an increase in the number of models can be expected. This increase is of course not limited to the Flanders region, but will also take place in the neighboring regions and countries (Wallonia, Brussels, France and the Netherlands), resulting in yet more (incompatible) models. To support the current evolution towards an integral and integrated catchment modelling and management in an adequate way, an integration of the different hydrodynamic models will be necessary. This integration should bear in mind the following important topics: the different intentions of the original models, the different time and space scales of the modelled processes, the accumulation of uncertainties of each (sub)component, the interaction and connection of the models and the different level of detail of each model. The currenct policy of most software suppliers, regarding integrated catchment modelling, is to allow and support this kind of modelling with the important limitation that the models have to be compatible with their specific current and future modelling products, which are mostly not compatible with the software of other suppliers. Based on the needs of the water managers and the problems concerning the interaction and compatibility of the different models, the main purpose of this study is: The development and application of a methodology for an integrated catchment modelling, based on the use of conceptual models. This methodology should allow to use the existing models of the different components of the water system (rivers and surface waters, groundwater, sewer systems and waste water treatment plants) at a maximum extent and has to be ‘open’ enough to allow an extension with other (sub)components. This study is limited to the current water management aspects: water quantity and physico-chemical water quality, but should allow future extensions to other aspects, like for example ecology. The current study consists of five subtasks or work packages. This report is situated in the first work package of the study, whose main objectives are the identification of different types of conceptual models that are (on a global scale) most commonly used for modelling the various components and aspects considered part of integrated catchment modelling; and a generalization of the listed conceptual models to make them applicable for the different components and aspects of the system. This identification and generalization is based on a literature review, which forms the main subject of this report. The second work package comprises the construction of the chosen model structures and calibration of their parameter values. This should result in a calibration tool that allows a semi-automatic calibration of the conceptual model components, based on the simulation results of fully detailed models. Semi-automatic implies that an interaction with the user is provided to account for expert knowledge. In a third work package the conceptual models and model structures are subjected to an uncertainty analysis. This should allow to quantify the reliability of the results of both the (sub)components and the conceptual integrated catchment model. The uncertainty analysis will focus on different types of uncertainties (input values, model parameter values, model structure, measurement errors, …) and their Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
6
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
propagation, on the way of quantifying and incorporating these uncertainties, and on the influence of one type of uncertainties on the overall model result. A decision support system (DSS) for integrated catchment modelling will be built in the fourth workpackage. This decision support system is regarded as an open software tool that incorporates the results of the previous three work packages: identification and calibration of the conceptual model structures, combined with the techniques for uncertainty analysis. The main intention of the DSS is to use the integrated catchment model for simulating the important components of the water system and hereby supporting the decisions of water managers. In the fifth and last work package the developed techniques, strategies and systems will be demonstrated for two cases: the river Zenne catchment and the river Dender catchment. Both case studies will also be used throughout the study during the execution of the other work packages.
2.2. Model types Different model types are available for simulating the different components of the water system. They range from physically-based models to simplified conceptual and empirical models. The main differences between these models are the level of detail and the physical basis of the model structure. This section gives a brief overview of these three model types. 2.2.1.
Empirical and data driven models
Empirical models attempt to find a relation between a model input x and the physical quantities and model outputs y that are the subject of the model. Because a physical basis for the relation F between x and y is missing for these models, they are often referred to as ‘black-box models’. The models are built and calibrated with simultaneous measurements of x and y, which therefore need to be representative and sufficiently long. Furthermore, the model structure might depend on the period that was selected for calibration, with the consequence that extrapolation outside the calibration range (e.g. to extreme events) might be very inaccurate. Application of empirical models is workable when the amount of data is sufficiently large and the model structure is very simple. (Willems, 2013) Although this is mostly not the case for hydrologic and hydraulic systems, some examples of empirical models for the different components of the water system can be found in literature. Examples are the simplified rainfall-runoff models of section 5.1.1, like the rational and SCS method, and the response surface method for water level prediction in tidal rivers discussed in section 6.4. In the past two decades a new type of black-box models has been increasingly applied: data driven models. These techniques make use of computational intelligence methods (particularly machine learning) in building models. Machine learning is an area of computer science that concentrates on the theoretical foundations of learning from data (Solomatine and Ostfeld, 2008). The technique takes a known set of input data and known responses to the data, and seeks to build a predictor model that generates reasonable predictions for the response to new data. Data driven models have the same problems as empirical models: they provide little physical insight into the system and are therefore likely to be less robust and possibly unreliable outside the calibration range. (Lees, 2000; Lekkas et al., 2001). An overview of the application of data driven models to the different components of the water system can be found in (Lockus et al., 2005) and (Solomatine and Ostfeld, 2008). A limited number of empirical and data driven models (further referred to as black-box models) are considered in this report, because they have been applied in the Flanders region before. The intention is, however, to disregard these black-box models in the next subtasks of the study because of the earlier mentioned drawbacks.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
7
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Table 2-1. Schematic overview of different model types. After (Willems, 2000).
Increasing level of physically-based modelling
Modelling types
2.2.2.
Physicallybased models
Conceptual models
Empirical models Data driven models
White box
Most model parameters can be measured
e.g. Hydrodynamic models
Grey box
Model parameters need calibration (e.g. using e.g. Linear reservoir measurements for model models output variables)
Black box
Also the model structure building depends on the measurements for the model output variables
e.g. Artificial Neural Networks
Physically-based models
Physically based models are the opposite of empirical models, because they do describe the physical processes between model inputs and outputs, with mathematical equations. Since it is impossible to describe and calculate each single one of these physical processes, the models are limited to the main processes that describe the largest part of the relation between x and y. Because of the link with the physical processes underlying the description of the input-output relation, the model structure is transparent. The models are therefore often called ‘mechanistic’ or ‘internally descriptive’ or ‘white box models’. (Willems, 2000; Willems, 2013) The physical processes modelled in physically-based models can be described mathematically on a macroscopic or microscopic scale. The former uses mathematical equations to describe the processes as they can be observed or perceived, whereas the latter tries to model the processes on a microscopic scale that are responsible of the macroscopic observations. The second approach will contribute to a better understanding of the macroscopic processes, but offers too few advantages when modelling on catchment scale and will have a large impact on the calculation times of the model. (Willems, 2013) Physically based models for hydrologic and hydraulic applications are based on two conservation laws, which govern the behavior of a fluid. These laws involve the conservation of mass or volume, known as continuity, and the conservation of momentum or conservation of energy (Zoppou, 2001):
∂A ∂Q + =0 ∂t ∂x
( 2.1 )
1 ∂u u ∂u ∂h q u ⋅ + ⋅ + = S0 − S f − ⋅ g ∂t g ∂x ∂x g h
( 2.2 )
With A the cross-sectional area, Q the discharge, u the average current velocity, g the gravitational acceleration, h the water height or depth, S0 the bed slope and Sf the friction slope. Equations ( 2.1 ) and ( 2.2 ) are also known as the ‘shallow water equations’ or the ‘de Saint-Venant equations’. The continuity equation ( 2.1 ) states that the rate of change in water depth (or cross-sectional area) with time in a slice of the channel equals the net inflow into the slice of channel. The momentum equation ( 2.2 ) expresses that the rate of change in momentum within a slice of the channel is equal to the sum of forces acting on the slice. (Zoppou, 2001) The ‘de Saint-Venant’ equations form a system of quasi-linear, hyperbolic partial differential equations. They can be solved fully or with a simplified form of the momentum equation, when some terms are neglected. It is important to notice that these simplified solutions are only accurate in situations where the Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
8
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
approximations are valid. The full hydrodynamic de Saint-Venant equations (thus with the complete momentum equation) provide an accurate solution under different boundary conditions and a wide range of possible configurations of river networks. (Willems, 2013) A number of different solution schemes exist for fhese full hydrodynamic solutions (implicit, explicit, 1D, 2D, quasi-2D), which all can be time consuming. This issue will become problematic in applications like real-time control or another type of optimization, uncertainty analysis and long term simulations. (Villazon, 2011) 2.2.3.
Conceptual models
Conceptual models attempt to simulate the most important perceived processes in a lumped way: they are aggregated in space and time into a number of key responses. This type of models can be regarded as physically-inspired or quasi-physical, in-stead of physically-based, since they have a structure that is based on a simplified representation of the physical process that take place in reality. (Knight and Shamseldin, 2006) Because the physical reality (and its underlying processes) are less transparent than for a detailed physically-based model, a conceptual model is also called a ‘grey-box model’ (Willems, 2000). Conceptual model parameters usually refer to a collection of aggregated processes and may cover a large number of subprocesses that are not covered by the model structure (Wagener et al., 2004). The parameters can be divided in two classes: those with a direct physical significance that may be determined by measurements or from general knowledge or experience; and those which can not be measured directly but must be calibrated by optimization, usually subjected to physical limits based on their interpretation (Knight and Shamseldin, 2006). The second type of parameters forms the vast majority. In reality, the number of observation points is insufficient to allow an adequate calibration of conceptual models. That is why in this study they will be calibrated based on the previously constructed physicallybased models, so that the information that is available within these models is used to a maximum extent. The suggested procedure to go from the real water system to a conceptual model of it can be summarized as follows (Vanrolleghem et al., 2005): I. II. III. IV. V.
Determine the system under study, its boundaries and the problem to be solved. Collect data on the system to calibrate a complex mechanistic model. Calibrate and validate the complex mechanistic model. Generate data with the complex model to calibrate the surrogate model. Calibrate and validate the surrogate model.
A scheme of the procedure is presented in Figure 2-1. In general, the approach has been applied to build conceptual and data driven models that replicate complex process-based models.
Figure 2-1. Scheme used for building a conceptual model based on the results of a process-based model. After (Velez et al., 2013)
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
9
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
3. Rainfall – runoff Conceptual rainfall-runoff models aggregate the hydrological processes occurring in a catchment in space and time, into a number of key responses (e.g. interception, evaporation, infiltration and groundwater and surface flows). Most of these models have a common global model structure that is based on the existence of a number of storage reservoirs. These reservoirs are used to conceptualise the storage of water at the surface, in the unsaturated and saturated zone or in the groundwater. The model parameters and functions describe the properties of these reservoirs and their mutual interactions. Despite the similarities in model structure of conceptual rainfall-runoff models, a large number of rainfall-runoff models exist, each differing in the number of hydrological processes they want to simulate and the according level of detail. (Wagener et al., 2004; Willems, 2013) The input of rainfall-runoff models consists of two basic meteorological data: rainfall and potential evapotranspiration time series. Rainfall data time series can originate from nearby rain gauges or from weather radar images, which might need some extra downscaling (Decloedt and Willems, 2013). Potential evapotranspiration series are measured in a couple of stations across Flanders. The rainfall-runoff models produce time series of rainfall-runoff discharges that can be used as upstream boundary conditions for the river model. This rainfall-runoff discharge is usually divided into three components: surface runoff, also called ‘overland flow’; a runoff component through the unsaturated zones of the soil, often called ‘interflow’; and a slow flow component via the groundwater, often referred to as ‘baseflow’. The determination of the relative proportion of these three subflows, which changes in time, is one of the main differences between the different conceptual rainfall-runoff models. (Willems, 2013) This chapter gives an overview and description of the conceptual rainfall-runoff models that are frequently used by water managers in Belgium. The overview is limited to these models, since a worldwide overview of conceptual rainfall-runoff models goes beyond the scope of this study. Reference is made to (Beven, 2001), (Moore and Bell, 2001) and (Wagener et al., 2004), among others, for a more extensive review of conceptual rainfall runoff models.
3.1. NAM The NAM – model (“Nedbør-Afstrømnings-Model”, which is Danish for precipitation-runoff-model) was originally developed by the Department of Hydrodynamics and Water Resources at the Technical University of Denmark (Nielsen and Hansen, 1973). It is now a part of the rainfall-runoff module of the MIKE 11 river modelling system (DHI, 2011), where it can be used independently or to represent a catchment that generates inflows for the river modelling system. NAM has been applied to a large number of catchments around the world, representing many different hydrological regimes and climatic conditions. In Belgium, the Hydrologic Information Centre (HIC), which is a part of Flanders Hydraulics Research, uses NAM as the hydrologic component of their flood prediction system (HIC, 2014). The land phase of the hydrological cycle, the rainfall-runoff process, is simulated in NAM with three or more different and interrelated storages that represent different physical elements of the catchment. These three storages are the surface storage, the lower or root zone storage and the groundwater storage. Optionally, extra storages can be defined for the snow storage and an extra baseflow component. In addition to these storages, NAM also allows the incorporation of artificial interventions in the hydrological cycle, such as irrigation and groundwater pumping. In this case time series have to be provided for both interventions. Figure 3-1 shows the structure of the NAM – model. (DHI, 2011)
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
10
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 3-1. Structure of the NAM – model. After (Madsen, 2000)
The four storages in NAM will be discussed here, starting with the surface storage, followed by the deeper ground zones (lower zone and groundwater storage). Snow storage is described last, since this is an optional component without an influence on the behavior of the other storages. The different components and parameters of the model are noted in the same way as they are depicted in Figure 3-1. The following description of the model is largely based on the MIKE 11 reference manual (DHI, 2011). The surface storage represents the water that is intercepted on the vegetation or trapped in depressions in the uppermost, cultivated part of the surface. The amount of water in this storage U increases with the rainfall P and is continuously diminished by evaporative consumption EP and by horizontal leakage, represented by the interflow QIF. When the upper limit Umax of the amount of water in the surface storage is exceeded, the excess water PN will give rise to overland flow QOF, whereas the remainder is diverted as infiltration into the lower zone and the groundwater storage. The lower zone storage describes the soil moisture in the root zone, a soil layer below the surface from which the vegetation can draw water for transpiration. The voids between the groundwater particles in this zone consist of a mixture of air and moisture. The moisture content in the lower zone L controls the amount of water that escapes by (actual) evapotranspiration Ea, that enters the groundwater as recharge G, and that contributes to the interflow and overland flow components. Evapotranspiration can distract water from both the surface storage and the lower zone storage. When the evapotranspiration demand EP is smaller than the surface storage U (in mm), the amount of water in this storage is diminished by the evapotranspiration demand. However, when the amount of water is insufficient (U < EP), the remaining fraction is assumed to be withdrawn by root activity from the lower zone storage at an actual rate Ea:
E a = (E p − U ) ⋅
L Lmax
≥0
( 3.1 )
Where Lmax (in mm) represents the upper limit of the soil moisture and L/Lmax is called the relative soil moisture content of the lower zone storage.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
11
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Overland flow arises when the surface storage spills (U > Umax). QOF denotes the part of the excess water PN that contributes to overland flow. It is assumed to be proportional to PN and to vary linearly with the relative soil moisture content of the lower zone storage:
L / Lmax − TOF CQOF ⋅ ⋅ PN QOF = 1 − TOF 0
for
L / Lmax > TOF
for
L / Lmax ≤ TOF
( 3.2 )
With CQOF the overland flow runoff coefficient and TOF the threshold value for overland flow, both are dimensionless constants between zero and one. The contribution to the interflow QIF is calculated in a very similar way:
L / Lmax − TIF (CKIF )−1 ⋅ ⋅U QIF = 1 − TIF 0
for
L / Lmax > TIF
for
L / Lmax ≤ TIF
( 3.3 )
With CKIF the time constant for interflow, and TIF the root zone threshold value for interflow (0 ≤ TIF ≤ 1). The overland flow and interflow are subsequently routed through two identical linear reservoirs in series (see section 4.1.2), with the reservoir constant CK12. In the case of interflow, CK12 is kept constant, while for the overland flow routing a variable reservoir constant is used:
CK 12 CK = CK 12
OF ⋅ OFmin
−β
for
OF < OFmin
for
OF ≥ OFmin
( 3.4 )
Where OF is the overland flow (in mm/hour), OFmin is the upper limit for linear routing (= 0.4 mm/hour) and β = 0.4 corresponds to the Manning formula for modelling overland flow. The groundwater storage represents the fully saturated part of the ground, beneath the root zone. This storage is filled by infiltrating water from the surface and root zone storage. A part of PN does not contribute to the overland flow QOF, but infiltrates into the lower zone storage. A portion DL (in mm), of the water available for infiltration (PN – QOF), is assumed to increase the moisture content in the lower zone storage. The remaining part G percolates deeper and recharges the groundwater storage:
L / Lmax − TG (P − QOF ) ⋅ G= N 1 − TG 0
for
L / Lmax > TG
for
L / Lmax ≤ TG
( 3.5 )
and:
DL = PN − QOF − G
( 3.6 )
With TG the root zone threshold value for groundwater recharge (0 ≤ TG ≤ 1). The groundwater gives rise to a baseflow BF, which is calculated as the outflow from a linear reservoir with the reservoir constant CKBF. An extra (lower) groundwater storage with a slower response can be defined to achieve a better description of the baseflow, that then consists of two components. When the snow melt component is activated, the precipitation is retained in the snow storage during cold periods and released later on as melt water when the temperature rises above a certain value. It is assumed that the precipitation falls as rain when the air temperature is above a certain base temperature level T0. This is typically somewhere around 0 °C. The snowmelt QS is calculated with a degree-day approach:
⋅ (T − T0 ) C QS = snow 0
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
for for
T > T0 T ≤ T0
WL2014R00_131_1
( 3.7 )
12
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
With Csnow the degree-day coefficient and T the actual temperature. The melted snow is retained in the snow storage as liquid water until this amount exceeds the water retention capacity of the snow storage. The excess melt water PS is routed to the NAM model where it contributes to the surface storage:
Q PS = melt 0
for for
WR ≥ C wr ⋅ S snow WR < C wr ⋅ S snow
( 3.8 )
Where WR is the water retention in the snow storage, Cwr is the water retention coefficient and Ssnow is the snow storage. Notice that the product of Cws and Ssnow yields the water retention capacity of the snow storage. The parameters of the NAM-model represent average values for the entire catchment and can be calibrated manually or with computer-based automatic procedures. Manual calibration implies a trial-and-error parameter value adjustment until satisfactory results are obtained. The automatic calibration routine is based on a multi-objective optimization strategy in which four objectives can be optimized simultaneously. These objectives are expressed by four performance measures: the overall volume error, the overall root mean square error (RMSE) and the average RMSE of peak and low flow events.
3.2. PDM The Probability Distributed Model, or PDM, was developed in the ‘80s of the previous century at the Institute of Hydrology in Wallingford, UK. The model is incorporated in the InfoWorks RS software package, where it can be used as the upstream boundary for the river network (Wallingford, 2012b). In Flanders the PDM is used by the Environment Agency as the rainfall-runoff component of their flood predicting and warning system (VMM, 2014). The PDM consists of three main components: a probability-distributed soil moisture storage component to make a distinction between direct runoff and subsurface runoff; a surface storage component to transform direct runoff into surface runoff; and a groundwater storage component that determines the groundwater component of the total runoff. Runoff production is represented in the soil moisture component as a saturation excess runoff process controlled by the absorption capacity of surface and soil. The model assumes furthermore that the spatial variation of the absorption capacity of different points in a catchment can be described by a probability distribution function. This principle will be explained in more detail, based on (Moore, 2007).
(a)
(b)
Figure 3-2. Structure of the PDM rainfall-runoff model (a) and storage representation of any point in the catchment (b). (Moore, 2007)
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
13
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
(a)
(b)
Figure 3-3. Catchment representation by storage elements of different depth and their associated probability density function (a); Direct runoff production from a population of stores (b). (Moore, 2007)
The starting point of the PDM is the idea that any point in a catchment can be conceptualized as a single storage with a certain capacity c’, that represents the absorption capacity of the soil column at that particular point. The storage is filled by the rainfall P and is emptied by evaporation E, until either the storage fills and spills, resulting in direct runoff q, or the storage dries up because of the evaporation. Figure 3-2 depicts such a single storage, together with the different water flows. The behaviour of the storage can be expressed as:
P − E − (c'− S 0 ) q' = 0
for
P > c'+ E
for
P ≤ c'+ E
( 3.9 )
With S0 the initial amount of water in the storage. Runoff production at every point in the catchment can now be described similarly, but with the assumption that each point has a different storage capacity. The storage capacity c at any point may then be considered as a random variable with probability density function f(c), so that the proportion of the catchment with capacities in the range (c, c+dc) will be f(c).dc. The next step is to construct a water balance for the whole catchment, whose storage capacity is distributed with such a probability density function. Suppose that all the storages, that together form the catchment, can be arranged in order of depth and with their open tops arranged at the same height. This results in a wedge shaped diagram as shown in Figure 3-3 (a). Consider now a rain event of intensity P, during a unit time interval Δt. If the catchment is initially dry so that all storages are empty, then the storages will fill to a depth P unless they are of smaller depth than P when they will fill and spill. At the end of the unit time interval all storages with a capacity lower than P are generating runoff. The upper dark blue triangle in Figure 3-3 (b) depicts the amount of runoff produced by storages of different depths over the unit time interval. The critical capacity below which all storages are full at some time t is denoted as C* ≡ C*(t) (for the simple example presented here, C* = P). The proportion of the basin containing storages of capacity less than or equal to C* is then:
prob(c ≤ C *) = F (C *) =
C*
∫ f (c ) ⋅ dc
( 3.10 )
0
Where F(c) is the distribution function of storage capacity, which is related to the density function f(c) through the relation f(c) = dF(c)/dc. This proportion is also the proportion of the catchment generating runoff, so that the contributing area at time t for a catchment of area A is:
Ac (t ) = F (C * (t )) ⋅ A
( 3.11 )
The instantaneous direct runoff rate per unit area from the catchment is the product of the net rainfall rate, π(t) and the proportion of the catchment generating runoff F(C*(t)):
q (t ) = π (t ) ⋅ F (C * (t ))
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
( 3.12 )
WL2014R00_131_1
14
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
The critical capacity C*, below which all storages are filled, is a function of time and will increase during every wet time interval. If rainfall and potential evaporation are constant during the i’th time interval, so that the net rainfall equals πi = Pi - Ei, then the critical capacity will increase over the time interval according to:
C * (t ) = C * (t ) + π i (t − t )
t ≤ t ≤ t + ∆t
( 3.13 )
During dry periods the water content in the storages decreases because of the potential evaporation. It is assumed during such periods that the water moves between storages of different depths, to equalize the depth of stored water at different points within the basin. Thus at any time all storages will have a water content C*, irrespective of their capacity, unless this is less than C* when they will be full. The evaporation loss is assumed to be dependent on soil moisture by the following function:
S E i' − S (t ) = 1 − max Ei S max
be
( 3.14 )
With Ei’ and Ei respectively the actual and potential evaporation, be an exponent that is usually equal to one or two, S(t) the water storage in the whole basin and Smax the total available storage, given by: ∞
∞
0
0
S max = ∫ c ⋅ f (c ) ⋅ dc =
∫ (1 − F (c )) ⋅ dc = c
( 3.15 )
The water stored in the whole basin is related to the critical capacity C*(t), and in turn to the instantaneous rate of basin production, via the following expression:
S (t ) = =
C *(t )
∞
0
C* t
∫ c ⋅ f (c ) ⋅ dc + C * (t ) ⋅ ∫( )f (c ) ⋅ dc
C *(t )
∫ (1 − F (c )) ⋅ dc
( 3.16 )
0
PDM was designed to provide a toolkit of model algorithms that can be configured to suit a wide variety of catchment response. Analytical solutions of the integrals occurring in the probability-distributed soil moisture storage component (e.g. equation ( 3.16 )) are listed in Appendix A of (Moore, 2007) for a number of possible distribution types. Recharge to groundwater can be introduced as an extra loss of the water in the storages. PDM offers three methods of incorporating this recharge. The first option assumes that the recharge to the groundwater is linearly related to the basin soil moisture content:
d i = k g−1 ⋅ (S (t ) − S t ) g b
( 3.17 )
With kg the drainage time constant, St the threshold storage below which no drainage occurs and bg the exponent of the recharge function (usually 1). In the second method the recharge is dependent on both the soil and groundwater storage, which is useful in catchments with important soil/groundwater interactions:
S (t ) S d i = q sat + max − q sat ⋅ f (t ) ⋅ ∆t S max
( 3.18 )
1 S mαx − S (t ) β g f (t ) = α ⋅ S gmαx 1
( 3.19 )
With:
for
g (t ) < α
for
g (t ) ≥ α
Where Sgmax is the maximum groundwater storage capacity, Sg(t) is the groundwater storage at time t, qsat is the outflow from the groundwater storage when Sg(t) equals Sgmax and f(t) is called the groundwater demand Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
15
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
factor. The third option does not consider soil drainage, but calculates the surface runoff as a fraction α of the direct runoff, while the remaining fraction 1-α goes to the groundwater storage. The total runoff of the catchment, which is the output of the PDM, is divided in surface storage originating from direct runoff and groundwater recharge from soil water drainage that comes from the subsurface storage. The former represents the fast response of the system, while the latter is used for the slow flow paths. Both flow systems can be defined by a variety of non-linear storage reservoirs or by a cascade of two linear reservoirs, with time constants k1 and k2. The Horton-Izzard equation is used for the nonlinear storage model:
dq = a ⋅ (u − q ) ⋅ q b dt
( 3.20 )
With u and q respectively the inflow and outflow of the reservoir. In case of surface runoff, u equals the direct runoff, whereas u will equal the drainage d when routing the slow flow component. Equation ( 3.20 ) is based on the classical continuity equation, combined with a nonlinear storage form of the momentum equation:
q = k ⋅Sm
( 3.21 )
Where S is the volume of water held in storage, k is the storage rate coefficient and m is the store exponent. A cubic form of ( 3.21 ), with m = 3, is usually considered most appropriate to represent the groundwater storage. Solutions of the Horton-Izzard equation ( 3.22 ) are listed in the appendix of (Moore and Bell, 2002). PDM has a total of 14 parameters (see table 1 of (Moore, 2007)), that can be calibrated within a generic Model Calibration Shell environment. This shell allows both automatic optimisation and informal visuallyinteractive parameter estimation. Methods are available to update the PDM with reference to observed flows for real-time flow forecasting applications. This is done by either state-correction (adjusting the values of variables like surface and groundwater stores) or by error prediction (predicting future errors based on an analysis of past errors).
3.3. SRM The two previously discussed rainfall runoff models that are used by the different water managers in the Flanders region (NAM and PDM), are calibrated based on rainfall, evaporation and discharge measurements. This is, however, not possible for ungauged catchments, since they lack discharge data. Furthermore, using the parameter values of a neighboring catchment might lead to strong differences and large errors. The Flemish Environmental Agency developed a new simple model concept that uses the results of a calibrated PDM of a nearby catchment, together with potential runoff coefficients to estimate the surface runoff of an ungauged catchment. This new model concept was called Simplified Runoff Model (SRM) and was implemented in the InfoWorks-RS software (Cabus et al., 2007). The simulation of surface runoff is based on two steps: determination of the effective rainfall and subsequently a routing through two linear reservoirs in series. The mean areal rainfall N is transformed in the effective rainfall N via:
N eff = N ⋅ ϕ ⋅ (1 − SMD )
( 3.22 )
With φ a constant potential runoff coefficient and SMD a soil moisture deficit factor. The value of the potential runoff coefficient is extracted from a table, which was constructed based on a large number of measurements. In this table, the potential runoff coefficient is defined by the soil use, gradient and soil characteristics (texture and geology) of the catchment. The effective rainfall is used as input for a simplified runoff model that ignores the soil moisture content, evaporation, infiltration and drainage to the groundwater. Routing of the runoff is done with two linear reservoirs in series. A time delay can be introduced to the system, to account for the translation time in the rainfall-runoff process (Cabus et al., 2007). Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
16
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
3.4. VHM The VHM-model (“Veralgemeend Hydrologisch Model”, which is a Dutch translation of “Generalized rainfallrunoff model”) was developed at the Hydraulics Section of the University of Leuven (Willems, 2014; Willems et al., 2014). VHM starts from a generalized lumped conceptual model structure, with a number of storage and routing elements (see Figure 3-4), and allows the user to choose and specify the model equations and parameters based on advanced time series analysis techniques. To do so, the response between rainfall and evapotranspiration inputs on the one hand and the river flow outputs on the other hand are first separated in a number of subresponses and submodel components that can be calibrated as independently as possible. The approach to identify the model structure and equations starts from a time series of river flow observations measured downstream of a catchment. A number of analysis techniques is applied to this time series: eliminating the hydraulic river influences, separation of the time series in its subflow components and dividing the series in nearly independent quick and slow flow events (Willems, 2009). The rainfall input is separated into different fractions that contribute to the different subflows by a ‘timevariable distributing valve’ (Figure 3-4). A first valve separates the rainfall input into a portion that goes to the slow flow runoff xsf, a portion that is stored as soil moisture storage xu and a portion that contributes to the quick runoff flow xQF. The quick flow portion can (optionally) be further divided by a second valve in an overland flow portion xOF and an interflow xIF portion.
xQF = f QF ⋅ P xOF = f OF ⋅ P x IF = f IF ⋅ P
x SF = f SF ⋅ P xTF = f TF ⋅ P xU = f U ⋅ P
( 3.23 )
With fi the time variable fractions of the rainfall input P, that all have values between 0 and 1. These fractions have to meet some restrictions to assure rain water continuity:
f QF + f SF + f U = 1 or
xQF + x SF + xU = P
f QF + f SF = f TF f OF + f IF = f QF
xQF + x SF = xTF xQF + x SF = xTF
or or
( 3.24 )
This implies that only three of the six time variable fractions need to be modelled, since the other three can be calculated from the equations of ( 3.24 ).
Figure 3-4. VHM model structure. (Willems, 2014)
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
17
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
VHM assumes, like many other conceptual rainfall-runoff models, that the time variability of the fractions fOF, fIF and fu can be related to the soil moisture content u, or to the relative soil moisture level u/umax, with umax the maximum soil moisture level. Two expressions are available to relate these two time variable quantities.
f i = a i ,1 + a i , 2 ⋅
u
( 3.25 )
u max
Or:
u f i = ai ,1 ⋅ exp ai , 2 ⋅ u max
( 3.26 )
With ai,1 and ai,2 the model parameters for fraction fi. These expressions are similar to the ones used in other rainfall-runoff models: equation ( 3.25 ) for example is identical to equation ( 3.2 ) of the NAM-model, whereas equation ( 3.26 ) can be found in the PDM when an exponential distribution for the absorption capacity is used. Calibration of the subflow fractions is based on the agreement between filtered and modelled peak discharges or total volumes during each subperiod, rather than looking at the agreement for all time steps. This is done to minimize random fluctuations caused by measurement errors or river influences, that may substantially influence the calibration (Willems, 2000). By allowing multiple expressions for the time-variability of the subflow proportions, the VHM-model allows the user to choose the optimal model equations to meet a wide variety of catchment responses. The moisture storage u represents the total surface and soil water storage in the whole catchment, which does not lead to either of the three subflows, but can be subjected to evapo(transpi)ration. Soil moisture storage volumes are estimated by cumulating the rest of the rainfall portions contributing to the soil storage and after subtracting the actual evaporation estimates. The actual evapotranspiration Ea is modelled as a fraction of the total evapotranspiration Ep, which is measured or modelled:
u ⋅ Ep E a = u evap E p
for
u < u evap
for
u ≥ u evap
( 3.27 )
With uevap the soil moisture storage when the evapotranspiration is maximal. Once the fractions of the different subflows are determined, they are transformed into runoff discharges yOF, yIF and ySF by linear reservoir models. One linear reservoir is taken as prior selection, but this can be replaced by a cascade of linear reservoirs to achieve a more detailed and accurate routing. The reservoir constants over these models were previously determined in a subflow filtering process (Willems, 2009). The total rainfall runoff y is then achieved by summing up the three subflows.
3.5. Hydromax Hydromax is a river flow forecasting system for extreme hydrological events in the Meuse river basin and its main tributaries, located in the Walloon part of Belgium. It was developed by SETHY (Service of hydrologic studies of the Walloon Public Administration), together with CESAME (Centre for Systems Engineering and Applied Mechanisms, of the university of Louvain-la-Neuve). The model can provide in real-time short-term predictions of river flows based on rainfall and past river flow measurements, and also long-term flood forecasting based on meteorological forecasts. Hydromax was set up for the first time during the flood of January 1995 and has been running since then, without interruption. During these years the model has proven to be very accurate (Bastin et al., 2009). The hydromax forecasting model consists of a conceptual and a statistical submodel (see Figure 3-5 for a schematic overview of the model components). The conceptual submodel is based on a storage reservoir model that computes the effective rainfall from the mean areal rainfall, while the statistical submodel is a linear ARX transfer function model that transforms the effective rainfall into river discharge predictions. The data are collected with a basic time step of one hour (Bastin et al., 2009).
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
18
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 3-5. Structure of the Hydromax forecasting model
The conceptual reservoir model is used to transform the mean areal rainfall PB into an effective rainfall PN which is supposed to reach the catchment outlet as surface runoff. This mean areal rainfall is estimated from a set of rainfall observations, by means of the Kriging method. The effective rainfall is determined with the mass balance in a reservoir during the time interval Δt (Bastin and Moens, 2000):
PB(t ) = PN (t ) + E1 (t ) + W (t )
( 3.28 )
Where E1(t) is the amount of rainfall that directly evaporates and W(t) represents the part of the precipitated water that will not participate in the runoff, but will be stored in the catchment under various forms: vegetation interception, depression storage, soil moisture, etc.. All the terms of equation ( 3.28 ) are expressed in mm. The storage in the catchment is represented by a reservoir model with inflow W(t) and can be described with another mass balance equation:
S (t ) = S (t − ∆t ) + W (t ) − I (t ) − E 2 (t )
( 3.29 )
Where S(t) denotes the water stored in the catchment, I(t) is the amount of water drained by percolation and E2(t) is the part of the stored water that evaporates during the time interval Δt. All terms are again expressed in mm. The amount of percolation is described by a linear function of the available water stock:
I (t ) = α ⋅ [S (t − ∆t ) + W (t )]
( 3.30 )
With α a percolation parameter. The evaporation terms E1(t) and E2(t) are calculated as:
E1 (t ) = min{PB(t ); ETP (t )} E 2 (t ) = max{0; min[ETP (t ) − PB(t ); S (t − ∆t ) + W (t ) − I (t )]}
( 3.31 )
Where ETP(t) is a periodic forcing input of the model which represents an estimate of the seasonal potential evaporation for the considered catchment (Wery, 1990). The water storage term W(t) is expressed as a function of S(t) and PB(t):
PB(t ) − E1 (t ) W (t ) = [S max − S (t )] ⋅ 1 − exp − β ⋅ S max − S (t )
( 3.32 )
With β a specific runoff parameter that lies between zero and one. Equation ( 3.32 ) is formulated this way so that the amount of stored water does not exceed a physical upper limit Smax and to verify the hydrological principle that the effective rainfall PN increases with both rainfall intensity PB and stored water S. Calibration of the conceptual submodel for the considered catchment thus only involves determining three positive parameters (Smax, α and β). Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
19
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
The statistical submodel computes a forecast
Qˆ (t + h )
at each time step, as a linear combination of past
river flow measurements and past effective rainfall values, with a linear regression (ARX) model (Bastin and Moens, 2000): n
n
i =1
j =1
Qˆ (t + h ) = ∑ ai ⋅ Q(t − (i − 1) ⋅ h ) + ∑ b j ⋅ PN (t − ( j − 1) ⋅ h )
( 3.33 )
With h the prediction horizon expressed as a number of time steps, Q(t - (i - 1).h) river flow measurements at previous time steps and PN(t - (j - 1).h) effective rainfall values calculated with the conceptual submodel and cumulated over h successive time steps. The values of the transfer function coefficients ai and bj are determined from experimental data by linear regression. The number of coefficients a and b are selected using classical statistical tools of system identification theory. The prediction horizon h must be smaller than the natural response time of the catchment and is mostly chosen somewhere between one fifth and one third of the peak time of the unit hydrograph. The statistical submodel is not only used for short-term predictions, but is also applicable to long term forecasting over prediction horizons that are significantly larger than the natural response time of the river basin. This obviously requires to anticipate the future rainfalls by meteorological informations. Such longterm forecasts are computed by iterating the short-term prediction model as: k −1
Qˆ (t + k ⋅ h ) = ∑ ai ⋅ Qˆ (t − (i − 1) ⋅ h ) + bi ⋅ Pˆ N (t − (i − 1) ⋅ h ) i =1
n
m
j =k
j =k
+ ∑ a j ⋅ Q(t − ( j − 1) ⋅ h ) + ∑ b j ⋅ PN (t − ( j − 1) ⋅ h ) Where
Qˆ
are successive iterated river flow forecasts and
( 3.34 )
Pˆ N are effective rainfall forecasts to be provided
by the user.
3.6. Conclusion Five conceptual rainfall runoff models were discussed in this chapter. All of them, except the VHM model, are (widely) used in Belgium by the different water managers and public services in flood and drought prediction and warning systems. An overview of conceptual rainfall-runoff models used in other parts of the world is not included in this chapter, since this goes beyond the scope of this study. The prediction systems of the Hydrologic Information Centre (using NAM), the Flemish Environment Agency (PDM) and the Walloon Public Administration (HYDROMAX) cover most of the Belgian territory. This implies that rainfall-runoff models have been calibrated for nearly all upstream catchments. It is therefore suggested to reuse these existing conceptual models for the rainfall-runoff component of the final integrated catchment model. It is the authors’ belief that there is no gain in both accuracy and simulation speed in developing and/or recalibrating the existing rainfall-runoff models.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
20
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
4. River Courses The use of simplified conceptual models for routing flows in river courses goes back to the early days of river modelling. Before the time of sophisticated hydrodynamic models, with a high physical information demand, these simplified models were the only possibility in estimating water levels and discharges in the near future. Many of these models have a simple model structure and a limited number of parameters, which makes their calibration and simulation very simple and hence computationally efficient. These two characteristics make them also very applicable in the conceptual modelling approach. This chapter discusses different modelling techniques for river courses and adjacent water bodies that can be used in the conceptual modelling approach. Many of these techniques were developed in the previous century, but still remain applicable. However, many of them have been modified and adapted recently to achieve better results. The techniques can be divided into two main classes: channel routing techniques for river reaches without important artificial influences and reservoir routing for river courses with significant backwater effects.
4.1. Channel flow routing Channel flow routing techniques basically transform an upstream inflow hydrograph into a downstream outflow hydrograph. Most of these techniques use a simple linear regression model in which the present outflow value is a function of in- and outflow values at previous and present time steps. Some multilinear and non-linear methods are available to account for the different responses of rivers at different states (water levels, discharges, …) This section describes a number of simple channel routing techniques that have been used in the past and whereof some are still being used at the moment. It will be shown that most of these models and methods can be generalized into one modelling technique that provides a powerful method to model the propagation of flood waves through rivers. 4.1.1.
Convex routing
The convex routing method is one of the oldest routing methods used in engineering practice. It was implemented in the 1965 version of the WinTR-20 software, a runoff and routing package developed by the National Resources Conservation Service of the US Department of Agriculture (NRCS, 1965). The method is based on the principle that when a natural flood flow passes through a natural stream channel, without significant local inflows or losses, a reach length L and a time interval Δt exist, such that the value of the downstream discharge Q(t) lies in the interval delineated by the up- and downstream discharge I(t-Δt) and Q(t-Δt) at the previous time step. (Mockus and Styner, 1972) This principle can be expressed as: If
I (t − ∆t ) ≥ Q(t − ∆t )
then
I (t − ∆t ) ≥ Q(t ) ≥ Q(t − ∆t )
If
I (t − ∆t ) ≤ Q(t − ∆t )
then
I (t − ∆t ) ≤ Q(t ) ≤ Q(t − ∆t )
The routing equation is derived from typical inflow and outflow hydrographs (e.g. in Figure 4-1). For a properly selected reach length L and, hence, a time step Δt, the downstream discharge Q(t) will fall somewhere on or between I(t-Δt) and Q(t-Δt) in magnitude, but not above or below them. The next step is to consider that I(t-Δt), Q(t-Δt) and Q(t) are members of a convex set, which means that if A and B are in the set, then all points on a straight line connecting A and B are also in that set (Bishop and Phelps, 1961). This theory of convex sets can now be applied for the routing equation as:
Q(t ) = (1 − C ) ⋅ Q(t − ∆t ) + C ⋅ I (t − ∆t )
( 4.1 )
Whit C a constant parameter between zero and one.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
21
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 4-1. Derivation of relationships for convex routing method. After (Mockus and Styner, 1972).
It follows from equation ( 4.1 ) that:
C=
Q(t ) − Q(t − ∆t ) I (t − ∆t ) − Q(t − ∆t )
( 4.2 )
The relationships between I(t-Δt), Q(t-Δt) and Q(t) form similar triangles (see upper right corner in Figure 4-1), so that:
Q(t ) − Q(t − ∆t ) I (t − ∆t ) − Q(t − ∆t ) = ∆t K
( 4.3 )
Where K, with units of time, is considered the reach travel time for a selected steady flow discharge of a water particle through the given reach. Combining equations ( 4.2 ) and ( 4.3 ) leads to:
C=
∆t K
( 4.4 )
From which comes the equation that defines Δt, the wave travel time and also the required routing time interval:
∆t = C ⋅ K 4.1.2.
( 4.5 )
Unit hydrograph approach
The technique of the unit hydrograph approach for stream flow prediction was developed in the 1930s. The first notions and proposals can be found in (Sherman, 1932). The unit hydrograph approach translates a time series of upstream discharges (or effective rainfall) into a corresponding time series of downstream outflows. As such, it influences the size of flow peaks and the duration of the flow after a peak. A unit hydrograph recession curve shows the stream flow response resulting from an impulse of one unit of input (Andrews et al., 2010). By using the superposition or convolution principle, a unit hydrograph can be transformed into a hydrograph for any given input. The process of converting an inflow hydrograph into an outflow hydrograph is a mixture of translation and reservoir action (Dooge, 1959). Both cases of pure reservoir action and pure translation are discussed below. The linear reservoir model is based on the continuity equation (Chow et al., 1988):
dS = I −Q dt Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
( 4.6 )
WL2014R00_131_1
22
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
It is assumed that the reservoir outflow discharge is directly proportional to the storage in the reservoir as:
Q = S /k
or
S = k ⋅Q
( 4.7 )
The parameter k (with units of time) is the equivalent of the mean residence time of the reservoir and is called the reservoir constant. After combining equations ( 4.6 ) and ( 4.7 ) the differential equation of the linear reservoir can be obtained:
k⋅
dQ = I −Q dt
( 4.8 )
Equation ( 4.8 ) can be solved analytically to obtain the instantaneous unit hydrograph of a single linear reservoir:
Q(t ) =
1 −t ⋅ exp k k
( 4.9 )
Solution of equation ( 4.8 ) with a simple explicit finite difference form of the continuity equation yields the outflow as:
−1 −1 Q(t ) = exp ⋅ Q(t − ∆t ) + 1 − exp ⋅ I (t ) k k
( 4.10 )
Where Q(t) and Q(t-Δt) are the reservoir outflows at time steps t and t-Δt, I(t) is the inflow at time t, and Δt is the time step or routing interval. Based on the backward shift operator z-n, that shifts the time series with a number (n) of time steps :
z − n ⋅ Q(t ) = Q(t − n ⋅ ∆t )
( 4.11 )
equation ( 4.10 ) can be rearranged and expressed as:
Q(t ) =
b ⋅ I (t ) 1 − a ⋅ z −1
( 4.12 )
with:
−1 a = exp k −1 b = 1 − exp k
( 4.13 )
The linear reservoir model discussed above only gives attenuation of the inflow hydrograph, without accounting for the translation time of the flood wave. In order to address the problem of translation, the concept of a linear channel can be introduced. A linear channel is defined as a reach in which at every point a linear relationship between discharge and area exists. This can be written as (Dooge, 1959):
A(t ) = d ′( x ) ⋅ Q
( 4.14 )
Combining this with the continuity equation:
∂Q ∂A + =0 ∂x ∂t
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
( 4.15 )
WL2014R00_131_1
23
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 4-2. Effects of attenuation by a linear reservoir (red line) and translation by a linear channel (green line) on an inflow hydrograph (blue line).
gives:
∂Q ∂Q + d ' (x ) ⋅ =0 ∂x ∂t
( 4.16 )
Which has the following solution:
Q[t − d ( x )] = C
( 4.17 )
With d(x) the translation time and C a constant. This solution indicates that a linear channel will translate any inflow hydrograph into an outflow hydrograph, without any change of shape, but shifting it in time. The linear reservoir and linear channel component can be combined in series into a model that incorporates translation and reservoir action (see Figure 4-2). The order in which they occur is of no importance, since the translation due to the linear channel merely involves a shifting of the time scale. The outflow of a linear reservoir and linear channel coupled in series is found by combining equations ( 4.9 ) and ( 4.17 ):
−1 −1 Q(t ) = exp ⋅ Q(t − ∆t ) + 1 − exp ⋅ I (t − d ⋅ ∆t ) k k
( 4.18 )
With k and d natural numbers. 4.1.3.
Att-Kin routing model
The convex routing method in the WinTR-20 software was replaced in 1981 by the Att-Kin routing model, which was developed by the Soil Conservation Service of the US Department of Agriculture. The Att-Kin routing model combines the features of a storage and a kinematic routing model. The storage routing component gives attenuation of the routed hydrograph without translation time, while the kinematic routing model incorporates translation but no attenuation. The Att-Kin model (so called because it gives attenuation and kinematic translation) solves both models in a way that satisfies the conservation of mass at the outflow peak (Theurer et al., 1981). In the first step, storage routing of the inflow hydrograph is done using the storage indication method, to solve the integral form of the conservation of mass equation. The equations which reflect attenuation are (SCS, 1986):
I −Q =
S (t 2 ) − S (t1 ) (t 2 − t1 )
( 4.19 )
Q(t ) = k ⋅ [S (t )]
m
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
24
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Where I and Q are the average inflow and outflow over the time interval t2 - t1 , S(t) is the storage in the reach at time t and k and m are the coefficient and exponent of the storage-discharge relation. The storage routed hydrograph is in the second step translated and distorted by kinematic routing. This routing technique solves the differential form of the conservation of mass equation by using an areadischarge relation:
∂Q ∂A + =0 ∂X ∂t
( 4.20 )
Q(t ) = x ⋅ [A(t )]
m
Where Q is the discharge at any distance X and time t, A is the area of the cross section at any distance and time, and x and m are the coefficient and exponent of the area-discharge relation. The exponent m has the same value in the storage-discharge relation ( 4.19 ) as in the area-discharge relation ( 4.20 ). The procedure of the original Att-Kin model is as follows: the inflow hydrograph is routed with assumed values for k, x and m, first with the storage component and secondly with the kinematic component. The routed hydrograph is then checked for conservation of mass between inflow, outflow and storage at the time to peak of the outflow hydrograph. Appropriate corrections are made to the assumed values of the coefficients to assure the conservation of mass. The entire process is repeated iteratively until the conservation of mass at the time to peak is satisfied (Theurer et al., 1981). The modified Att-Kin model (Theurer et al., 1981) is derived in a similar way as the convex routing model, but with one major exception: an arithmetic average is assumed for the outflow in the discrete form of the continuity equation:
I (t − ∆t ) −
Q(t − ∆t ) + Q(t ) S (t ) − S (t − ∆t ) = ∆t 2
( 4.21 )
Assuming that S = K.Q, the above equation can be rewritten as:
I (t − ∆t ) −
Q(t − ∆t ) + Q(t ) K = ⋅ (Q(t ) − Q(t − ∆t )) 2 ∆t
( 4.22 )
which yields:
Q(t ) = C R ⋅ I (t − ∆t ) + (1 − C R ) ⋅ Q(t − ∆t )
( 4.23 )
with
CR =
2 ⋅ ∆t 2 ⋅ K + ∆t
( 4.24 )
In the modified Att-Kin model a direct solution is used, instead of an iterative structure like the original AttKin model. The value of K is derived from the slope of the storage-outflow curve. The best straight line estimate of the slope is the tangent to the curve at the peak inflow as this represents the wave celerity or kinematic wave velocity (SCS, 1986).
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
25
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Slope of line = K
Storage FT³
SS
Discharge CFS
Q pi
Figure 4-3. Storage-outflow slope used in the modified Att-Kin routing method (SCS, 1986).
4.1.4.
Lag-and-route model
The lag-and-route model was initially developed as a graphical routing technique, developed by Linsley et al. (1975). The method consists of two components: a routing component which transforms an upstream flow into a downstream flow and a lag component to account for the travel time of the flood wave. Different model forms are available (e.g. parameter values can be either constant or variable; the relation between storage and discharge can be linear or quadratic; …), which makes the model very flexible (Bentura and Michel, 1997). The linear lag-and-route model is based on the continuity equation and the assumption that the outflow of a reservoir can be related to the storage in the reservoir by a linear relationship: equations ( 4.6 ) and ( 4.7 ). Although the base for this model is the same as that of the linear reservoir model, a different solution is obtained. It is assumed that the upstream discharge is described by straight line segments between data ordinates, instead of constant values during one time step:
t I (t ) = I 1 + (I 2 − I 1 ) ⋅ ∆t
( 4.25 )
With I(t) the inflow at time t, Δt the time step and I1 and I2 the upstream flow rates at the beginning and end of the time step. Substitution of equation ( 4.25 ) in ( 4.8 ) yields:
k⋅
dQ t + Q = I 1 + (I 2 − I 1 ) ⋅ dt ∆t
( 4.26 )
This differential equation is then analytically solved. The solution is the sum of the general solution of the linear equation (with a zero right hand side of equation ( 4.26 )) and a particular solution of the complete equation:
−t Q(t ) = [Q1 − I 1 + (I 2 − I 1 )k / t ] ⋅ exp + (I 2 − I 1 )t / t + I 1 − (I 2 − I 1 )k / t k
( 4.27 )
The outflow at the end of the time step can then be calculated by putting t = Δt. After rearranging and substituting the coefficients C0, C1 and C2, equation ( 4.27 ) can be written as:
Q(t ) = C 0 ⋅ Q(t − ∆t ) + C1 ⋅ I (t − ∆t ) + C 2 ⋅ I (t )
( 4.28 )
with:
− ∆t C 0 = exp k
( 4.29 )
k C1 = ⋅ (1 − C 0 ) − C 0 ∆t
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
26
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
C2 = 1 −
k ⋅ (1 − C 0 ) ∆t
The next step is to consider the lag component. The outflow at the end of the channel reach will be Q(t + T), evolving from Q1 at time T and Q2 at time Δt + T. Hence, the linear lag-and-route model has just two parameters: k and T, both with units of time. In the quadratic lag-and-route model the relation between the outflow of the reservoir and storage in the reservoir is assumed to be quadratic, instead of linear (Bentura and Michel, 1997):
Q = [S / k ]
2
( 4.30 )
or, more generally:
Q = [S / k ]
α
( 4.31 )
Substituting this non-linear relation in the continuity equation ( 4.6 ) yields:
dS S = I − dt k
α
( 4.32 )
This differential equation can be solved analytically for the quadratic case, where α = 2. For higher values of α this becomes impossible. The analytical solution presented in (Bentura and Michel, 1997) is only possible when the upstream inflow is taken to be constant over the time step. Therefore the time step Δt is divided into N subintervals δt, with a constant flow rate I* in each subinterval. The analytical solution for differential equation ( 4.31 ) then becomes:
I * ⋅δt 0.5 Q(t ) + k Q(t + δt ) = 0.5 Q(t ) ⋅ δt + 1 k
2
( 4.33 )
To calculate the downstream outflow at time t + Δt equation( 4.32 ) needs to be iterated N times. As for the linear model, allowance has to be made for a lag T, where T is again the second parameter of the model. 4.1.5.
Kalinin – Milyukov model
The Kalinin – Milyukov model is named after two Russian hydrologists who developed the model in 1958. The idea behind the model is that a river reach can be divided into a number of elementary reaches (with a characteristic length Lc) in such a way that each individual reach can be modeled with a linear reservoir model. The unsteady flow in the river is then approximated by a quasi-stationary flow in every single reservoir (Voll, 2005). During flood conditions, the assumption of a stationary flow is not valid and a hysteresis between outflow and water level exists (see Figure 4-4). Considering this hysteresis curve, it is possible to assign three different water levels to one discharge: one on the rising and falling limb of the non-stationary flow and one for the case of stationary flow. The water levels during the rising and falling limb of the flood, will arise at a certain time step Δt before those of the stationary case. A clear relation between discharge and water level can then be found by using water levels that arise at a time Δt after the discharges. In the same way the water levels in a downstream section can be used to compose a relation with discharges in a certain cross section (e.g. water levels in section r and discharges in section m of Figure 4-5). Extending on this idea the characteristic length can be derived from the stationary (Q, h)-relation as resented in equation ( 4.34 ). A full mathematical derivation of this formula can be found in TUHH (2006)
Lc =
Q dh ⋅ S f dQ(h )
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
( 4.34 )
WL2014R00_131_1
27
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
The characteristic length calculated with equation ( 4.34 ) isn’t constant over the entire reach of the (Q,h)relation. The actual value should be taken as the average of all (or a great number of) values of the characteristic length Lc. The number of linear reservoirs (n) to model a given river reach with length L is easily calculated as:
n=
L Lc
( 4.35 )
The recession constant k of all the linear reservoirs follows from the linear reservoir principle:
k=
dV ∆V ≅ dQ ∆Q
( 4.36 )
Figure 4-4. Hysteresis on rating curve
Figure 4-5. Characteristic length for Kalinin – Milyukov method
The Kalinin – Milyukov routing method has two routing equations: a discrete one for one linear reservoir and a continuous one for a series of linear reservoirs. In discrete form the following equation is used, where the outflow from one reservoir is the inflow of the next one.
Q(t ) = Q(t − ∆t ) + [I (t − ∆t ) − Q(t − ∆t )] ⋅ C1 + [I (t ) − I (t − ∆t )] ⋅ C 2
( 4.37 )
− ∆h C1 = 1 − exp k
( 4.38 )
With:
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
28
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
C2 = 1 −
k ⋅ C1 ∆h
The total reach can be understood as a unit, if the reservoir chain is represented as a summation of the individual reservoirs. This can be expressed by a continuous form of the routing equation, which is equivalent to the Gamma-distribution:
∆t 1 Q(n ⋅ ∆t ) = ⋅ I (n ⋅ ∆t ) ⋅ k ⋅ (n − 1)! k
n −1
⋅ e − n⋅∆t / k
( 4.39 )
Note also that the concept of the Kalinin – Milyukov model, representing a river reach by a series of linear reservoirs, is analogous to the Nash-Cascade model (Nash, 1960). Both models were developed more or less simultaneous, but independently. 4.1.6.
Muskingum routing method
The Muskingum method of channel routing is a well-known and regularly used procedure. The credit for developing the method is given to McCarthy (1938) of the US Corps of Engineers. The method is based on the assumption that the storage in a river reach consists of two parts: prism storage and wedge storage. Prism storage is essentially the storage under the steady flow water surface profile, while wedge storage is the additional storage under the actual water surface profile. As it is shown in Figure 4-6, the wedge storage is positive during the rising stages of the flood wave and added to the prism storage. During the falling stages of a flood wave, the wedge storage is negative and subtracted from the prism storage (USACE, 1994). The principle of dividing the total storage in a prism storage and a wedge storage can be expressed as:
S = K ⋅ Q + K ⋅ X ⋅ (I − Q ) = K ⋅ [X ⋅ I + (1 − X ) ⋅ Q ]
( 4.40 )
Where S is the storage, I the upstream inflow, Q the downstream outflow, K the travel time of the flood wave through the reach and X a dimensionless weighting factor between 0 and 0.5. The Muskingum routing equation is obtained by combining equation ( 4.33 ) with the continuity equation. Therefore an approximate form of the differential equation ( 4.6 ) is used, by integrating over the routing interval Δt :
Q(t − ∆t ) + Q(t ) I (t − ∆t ) + I (t ) S (t ) − S (t − ∆t ) − = ∆t 2 2
( 4.41 )
Combining equations ( 4.40 ) and ( 4.41 ) yields:
Q(t ) = C1 ⋅ I (t − ∆t ) + C 2 ⋅ I (t ) + C 3 ⋅ Q(t − ∆t )
( 4.42 )
∆t + 2 ⋅ K ⋅ X ∆t + 2 ⋅ K ⋅ (1 − X ) ∆t − 2 ⋅ K ⋅ X C2 = ∆t + 2 ⋅ K ⋅ (1 − X ) 2 ⋅ K ⋅ (1 − X ) − ∆t C3 = ∆t + 2 ⋅ K ⋅ (1 − X )
( 4.43 )
with:
C1 =
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
29
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 4-6. Muskingum prism and wedge storage concept.
The parameters K and X have no real physical meaning and should be calibrated based on measurements or on the results of a detailed hydrodynamic model. However, Cunge (1969) showed that it is possible to estimate the values of K and X from hydraulic properties of the river reach:
X =
Q0 1 ⋅ 1 − B ⋅ S 0 ⋅ c ⋅ ∆x 2
( 4.44 )
∆x K= c
Where Q0 is the discharge under steady flow conditions, B the river width, S0 the bed slope, c the theoretical wave celerity and Δx the distance increment. K can now be interpreted as the travel time of the flood wave. The accuracy of the estimates for X and K with equations ( 4.44 ) is tested in (Ponce et al., 1996), where a good agreement between the numerical solution and the Muskingum-Cunge routing method was found. 4.1.7.
Kinematic wave model
The kinematic wave approximation is the most simplified way of solving the ‘de Saint-Venant’ equations. It is assumed that the local, convective and pressure slopes in the momentum equation can be neglected, which results in the assumption that the friction slope balances the bed slope, so that S0 = Sf. Flows of this nature will not be accelerating appreciably and the flow will remain approximately uniform along the channel. Under these conditions the discharge can be described as a function of water depth (MacArthur and De Vries, 1993; Zoppou, 2001). Smith (1980) and Biesenthal (1974) employed a finite-difference technique to obtain a solution of the continuity equation. This equation can, with reference to the rectangular space-time grid of Figure 4-7, be expressed in finite difference form as:
∆Q ∆A + =0 ∆x ∆t
( 4.45 )
With Q the discharge and A the cross-sectional area of the main channel, at location x and time t. In Figure 4-7, the conditions (discharge and cross-sectional area) at points 1, 2 and 3 are known and will be used to compute the conditions at point number 4. Two weighting factors α and β, respectively in the space- and time-direction, are introduced to define the finite differences between the different points. A quantity ø and the spatial and temporal gradients Δø / Δx and Δø / Δt may then be calculated in terms of α and β and the four values of ø at the corners of the space-time molecule.
φ m = α ⋅ (1 − β ) ⋅ φ1 + α ⋅ β ⋅ φ 3 + (1 − α ) ⋅ (1 − β ) ⋅ φ 2 + (1 − α ) ⋅ β ⋅ φ 4
( 4.46 )
1 ∆φ = ⋅ [(1 − β ) ⋅ φ 2 + β ⋅ φ 4 − (1 − β ) ⋅ φ1 + β ⋅ φ 3 ] ∆x ∆x
( 4.47 )
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
30
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
∆φ 1 = ⋅ [α ⋅ φ 3 + (1 − α ) ⋅ φ 4 − α ⋅ φ1 − (1 − α ) ⋅ φ 2 ] ∆t ∆t Applying this to the terms of continuity equation ( 4.45 ) yields:
∆Q 1 = ⋅ [β ⋅ (Q4 − Q3 ) + (1 − β ) ⋅ (Q2 − Q1 )] ∆x ∆x ∆A 1 ⋅ [α ⋅ ( A3 − A1 ) + (1 − α ) ⋅ ( A4 − A2 )] = ∆t ∆t
( 4.48 )
Substitution of these expressions in the continuity equation gives a finite-difference approximation of the continuity equation centered around a point P in the space-time grid, also called nucleus. This point may lie anywhere within or on the bounding edges of the rectangular molecule defined by the space and time increments Δx and Δt, depending on the values of α and β, which both lie in the interval between zero and one. Re-arrangement of the resulting equation leads to:
∆x ∆x = β ⋅ Q2 + (1 − α ) ⋅ A2 ⋅ − Q2 − ∆t ∆t ∆x ∆x + Q1 + β ⋅ Q3 − α ⋅ A3 ⋅ β ⋅ Q1 + α ⋅ A1 ⋅ ∆t ∆t
β ⋅ Q4 + (1 − α ) ⋅ A4 ⋅
( 4.49 )
The two unknowns Q4 and A4 are grouped in the left-hand side of equation ( 4.49 ). Because of the kinematic wave approximation a solution can be obtained by assuming a unique relation between the discharge and the cross-sectional area. Equation ( 4.49 ) can be rewritten in terms of the following functions:
f (Q ) = β ⋅ Q − α ⋅ A
∆x ∆t
( 4.50 )
∆x g (Q ) = β ⋅ Q + (1 − α ) ⋅ A ∆t as:
g (Q4 ) = g (Q2 ) + f (Q3 ) − f (Q1 ) + Q1 − Q2
( 4.51 )
Knowledge of the stage-discharge curve at the upstream and downstream sections allows the f(Q) and g(Q) functions to be computed. For known discharge values at points 1, 2 and 3 the value of g(Q4) is obtained and thus Q4 by back-interpolation. For stability reasons the values of Δx and Δt are limited, which might restrict the use of the kinematic wave model as a conceptual model. Smith (1980) shows that the condition for stability can be expressed as:
α 1−α ≤ Cr ≤ β 1− β
( 4.52 )
With Cr the Courant stability number:
Cr =
∆t ⋅ c ∆x
( 4.53 )
where c is the wave velocity. For example when α and β are equal to 0.5 or when the point P lies in the center of the molecule in Figure 4-7, Cr must be equal to one.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
31
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 4-7. Definition sketch of molecule in a space-time grid (Smith, 1980).
4.1.8.
Unified coefficient routing model
Over the years many different simplified routing models have been developed, whereof the most widely used are listed above. Fread (1983) shows that many of these models (convex routing, unit hydrograph approach, lag-and-route, Kalinin-Milyukov, Muskingum and kinematic wave routing) can be unified into one single model. This unified coefficient routing model has the following form:
Q(t ) = C1 ⋅ I (t − ∆t ) + C 2 ⋅ I (t ) + C 3 ⋅ Q(t − ∆t ) + C 4
( 4.54 )
In which C1, C2 and C3 are routing coefficients which may be derived empirically or evaluated from the hydraulic characteristics of the channel reach. The coefficient C4 accounts for the effect of lateral inflows along the routing reach. This last term could also be neglected by adding the lateral inflows to the upstream discharges I. The unified coefficient routing model can be used in two different ways: as a linear model or as a nonlinear model. In the linear case the coefficients are considered to be constant for each routing reach and throughout the duration of the routing computation. However, in the nonlinear case of the unified coefficient routing model the coefficients can vary with each routing reach and with time. The value of the coefficients then depends on the average discharge or water depth in the previous routing interval (Fread, 1983). 4.1.9.
Transfer functions
A more general formulation of the linear reservoir model ( 4.8 ) is based on the following linear function (Chow et al., 1988):
d 2Q d n −1Q dQ + a 3 ⋅ 2 + ... + a n ⋅ n −1 dt dt dt dI d 2I d m −1 I + b1 ⋅ I + b2 ⋅ + b3 ⋅ 2 + ... + bm ⋅ m −1 dt dt dt
S = a1 ⋅ Q + a 2 ⋅
( 4.55 )
where a1,..,an and b1,..,bm are time invariant coefficients. Finding the derivative of the reservoir storage S with respect to time and replacing them into the continuity equation ( 4.6 ) yields a continuous linear system expressed in differential operator form:
dn d n −1 d2 d a n ⋅ n + a n −1 ⋅ n −1 + ... + a 2 ⋅ 2 + a1 ⋅ + 1 ⋅ Q = dt dt dt dt m −1 m 2 d d d d 1 − b1 ⋅ − b2 ⋅ 2 − ... − bm −1 ⋅ m −1 − bm ⋅ m ⋅ I dt dt dt dt Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
( 4.56 )
32
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
d d2 d m −1 dm 1 − b1 ⋅ − b2 ⋅ 2 − ... − bm −1 ⋅ m −1 − bm ⋅ m dt dt dt dt M (D ) Q= ⋅I = ⋅I n n −1 2 N (D ) d d d d a n ⋅ n + a n −1 ⋅ n −1 + ... + a 2 ⋅ 2 + a1 ⋅ + 1 dt dt dt dt
( 4.57 )
The combined operator M(D)/N(D) is called the system’s transfer function and is used to relate the system inputs with the outputs. Comparison between equations ( 4.13 ) and ( 4.56 ) shows that the linear reservoir can be regarded a first order transfer function. The discrete form of the transfer function given by ( 4.56 ) expressed in terms of the backward shift operator z-1 is:
Q(t ) =
b0 + b1 ⋅ z −1 + ... + bm ⋅ z − m 1 − a1 ⋅ z
−1
− ... − a n ⋅ z
−n
⋅ z −δ I (t )
( 4.58 )
Figure 4-8. Arrangements of linear reservoirs that represent the transfer function structures: a) [2 3 0] and b) [3 1 0] (Villazon et al., 2011).
Or: n
m
i =1
j =0
Q(t ) = ∑ ai ⋅ Q(t − i ⋅ ∆t ) + ∑ b j ⋅ I (t − ( j + δ ) ⋅ ∆t )
( 4.59 )
where m and n are called the transfer function ‘orders’ and δ is the time delay for the input series. One of the most important characteristics of the transfer function models is that it is possible to represent them as an arrangement of linear reservoirs. For example, the transfer function of two reservoirs in series can be obtained as the product of the individual linear reservoirs; for two reservoirs in parallel the two first order components have to be added, which yields a second order transfer function (Willems, 2000). Comparison of equations ( 4.54 ) and ( 4.59 ) proves that the unified coefficient routing model can be expressed as a transfer function of first order, with n = 1 and m = 1. Hence, most of the conceptual routing models discussed above can be considered as some sort of transfer function, with a predefined structure. 4.1.10. Non-linearity It is an accepted fact that the propagation of a flood wave is a non-linear process (e.g. because of the inundation of floodplains) and therefore application of the previously described linear models with timeinvariant parameters may result in crude approximations of actual waves. (Camacho and Lees, 1999). It was pointed out that if the model is calibrated around high discharges the low flows arrive to soon and are over-damped. On the other hand, if the model is calibrated around low discharges, the peaks will arrive too
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
33
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
late and are under-damped (Perumal, 1992). Several techniques are proposed to overcome this problem, whereof the most widely used are mentioned here. The first possibility in dealing with the non-linear behaviour of a flood wave is the use of non-linear reservoirs, instead of linear ones. The outflow of such non-linear reservoirs is related to the storage of the reservoir by a power relation (Hughes and Murrell, 1986):
S = k ⋅ Qn
( 4.60 )
with k > 0 and 0 < n < 1 model parameters, or by a logarithmic relation:
S = k ⋅ ln(Q )
( 4.61 )
In the case of the power relation no analytical solution to this problem exists unless n is equal to ½ or 1. Hughes and Murrell (1986) present some numerical solutions to the system of equations ( 4.6 ) and ( 4.60 ), based on first-order approximations of one or two of these equations. Because of the iterative nature of these solutions the computational effort can become quite large, which restricts the use of these techniques in conceptual channel routing models, where a straight-forward solution is desirable.
Figure 4-9. Computational procedure of a multilinear model. After (Kim, 2011).
The second possibility is the use of so called multilinear or multiple input linear models, which describe a nonlinear system with a set of linear models (see Figure 4-9). The principle of this technique is to distinguish different components on the input hydrograph, each of which is subsequently routed through a linear model (Perumal, 1992). Combining the responses Qj(t) of each linear sub-model to the imposed inflow hydrographs Ij(t), gives the overall outflow hydrograph Q(t). The multilinear principle has previously been applied to the Muskingum method (Perumal, 1992), the Kalinin-Miljukov model (Szolgay et al, 2006) and the unit hydrograph approach (e.g. Perumal, 1994; Camacho and Lees, 1999). Note that the principle of the multilinear models is also included in the transfer function methodology, since these models can be represented as an arrangement (in series or in parallel) of linear reservoirs (Willems, 2000; Beven, 2001). A multilinear transfer function model consists then of a number of linear reservoirs in parallel. The next option is a non-linear transformation of the input, to account for the differing responses at different river states. Previous studies (Young, 2001; Romanowicz et al., 2008; Beven et al., 2009) show that such a transformation is able to successfully represent the non-linear behaviour of rainfall-runoff or the flood effect along rivers. The principle of the input transformation can be expressed as (Villazon, 2011):
I ′j (t ) = f j [I (t )] ⋅ I (t )
( 4.62 )
Where I’j(t) is the transformed input at location j and time t, which is the product of the input I(t) and the nonlinear transformation function fj[I(t)], whose value also depends on the input state. The transform I’j(t) is used as the input to the linear transfer function model. This technique was successfully applied by Romanowicz et al. (2008) to the river Severn case (UK) and by Villazon (2011) to the river Dender case (Belgium). Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
34
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Storage S
Linear Reservoir
Non-linear Reservoir
k1 k2 Discharge Q
Figure 4-10. Storage-discharge relations: linear versus non-linear reservoirs.
The last possibility is to consider that the nonlinear model can be represented as a small perturbation of linearized models and that the linearized model parameters can be related to state variables. This technique is called the state parameter dependency method (Young, 1998) and allows the parameters a1,..,an and b0,..,bm of the linear transfer function ( 4.58 ) to vary, according to the current upstream or downstream flow. In the case where only the first parameter (b0) is estimated as a state dependent parameter with the remaining parameters fixed, the state parameter dependency method simplifies to the non-linear input form. The technique has been applied to rainfall-runoff modelling and to flood routing in the river Trent catchment (Lees, 2000; Lekkas et al., 2001). A link can be made to the non-linear reservoirs of equation ( 4.60 ) and the concept of relating linearized model parameters with state variables. For a linear reservoir, it is assumed that outflow and storage are directly proportional, by a constant k, also called the reservoir constant. This constant is also the only parameter of equation ( 4.10 ), which allows to calculate the outflow of the linear reservoir. In the case of the non-linear reservoir a linear relation can be established between any combination of storage and discharge on the non-linear relation and the origin. See for example k1 and k2 in Figure 4-10. Because of the non-linear character of equation ( 4.60 ), the value of the reservoir constant will not be constant, but will depend on the discharge or the storage. Linear transfer functions with state dependent parameters thus still have a physical basis, which justifies their use in conceptual models. 4.1.11. Conclusion Channel routing techniques transform an upstream inflow hydrograph into a downstream outflow hydrograph. Several simplified channel routing methods, suitable for the conceptual modelling approach, were discussed. Most of these methods start from the conservation of mass equation ( 4.6 ) or the continuity equation ( 4.15 ) and use a certain relation between discharge and storage to calculate a downstream discharge. This results in a linear regression model in which the present outflow value is a function of in- and outflow values at previous and present time steps. It was shown that these methods can be considered as special cases of a general model structure: the unified coefficient routing model, which covers a wide range of possible model structures. This unified coefficient routing model is on his turn a specific case of a more general method that uses transfer functions to calculate a downstream discharge based on in- and outflow values at one or more previous time steps. Furhtermore, this model structure can, be seen as an arrangement of linear reservoirs, which provides a physical basis to the model. Drawback of linear regression models is the difficulty to describe the non-linear behaviour of flood propagation in natural channels. A number of methods to account for these non-linearities were described and it can be concluded that the transfer function model, extended with state dependent parameters, provides the most general model structure and is therefore probably the most powerful. Using this technique will probably lead to a methodology capable of modelling almost all channel routing problems. It is important to note that the use of channel routing techniques is restricted to river courses without artificial influences like regulation structures. These artificial influences will cause significant backwater Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
35
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
effects, which the channel routing techniques can not handle explicitly. In all situations where (important) backwater effects arise reservoir routing techniques should be used.
4.2. Reservoir routing 4.2.1.
Introduction
Unlike channel routing techniques, reservoir routing methods can handle the backwater effects encountered in different parts of river courses and their adjacent water bodies. These techniques use numerical solutions for the storage equation that governs the rate of change of reservoir storage volume as a function of inflow and outflow rates (Fenton, 1992):
dS = I (t ) − Q(t , S ) dt
( 4.63 )
With S the volume of water stored in a reservoir, t the time, I the inflow discharge and Q the outflow. To solve equation ( 4.63 ) it is necessary to relate the outflow to the storage volume, usually via the dependence of both with the elevation of the water surface. In addition to the storage formulation, other forms of the differential equation ( 4.63 ) can be obtained by considering that storage and water level can be related as:
dS = A(h ) ⋅ dh
( 4.64 )
where A(h) is the plan area of the water surface at elevation h. This yields an equivalent form of the storage equation:
dh I (t ) − Q(t , h ) = dt A(h )
( 4.65 )
This is a differential equation of the surface elevation itself. To solve this differential equation, a relation between the elevation of the water surface and the outflow is again needed. The advantage of the formulation of ( 4.65 ) over that of ( 4.63 ) is that it makes no use of the storage volume S, which then does not have to be calculated for routing purposes. Drawback is that it introduces an extra nonlinearity (or unknown) in the storage equation, since the plan area of the water surface is also dependent on the elevation of the water level. Both forms of the storage equation can be generalized by the expression (Fenton, 1992):
dy = F (t , y ) dt
( 4.66 )
where y is the dependent variable, either S or h. In the case of equation ( 4.63 ), y = S and F(t,S) = I(t) – Q(t,S), while for equation ( 4.65 ), y = h and F(t,h) = [I(t) – Q(t,h)]/A(h). Numerical solutions of the storage equations can be grouped into two classes: the first uses a numerical approximation of the storage equation with the help of the trapezoidal rule, whereas the second solution scheme uses explicit numerical solutions by employing a decomposed form of the storage equation. The former is more accurate but slower, while the latter is straight-forward and less time consuming but possibly less accurate. Different solution schemes for both classes will be discussed below. 4.2.2.
Numerical approximation by using the trapezoidal rule
The traditional method of solving equation ( 4.66 ), described in almost all books on hydrology (e.g. Mockus and Styner, 1972; Chow et al., 1988 and USACE, 1994), uses the following approximation of the differential equation:
y (t ) − y (t − ∆t ) 1 = ⋅ [F (t , y (t )) + F (t − ∆t , y (t − ∆t ))] ∆t 2 Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
( 4.67 )
36
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Which translates equation ( 4.63 ) into:
S (t ) − S (t − ∆t ) I (t ) + I (t − ∆t ) Q(t ) + Q(t − ∆t ) = − ∆t 2 2
( 4.68 )
Equation ( 4.68 ) is in fact a discretization of equation ( 4.63 ) by using the average values of inflow and outflow over the routing interval Δt, which is a finite step in time. The equation contains four known values: the inflow hydrograph, represented by I(t) and I(t-Δt),and the storage S(t-Δt) and outflow Q(t-Δt) at the start of the routing interval. The two unknown values are the storage S(t) and the outflow Q(t) at the end of the routing interval. Outflow and storage are assumed to be related in a nonlinear way via the water level, so that solving equation ( 4.68 ) is equivalent to finding the root of a single nonlinear equation. Many solution schemes (iterative or not) exist for this non-linear problem. A number of them, which have previously been applied to reservoir routing are described in this section. All of these solution schemes are said to be unconditionally stable, which explains the success and continued use of the traditional methods, despite their complexity (Fenton, 1992). 4.2.2.1
Iterative methods
Numerous iterative methods are available to find the root of a single nonlinear equation. Chapra (2012) describes a few of these methods, like simple fixed-point iteration, the Newton-Raphson or tangent method, the secant method and Brent’s method. All of them, except the last, have been used before to solve equation ( 4.68 ) and are described below together with an extra method: the Laurenson-Pilgrim method. The fixed point iteration scheme (also called one-point iteration) presented in (Mockus and Styner, 1972) contains four calculation steps. Equation ( 4.67 ) is rearranged so that y(t) is on the left hand side of the equation (Chapra, 2012):
y (t ) = f ( y ) = y (t − ∆t ) +
∆t ⋅ [F (t , y (t )) + F (t − ∆t , y (t − ∆t ))] 2
( 4.69 )
The advantage of this formulation is that it provides a formula to predict a new value of y as a function of an old value of y. Thus, given an initial guess for the root, yi, equation ( 4.71 ) can be used to compute a new estimate, yi+1:
yi +1 = f ( yi )
( 4.70 )
With i an iteration counter. In the case of the method described in (Mockus and Styner, 1972), the solution starts with assuming a value for the outflow Q(t) at the end of the routing interval. Based on this assumption a value for the reservoir storage S(t) can be computed with equation ( 4.68 ). The reservoir storage is subsequently transformed into a water level using an elevation-storage curve, which on his turn can be used to compute the instantaneous outflow Q’(t). In the last step the calculated and assumed discharges are compared. The routing operation is completed when the difference between Q(t) and Q’(t) lies within the specified degree of accuracy. If the two values do not agree well enough the procedure is repeated once more. It can be shown (Chapra, 2012; example 6.1) that fixed-point iteration has a linear convergence. This means that the relative error for each iteration is roughly proportional to the error from the previous iteration, resulting in a convergence that is rather slow. The second iterative solution scheme, using the Newton-Raphson method, is described in (Fread and Hsu, 1993). Equation ( 4.67 ) can be rewritten as:
I (t ) + I (t − ∆t ) Q (t ) + Q(t − ∆t ) [A(t − ∆t ) + A(t )] ⋅ [h(t ) − h(t − ∆t )] − − =0 2 2 2 ⋅ ∆t
( 4.71 )
The unknowns in equation ( 4.71 ) consist of h(t), A(t) and Q(t). However, A(t) and Q(t) are both known nonlinear functions of h(t), so in fact h(t) is the only unknown. Equation ( 4.71 ) can be solved for h(t) by the iterative Newton-Raphson method, in which:
hi +1 (t ) = hi (t ) −
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
f [hi (t )] f ′[hi (t )]
( 4.72 )
WL2014R00_131_1
37
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Where f(hi) is the left-hand side of equation ( 4.71 ) evaluated with the first estimate for hi(t), which for i = 1 is either h(t-Δt) or a linear extrapolated estimate of h(t), and f’(hi) is the derivative of equation ( 4.71 ) with respect to h(t). In most cases of reservoir routing an analytical derivative of the left-hand side of equation ( 4.71 ) is very hard or even impossible to calculate. An alternative approach is to calculate a numerical derivative by using a fractional perturbation of the independent variable (Chapra, 2012):
f ′[hi (t )] ≅
f [hi (t ) + δhi ] − f [hi (t )] δhi
( 4.73 )
With δhi a small perturbation factor (e.g. 1 cm). This approximation can then be substituted in equation ( 4.72 ) to yield the following iterative equation:
hi +1 (t ) = hi (t ) −
δhi ⋅ f [hi (t )] f [hi (t ) + δhi ] − f [hi (t )]
( 4.74 )
This is called the secant method (Kožar et al., 2010). Unlike the fixed-point iteration scheme, which has a linear convergence, the Newton-Raphson method and the secant method have a quadratic convergence, resulting in an error that is roughly proportional to the square root of the previous error (Chapra, 2012). According to Fread and Hsu (1993) one or two iterations usually solve equation ( 4.71 ) for h(t). Once h(t) is obtained A(t) and Q(t) can easily be computed with the according nonlinear relations. An alternative method, that uses the regula-falsi (or false position) method, was introduced by Laurenson and Pilgrim. The LP-method is represented by the general equation (Fiorentini and Orlandini, 2013):
Yi +1 (t ) = Yi +1 / 2 (t ) − Y (Yi +1 / 2 (t )) ⋅
Yi +1 / 2 (t ) − Yi (t ) Y (Yi +1 / 2 (t )) − Y (Yi (t ))
( 4.75 )
Where Y can either be the outflow discharge Q or the storage S, i is an iteration counter and Ψ(Y) is a function of Q incorporating the direct S-Q function S = F(Q), as given by:
∆t ∆t + F (Qi (t )) + Q(t − ∆t ) ⋅ 2 2 ∆t − (I (t − ∆t ) + I (t )) ⋅ − F (Q(t − ∆t )) 2
Y (Yi (t )) = Qi (t ) ⋅
( 4.76 )
or a function of S incorporating the inverse S-Q function Q = F-1(S) , as given by:
∆t ∆t + S (t ) − F −1 (S i (t − ∆t )) ⋅ 2 2 ∆t − (I (t − ∆t ) + I (t )) ⋅ − S (t − ∆t ) 2
Y (Yi (t )) = F −1 (S i (t )) ⋅
( 4.77 )
The function Ψ(Y) is clearly defined on the basis of equation ( 4.68 ). When this equation is verified, then (and only then) is Ψ(Y) equal to zero. Equation ( 4.75 ) provides an approximation of zero of the function Ψ(Y) obtained by considering the secant for the points (Yi(t), Ψ(Yi(t)) and (Yi+1/2(t), Ψ(Yi+1/2(t)) and the approximate condition Ψ(Yi+1(t)) = 0. At each time interval the LP method starts the iterations with equations ( 4.75 ) and ( 4.77 ) by using S(t- Δt) and S(t- Δt) + [I(t- Δt) – Q(t- Δt)].Δt as initial guesses for Si(t) and Si+1/2(t) respectively. Further details on this method can be found in (Pilgrim, 1987). 4.2.2.2
Non-iterative methods
To avoid the previously described iterative methods, which might be too time consuming, two non-iterative methods are suggested: the storage indication method and a graphical or tabular method. The storage-indication method (Mockus and Styner, 1972), also called modified pulse reservoir routing (USACE, 1994) or level-pool routing method (Chow, 1988) is the most widely used solution scheme, because of its simple approach and accurate solution. Equation ( 4.68 ) can be rearranged to get both unknown variables on the left-hand side and all the known values on the right-hand side: Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
38
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
2 ⋅ S (t ) 2 ⋅ S (t − ∆t ) + Q(t ) = + I (t ) + I (t − ∆t ) − Q(t − ∆t ) ∆t ∆t
( 4.78 )
The right-hand side of equation ( 4.78 ) can easily be calculated, since all of these values are known. A “storage-indication” curve, which relates the outflow Q(t) with the storage indication (2.S(t)/Δt + Q(t)), is then used to transform the left-hand side into an outflow value. Based on the outflow value the corresponding storage in the reservoir can be derived. The major advantages of this method are the straightforward solution scheme and the fact that it only involves a calibration of the storage-indication curve. However, this storage-indication curve is only valid for uncontrolled reservoirs, with fixed outlet levels, and where the outflow of the reservoir is solely determined by the water level in the reservoir. When the reservoir is bounded by weirs with adjustable gates and/or the weir can operate in drowned conditions, the outflow is determined by a nonlinear combination of the upstream and downstream water level and the gate level. This implies that a storage-indication curve is needed for every possible combination of these three elevations, which will result in a three dimensional matrix with a massive amount of calibration data. This problem might be overcome by a next solution scheme that solves the problem in a graphical way or by interpolating in a table of value pairs. The outflow of the reservoir can be expressed as:
Q(t ) = f [hUS (t ), hDS (t ), h0 (t ),...]
( 4.79 )
Where hUS is the upstream water level and a function of the storage in the reservoir, hDS is the downstream water level and h0 is the gate level of the downstream structure. Equation ( 4.68 ) can be rearranged in a similar way so that the outflow becomes a function of the storage:
Q(t ) = I (t ) + I (t − ∆t ) − Q(t − ∆t ) −
2 ⋅ [S (t ) − S (t − ∆t )] ∆t
( 4.80 )
The solution of the routing problem is situated at the intersection of equations ( 4.79 ) and ( 4.80 ) in a (S,Q)-plane. This intersection can be found in a graphical or tabular way, using an optimal numerical technique. Assuming that the outflow and storage will not change very much during the routing interval, the number of (Q,S)-points to calculate and hence the computation time can be limited. 4.2.3.
Explicit numerical solutions
Unlike the methods mentioned in the previous section, the explicit numerical solutions solely use storage and discharge values of the previous time step to calculate storage and outflow at the present time step. These numerical solutions for equation ( 4.66 ) can mathematically be expressed as (Chapra, 2012):
y (t ) = y (t − ∆t ) + ϕ ⋅ ∆t
( 4.81 )
With Δt the time interval and φ the slope, also called increment function. Thus, the slope estimate of φ is used to extrapolate from an old value of y at time step t-Δt to a new value at time step t. Two main classes for estimating the increment function can be distinguished: one-step methods where the value of the slope is based on information at a single point, and multi-step methods that use information from several previous points to estimate the value of the slope (Chapra, 2012). So far, only the one-step methods seem to be applied to reservoir routing. Euler’s method, the simplest but least accurate of all these models, being of first order accuracy only, is given by (Fenton, 1992):
( )
y (t ) = y (t − ∆t ) + ∆t ⋅ F [t − ∆t , y (t − ∆t )] + O ∆t 2
( 4.82 )
With the Landau order symbol O(Δt2) meaning that neglected terms vary like the second or higher powers of the time step Δt. Applying equation ( 4.82 ) to ( 4.63 ) yields:
( )
S (t ) = S (t − ∆t ) + ∆t ⋅ [I (t − ∆t ) − Q(t , S (t − ∆t ))] + O ∆t 2
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
( 4.83 )
39
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
The evaluation on the right-hand side of equation ( 4.83 ) to obtain S(t) is a much simpler process than the numerical solution of the transcendental equation ( 4.68 ) and has to be evaluated only once per time step. The advantage of Euler’s method is its simplicity, but at the mean time this is also the main drawback. The accuracy of the method can be rather poor, which can result in an unstable conceptual model. To overcome this instability problem and to achieve a more accurate model, higher order decomposition models can be used. For example the Runge-Kutta second order method, which is expressed as (Fenton, 1992):
y (t ) = y (t − ∆t ) +
k1 + k 2 + O ∆t 3 2
( )
( 4.84 )
Where:
k1 = ∆t ⋅ F (t − ∆t , y (t − ∆t )) k 2 = ∆t ⋅ F (t , y (t − ∆t ) + k1 )
( 4.85 )
For equation ( 4.63 ) this becomes:
S (t ) = S (t − ∆t ) +
∆t ⋅ [I (t − ∆t ) − Q(S (t − ∆t ))] + 2
∆t ⋅ [I (t ) − Q(S (t − ∆t ) + ∆t ⋅ [I (t − ∆t ) − Q(t − ∆t )])] + O ∆3 2
( )
( 4.86 )
This solution scheme can be explained as consisting of three steps. In the first step the outflow value is approximated with equation ( 4.83 ) and the relation between storage and discharge. This approximation Q’(t) is then used to calculate the storage value S(t) with equation ( 4.68 ), which in this case can be rewritten as:
S (t ) = S (t − ∆t ) +
∆t ⋅ [I (t − ∆t ) − Q(t − ∆t ) + I (t ) − Q ′(t )] 2
( 4.87 )
The storage S(t) can subsequently be transformed into the actual discharge Q(t) with the same storage – discharge relation. A comparison of the calculation schemes of Euler’s method and the Runge-Kutta second order method is depicted in Figure 4-11. Euler’s method uses the slope (equal to the first derivative) at the current time step to predict a new value at the next timestep, by linear extrapolating over the time step. The second order Runge-Kutta is two-step process: Euler’s method is used to obtain a first prediction (predictor step). This point provides an estimate that allows the determination of the slope at the end of the time interval. The two slopes are then averaged to obtain an average slope for the interval, which is again used to extrapolate linearly to a new point (corrector step) (Chapra, 2012). The explicit numerical solution methods are not restricted to first and second order solution. Higher order Runge-Kutta methods, like the third order method described in (Chow et al., 1988) or the fourth order solution described in (Li and Wang, 2007) are also possible. These methods obtain the solution at time t in terms of the solution at t - Δt and the evaluation of F(t,x) at those two times and at intermediate points (Fenton, 1992; Fiorentini and Orlandini, 2013).
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
40
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 4-11. Schematic overview of Euler’s method (a) and Runge-Kutta second order method (b: predictor and c: corrector). (Chapra, 2012)
4.2.4.
Discussion
Several methods for reservoir routing, with different complexity and accuracy, were discussed in the previous two sections. The most favorable of these methods should have an ideal balance between accuracy (and thus avoiding the growth of instabilities) and computational effort, which is defined by the complexity of the method and the time step for routing. For example: Euler’s method is very simple and straightforward, but might need a small time step for stability reasons, whereas the iterative NewtonRaphson technique is more complicated and accurate, which might increase the time step or routing interval. A comparison of these two requirements was previously carried out by Fenton (1992) and Fiorentini and Orlandini (2013).
Figure 4-12. Results obtained by the selected computation schemes (left) and comparison of the accuracy of various schemes and the dependency on the time step used (right) (N is the number of time steps over a total time of 6000 s.) (Fenton, 1992).
Fenton (1992) made a comparison between Euler’s method, the Runge-Kutta second order method, the traditional storage-indication method and the analytical solution of a problem described in table 3 of (Yevdjevich, 1959). This is a very simple problem that describes a detention reservoir with a single sharp crested weir at the outlet. The inflow hydrograph is described by an exponential function, similar to the inflow hydrograph in Figure 4-2. The results of the comparison are shown in Figure 4-12, where the results of the four methods are plotted versus time and the mean error over 20 time steps is plotted versus the number of time steps. It is clear that Euler’s method is the least accurate, since it shows the biggest deviation of the analytical solution, while the difference of the other methods with the analytical solution is hardly visible. Euler’s method is behaving as a first order method, such that halving of the time step halves the error. The error of the second order Runge-Kutta method and the traditional method vary like 1/N², clearly resulting in a more accurate solution. The right-hand side of Figure 4-12 suggests that for small time steps the accuracy of the Runge-Kutta second order method is quite similar to that of the traditional method, while the latter is more time demanding to solve. Fiorentini and Orlandini (2013) compared the robustness, accuracy and efficiency of the Laurenson-Pilgrim method (LP), a fourth order Runge-Kutta method (RK) and two slightly modified versions of the fourth order Runge-Kutta method (CK and RK-B). A real situation, a river flood control reservoir with an observed inflow Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
41
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
hydrograph, was used as case study for their experiments. The smallest CPU-times (a measure for the efficiency of the algorithm) are required by the fourth order Runge-Kutta method for all possible time steps between 1 and 900 s. The CPU-time of the Laurenson-Pilgrim method is always higher, because of the iterative nature of this technique, but the ratio to the CPU-time of the Runge-Kutta method is more or less constant for all time steps. For all methods an inverse proportionality exists between the CPU-time and the time step, or a direct proportionality between CPU-time and the number of time steps. Comparison of the relative error on the peak water level εH again favors the Runge-Kutta method over the Laurenson-Pilgrim method: the former produces the most accurate results as long as the time step is not too large (720 seconds for this case study). For larger time steps the Runge-Kutta method fails to provide a result when the computed points fall outside the domain of the calibrated Q-H relation. This can be solved by using the modified Runge-Kutta method with backstepping (RK-B). Both papers seem to advice the Runge-Kutta method (regardless of the order) as the most favorable algorithm to be used for reservoir routing, with an optimal balance between accuracy and efficiency. Two remarks can be given to this conclusion. The first one is that both papers are incomplete: many more methods for reservoir routing are available and all of them should be evaluated in terms of accuracy and efficiency. The second reason is that both papers use a rather large reservoir with a low gradient of Q with S for their case study, which makes it less sensitive to instability. Small storage reservoirs with high Q-S gradients that can arise in conceptual modelling should be tested as well.
Figure 4-13. Comparison between Laurenson-Pilgrim method (LP) and fourth order Runge-Kutta method (RK) in terms of CPU time (left) and relative error on peak water level (right), for different values of the time step. (Fiorentini and Orlandini, 2013)
4.2.5.
Conclusion
Reservoir routing is the designated method for conceptually modelling water bodies with (important) backwater effects. Examples of this are river reaches bounded by a regulation structure at the downstream end, holding basins and floodplains. All of them have a water level that can be related to the storage in the reservoir. Reservoir routing techniques solve the storage equation ( 4.63 ) with a discrete approximation and a nonlinear relation between the outflow of the reservoir and the storage in the reservoir. Two main categories are available to do so: numerical approximation methods that use a trapezoidal rule and explicit numerical solutions. A number of methods, with different complexity and accuracy, were discussed for both categories. Two papers that compare a couple of possible solution schemes for reservoir routing, were discussed in terms of accuracy and efficiency. It was found however that these papers are incomplete and not completely focused on the problems encountered in conceptual modelling (e.g. the growth of instabilities due to interaction of two adjacent reservoirs). Further research on the different techniques is needed to come to a decision on what technique to use.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
42
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
4.3. Application The concepts of channel routing and reservoir routing were used in two previous studies (Villazon et al., 2012; Meert et al., 2012) to develop two strategies for the construction of conceptual models for river courses and their adjacent water bodies, like floodplains and retention basins. Within these two approaches the river network is schematized by means of ‘building blocks’ that can be connected with each other, to cover a wide range of possible situations. The two strategies are outlined below, together with a guideline on how and where to use them. A first approach is based on the use of transfer functions and is therefore referred to as ‘transfer function methodology’ or ‘TF-method’. The different elements of the method are depicted in Figure 4-14. A transfer function is used for routing the upstream discharge, which might have been subjected to a non-linear transformation, through the segment to achieve a downstream discharge value. Use is made of a (Q,h)relation, also called ‘rating curve’, to calculate the water levels at a point of interest hi (e.g. near floodplains). This can be done from the discharge at the up- or downstream boundary of the river segment, depending on which of these two gives the most univocal relation. The water that flows from the river over the dike or river bank into the floodplain is modelled as a flow over a weir. This flow is sent back to the transfer function, where it is regarded as a negative input value. The floodplain itself is modelled by means of a storage reservoir, with a relation between the volume and the water level in the reservoir, through a so called ‘hypsometric curve’. (Villazon et al., 2012)
Figure 4-14. Transfer function methodology for conceptual models of river courses, consisting of a river reach (in blue) and a floodplain component (in red).
The transfer function methodology is applicable for river reaches without important artificial influences (e.g. flow regulating structures), where the linear reservoir concept is applicable. This approach is no longer valid when the river segment is bounded by a structure whose gate level can change in time. Modelling structures with a constant cross section (e.g. culverts and bridges) is possible with the TF-method, though, since their influence can be considered as invariable. A second approach was set up to allow the incorporation of flow regulating structures in conceptual models. The approach was therefore called ‘hydraulic structures methodology’, or ‘HS-method’. The river reach upstream of the hydraulic structure is modelled with a storage reservoir (based on equation ( 4.83 )), so that the water level at the upstream side of the structure is computed with a hypsometric curve, which is independent of the gate level. The flow over or under the structure is then calculated with the equations that are used in InfoWorks RS (Wallingford, 2012b). These equations need three input values: up- and downstream water level and the gate level. Another hypsometric curve is used to determine the water level at another point of interest in the river reach. The flows into a floodplain and the water level in that floodplain are calculated in the same way as in the transfer function methodology: with the weir equations and a storage reservoir. Figure 4-15 depicts the different components of the hydraulic structures methodology. (Meert et al., 2012) The hydraulic structures methodology is to be used in case of important backwater effects or for storage reservoirs like floodplains and retention basins, where a clear relation between storage and water level exists. This approach is thus advised when dealing with weirs with adjustable gate levels, certainly when they are used in real-time control applications or sensitivity analysis.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
43
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 4-15. Hydraulic structures methodology for conceptual models of river courses, consisting of a river reach (in green) and a floodplain component (in red).
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
44
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
5. Urban drainage systems Urban drainage systems usually consist of the following four major components: an upstream urban catchment, a sewer network, a wastewater treatment plant and the receiving water body. The catchment component estimates the inflow to the sewer system, which contains wastewater from domestic, commercial and industrial uses and/or, depending on the type of sewer network: combined or separated, overland flow caused by precipitation. The sewer network subsequently transports this water to a downstream point where it can be treated in the wastewater treatment plant and discharged into the receiving water. In cases of sewer discharges higher than the capacity of the treatment plant the superfluous amount of water is discharged directly to the receiving water bodies through combined sewer overflows (Saagi, 2013). This chapter will focus on the three first components of the urban drainage system, since the fourth component (the receiving water) was treated in the previous chapter. It is rather obvious that the conceptual models handled in this chapter will show some similarities with the methods and models mentioned in the previous chapter. The outline of this chapter is as follows: The first section gives an overview of the different modelling possibilities for the three major components of the urban wastewater system. In the second section a couple of software packages (commercial and academic) that model the entire system are described. Finally, a link is made to the previous chapter in an attempt to come to a general modelling approach.
5.1. Components of the urban water system 5.1.1.
Urban Catchment
The urban catchment component calculates the amount of water flowing into the sewer network. Two types of inflow can be considered: the dry weather flow caused by the wastewater of domestic and industrial usage and the overland flow generated by a rain event in the urban catchment. Both flows will be discussed separately. 5.1.1.1
Dry weather flow
Dry weather flow is the flow in the system originating from domestic, commercial and industrial wastewater. It plays an important role in the design of combined and separated sewer networks and waste water treatment plants. The main factors that define the amount of dry weather flow are the number of inhabitants living in the catchment and their mean water consumption, the land use in the catchment (domestic, industrial or commercial) and the time of the day (Berlamont, 2004). The two main methods to model this diurnal variation in flow are the use of predefined normalized variations and the fitting of Fourier series to observed data. The most basic form of the former model type is used in the German KOSIM model (ITWH, 2000). The user must specify a mean daily quantity of waste water produced per population equivalent, which is then multiplied by hourly contribution factors to generate a daily dry weather hydrograph. The model provides three standard patterns for domestic (for different population numbers) and one for commercial waste water, but hourly factors can be modified easily for a given case study (Solvi, 2007). In the SEWSYS-model (Ahlman and Svensson, 2005) the dry weather flow consists of two components: a diurnal variation which gives an average value for each hour and a random variation within the hour. The random number lies between 0.7 and 1.2, with a mean value of 1 and is multiplied with the hourly value to become an individual variation within the hour. A more elaborated form of this model type, generating dynamic influent pollutant disturbance scenarios, was developed by Gernaey et al. (2011). The generation of the influent flow rate is achieved by combining the contributions from households, rainfall and infiltration. The influent flow rate caused by households and industry is based on three normalized user-defined data files containing: a diurnal profile, a weekly household flow pattern including the weekend effect and a yearly pattern including the holiday effect. The resulting normalized pattern is then multiplied by the mean flow rate per population equivalent and by the number of population equivalents in the urban catchment. Optionally, zero mean white noise can be added to the flow rate signal. The amount of water infiltrating in the soil is implemented by means of a sine function with higher infiltration rates during the cold season and reduced infiltration during the warm season. The Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
45
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
infiltration is then used as input for a variable volume tank model that describes the soil. When the water level in the groundwater storage tank reaches a certain level, infiltration to the sewer system will occur. The total influent flow rate can be obtained by summing up the contributions of households, industry and infiltration. Various authors have tried to generate dry weather flow rates by fitting Fouries series to observed dry weather flow data (e. g. Langergraber et al., 2008; Rodríguez et al., 2013). The main goal of these studies is mostly to generate realistic influent data and not to provide detailed predictions of the influent of a waste water treatment plant. The dry weather flow patterns are defined by a mean value and a set of hourly multiplier factors, that can be expressed with an nth order Fourier series as (Rodríguez et al., 2013): n
X (t ) = a 0 + ∑ ai ⋅ cos(i ⋅ t ⋅ ω ) + bi ⋅ sin (i ⋅ t ⋅ ω )
( 5.1 )
i =1
Where a0 is a constant (with a value of 1 due to the applied normalization), n is the Fourier series order, ai and bi are constant parameters for each series order (i = 1, 2, …, n) and ω = 2 π / T is the fundamental frequency of the signal, with T = 1 day. Two approaches are now possible to calibrate the parameters of equation ( 5.1 ): assuming that only the total flow can be represented by a Fourier series or disaggregating the total flow in a number of different waste water streams that can each be modeled as a Fourier series (e.g. Langergraber et al., 2008; Rodríguez et al., 2013). The former has less parameters to calibrate and is less complex, while the latter is more accurate and applicable to ungauged catchments. Some models (e.g. De Keyser et al., 2010; Rodríguez et al., 2013) add an extra stochastic component to model sources of variability that cannot be represented entirely by the deterministic Fourier series. 5.1.1.2
Rainfall – runoff generation
In urban drainage systems flows originating from rainfall are of major importance. The amount of surface runoff generated due to a rain event is calculated by a rainfall-runoff model. The level of detail of rainfallrunoff models for urban catchments is in most cases less complicated than those for river catchments, discussed in chapter 3. The calculation process consists of two major steps: determining the net rainfall that will contribute to surface flow and subsequently routing the net rainfall to a point in the watershed where it enters the sewer network. Net rainfall The first step of all rainfall-runoff models is the determination of the net rainfall that contributes to overland runoff. This involves the calculation of losses due to wetting, interception, depression storage, infiltration to the soil layer and evapotranspiration. Various approaches exist to incorporate these losses. Elementary empirical models do not calculate infiltration, evaporation or depression storage, but assume that these values are more or less constant during a rainfall event. In the rational method all rainfall losses are taken into account by a constant runoff-coefficient φ (Berlamont, 2004):
ϕ=
Pnet Ptotal
( 5.2 )
With Pnet the net rainfall and Ptotal the total or measured rainfall. The value of the runoff-coefficient is a function of the land use in the catchment, the permeability and slope of the catchment, the time of the year and many more (Meert and Willems, 2013). The value can be calculated with discharge and rainfall measurements, or roughly estimated from experience or tables, if no measurements are available. Examples of such tables can be found in (Vandekerckhove et al., 2001) and (Meert and Willems, 2013). For catchments with different types of land use an average runoff-coefficient is calculated by a weighted mean over the areas of the different types of land use:
∑ϕ ⋅ A ϕ= ∑A i
i
i
i
( 5.3 )
i
Applications of the rational method can be found in (Mehmood,1995), (Calbró, 2001) and (Gernaey et al., 2011). Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
46
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
A similar approach is used in the SCS-method, developed by the Soil Conservation Service of the U.S. Department of Agriculture (Mockus and Hjemlfelt, 2004). An empirical method is used to determine the maximum potential retention S:
1000 − 10 CN 25400 S= − 254 CN
( 5.4 )
S=
( 5.5 )
Equation ( 5.4 ) is used when S is expressed in inches, while equation ( 5.5 ) expresses S in terms of millimeters. In both equations CN is called the curve number, whose value can be found in tables as a function of soil type and land use. The net rainfall R is then calculated as:
R=
(P − 0.2 ⋅ S )2 (P + 0.8 ⋅ S )
R=0
if
P > 0.2 ⋅ S
if
P ≤ 0.2 ⋅ S
( 5.6 )
With P the rainfall depth, in mm. The derivation of formula ( 5.6 ) can be found in (Mockus and Hjemlfelt, 2004). The SCS-method was developed based on daily rainfall and discharge values, which limits the applicability for lower rainfall events with a (much) shorter duration. Meert and Willems (2013) adapted the method to circumvent this problem:
N − 0.2 ⋅ S R = ϕ ⋅ P * = ϕ ⋅ (P − I a ) = ⋅ (P − 0.2 ⋅ S ) N + 0.8 ⋅ S
( 5.7 )
With N the daily average rainfall volume and P the actual amount of rainfall. Application of the SCS-method can be found in the German SMUSI-model (Kamradt, 2008) and the U.S. SWMM-model. Another, a bit more detailed modelling approach is to calculate the net rainfall by subtracting initial and permanent losses from the total rainfall input. Wetting of surfaces occurs at the beginning of a rain event resulting in a temporary storage of precipitation. After the wetting, depression losses arise due to the storage of rainfall in catchment depressions. The stored water is (at a later point in time) reduced by either evaporation or infiltration. Both wetting and depression losses can be considered as initial losses (Pi), since overland flow will only occur when the cumulated rainfall exceeds a certain threshold value that is defined by the wetting and depression storage. Another component that reduces the total rainfall are permanent losses (Pp) like infiltration and evapotranspiration. During rainfall events in urban catchments the latter is considered to be negligible when compared to infiltration, but it remains important for reducing the rainfall stored in depressions. The easiest way to incorporate the losses due to infiltration is to consider them as a constant percentage of the rainfall. Figure 5-1 illustrates the constant initial and permanent losses from a total bloc of rain, as it is used in the CITY DRAIN-package and the SEWSYS-model. (Achleitner, 2006; Achleitner, 2007a; Ahlman and Svensson, 2005; Kamradt, 2008)
Figure 5-1. Loss model with initial (Pi) and permanent (Pp) losses and net rainfall (Pnet). After (Achleitner, 2006). Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
47
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
In more detailed models (e.g. Vaes, 1999; Solvi, 2007), the evolution in time of the losses due to depression storage and infiltration are calculated. The initial losses are now restricted to the wetting losses. The filling of surface depressions will start once the rainfall exceeds the maximum storage volume available for wetting. An exponential filling is assumed for the depression storage model:
S (t ) = S max ⋅ [1 − exp(− c ⋅ I (t ))]
( 5.8 )
With the following runoff discharge:
Q(t ) =
S S max
⋅ I (t ) = [1 − exp(− c ⋅ I (t ))] ⋅ I (t )
( 5.9 )
Where Smax is the maximum depression height, c is the rate of storage loss and I is the rainfall height, fallen down after wetting. Equation ( 5.8 ) expresses that, the more rain has fallen, the less water can be stored in depressions. For pervious surfaces some extra losses need to be incorporated. The rainfall not only wets the pervious surfaces, but is also intercepted by vegetation. This combination can modeled as initial losses, but with different parameter values than for the impervious surfaces. Furthermore, the interception by vegetation is not constant but shows a seasonal variation. (Solvi, 2007) Another important factor for pervious surfaces is the infiltration of rainfall into the soil. Infiltration is a complex process where many parameters are involved (e.g. the grain size and saturation of the soil), but Horton’s equation can be used as a basic (and widely used) model for the infiltration capacity (Beven, 2004):
f (t ) = f min + ( f max − f min ) ⋅ exp(− k i ⋅ t )
( 5.10 )
With fmin the minimum infiltration capacity, fmax the maximum or initial infiltration capacity and ki a regression constant. An overview on a number of other equations to calculate the infiltration capacity can be found in (Borah, 2011). The infiltration process is initiated only if the rainfall intensity is greater than the steady state infiltration rate. In case of very low intensity rains infiltration equals the rainfall intensity. Depressions filling starts after the infiltration process, more precisely once the infiltration capacity is exploited. Figure 5-2 illustrates the different losses from a total bloc of rain falling onto a pervious surface. Optionally a runoff coefficient can be inserted to account for extra permanent losses like rainwater reuse (Vaes, 1999; Solvi, 2007; Saagi, 2013). The determination of net rainfall as described above is used in the KOSIM-model (Solvi, 2007), the REMULI-package (Vaes, 1999) and the QUURM-model (Watt and Kidd, 1975).
Figure 5-2. Runoff from pervious surface resulting from a block rain event (Solvi, 2007). Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
48
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Overland runoff routing Once the rainfall excess is determined, a hydrograph is developed by routing the excess rainfall to a point in the watershed. This component is primarily responsible for accounting the dynamics of an urban catchment in a model. (Borah, 2011) Different flow routing procedures are available to calculate the overland runoff flow and most of them are very familiar to the techniques discussed in the previous section on channel flow routing for river courses. The most widely used technique is that of the unit hydrograph. A unit hydrograph is defined as the hydrograph resulting from direct runoff produced by a unit of excess rainfall over a catchment. Hydrographs for storms of different intensity and duration are constructed using the principle of superposition (Zoppou, 2001). The most famous unit hydrograph is that of the linear reservoir model, where the catchment acts as a reservoir and the outflow is a linear function of the storage in that reservoir (Chow et al., 1988):
dS = I −Q dt
Q = S /k
( 5.11 )
or
S = k ⋅Q
The instantaneous unit hydrograph of a single linear reservoir is an exponential function:
Q(t ) =
1 − 1 ⋅ exp k k
( 5.12 )
A large catchment can be subdivided into equal sub-catchments with each considered as a separate linear storage. The instantaneous unit hydrograph for a cascade of n linear reservoirs is given by (Nash, 1960):
1 t ⋅ Q(t ) = k ⋅ (n − 1)! k
n −1
t ⋅ exp − k
( 5.13 )
Examples of the linear reservoir model or the Nash-cascade for simulating overland flow in urban catchments can be found in the QUURM model (Watt and Kidd, 1975), the FLUPOL model (Bujon, 1988), the SEWSIM model (Ruan and Wiggers, 1998), the HORUS model (Zug and Phan, 1999) the Cosmoss model (Calabró, 2001) and the KOSIM-software (Kamradt, 2007). The InfoWorks CS software (Wallingford, 2012a) that is used by most sewer managers in Flanders has implemented a double linear reservoir model (two identical linear reservoirs in series) to transform the net rainfall of each (sub-)catchment into an inflow hydrograph at each node. The reservoir constant k is expressed as a function of the catchment characteristics and the rainfall intensity:
k = C ⋅ i *−0.39
( 5.14 )
With:
C = 0.117 ⋅ S −0.13 ⋅ A 0.24 i* = 0.5 ⋅ (1 + i10 )
( 5.15 )
Where S (m/m) is the average slope of the catchment, A (m2) is the area of the surface contributing to the runoff and i10 (mm/min) is the running ten minute average of rainfall intensity. This last term is calculated in InfoWorks CS as follows:
i (t − 1) ⋅ RNIFAC − 1 + i (t ) i10 (t ) = max 0.03; 10 RNIFAC 0.3 ⋅ ∆t RNIFAC = max 1;4 − 60
( 5.16 )
With Δt the timestep, expressed in seconds. Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
49
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Another widely used unit hydrograph is that of the rational method, which has a trapezoidal shape (Berlamont, 2004). It is assumed that a fraction of the catchment contributes to the runoff, proportional in time from the start of the rainfall until the concentration time is reached. At this time the runoff discharge is equal to the peak discharge:
Qmax = ϕ ⋅ A ⋅ i M
( 5.17 )
With φ the runoff-coefficient, A the area of the catchment and i the rainfall intensity. If the rainfall lasts longer than the concentration time of the catchment, the peak discharge is held constant until the end of the rainfall. After this, the runoff discharge diminishes, again proportional with time. If the duration of the rainfall is shorter than the concentration time a trapezoidal unit hydrograph with a lower peak discharge is used (see Figure 5-3). Many other unit hydrograph models with similar principles are available, for example in the SCS-method (Mockus and Hjemlfelt, 2004), the Snyder unit hydrograph, Clark’s unit hydrograph (Borah, 2011) and the Colorado urban hydrograph procedure (UDFCD, 2008). Reference is made to literature for further information on these unit hydrographs. An extension of the rational method (or in fact any other unit hydrograph method) is the so-called time-area method. The time-area diagram indicates which part of the area of the catchment contributes to the total runoff at a certain time (see Figure 5-4). The peak discharge is then equal to the sum of flow contributions from subdivisions of the catchment defined by time contours (also called isochrones). This can be expressed as (Berlamont, 2004): t
Q(t ) = ∑ i r ⋅ At +1− r
( 5.18 )
r =1
With ir the net rainfall intensity and Ar the area of a subdivision r of the catchment. To calculate the overall peak flow it is assumed that the whole catchment is contributing to the flow when t is equal to the time of concentration of the catchment (Mehmood, 1995).
Figure 5-3. Unit hydrograph of the rational method (tp: duration of rainfall, tc: concentration time)
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
50
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 5-4. Time-area diagram and areas contributing to the total flow at the outlet point. (Berlamont, 2004)
In storage reservoir models the spatial derivative in the mass balance equation ( 5.11 ) is replaced by finite differences, so that:
S (t ) − S (t − ∆t ) = I −Q ∆t
( 5.19 )
S (t ) − S (t − ∆t ) I (t ) + I (t − ∆t ) Q(t ) + Q(t − ∆t ) = − ∆t 2 2
( 5.20 )
Or:
The outflow is again related to the storage in the reservoir, but with a nonlinear relation like for example the Manning equation (Borah, 2011). Different solution techniques for this system of non-linear equations were discussed in the previous chapter, whereof the storage-indication method (or level-pool routing method) seems the most appropriate in this case, since the relation between storage and outflow is most likely to be invariable over time, or during a rainfall event. The storage indication method rearranges equation ( 5.20 ) so that all unknown variables are on the left hand side:
2 ⋅ S (t ) 2 ⋅ S (t − ∆t ) + Q(t ) = + I (t ) + I (t − ∆t ) − Q(t − ∆t ) ∆t ∆t
( 5.21 )
The right hand side of equation ( 5.21 ) can be calculated with the results of the previous time step. A “storage-indication” curve, that relates the outflow Q(t) with the storage indication (2.S(t)/Δt + Q(t)), is then used to transform the left-hand side into an outflow value, from which follows the corresponding storage in the reservoir. A third reservoir model that is used for routing overland flow is the Muskingum reservoir. In the Muskingum reservoir, it is recognized that the storage in a river or reservoir depends on the inflow, as well as the outflow. It is assumed that the storage is a linear function of inflow and outflow, such that (Zoppou, 2001):
S (t ) = K ⋅ Q(t ) + K ⋅ X ⋅ [I (t ) − Q(t )] = K ⋅ [X ⋅ I (t ) + (1 − X ) ⋅ Q(t )]
( 5.22 )
In which K is the travel time of the flood wave through the reach and X is a dimensionless weighting factor between 0 and 0.5. The values of these parameters can be determined by trial and error or by calibration with observed hydrographs. Equation ( 5.22 ) can be substituted in the mass balance equation ( 5.20 ), which yields the following solution:
Q(t ) = C1 ⋅ I (t − ∆t ) + C 2 ⋅ I (t ) + C 3 ⋅ Q(t − ∆t )
( 5.23 )
Where C1, C2 and C3 are calculated as stated in the equations of ( 4.43 ). 5.1.2.
Sewer network
The sewer network collects the water from the urban catchment and transports it to a downstream point, which is either a waste water treatment plant and/or a receiving water body. Wastewater and stormwater Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
51
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
can be transported together in a combined or single pipe network or in a separate sewer network where wastewater and stormwater are not mixed. Three types of flows can be distinguished in the sewer network: gravitational flow with free surface, pressurized flow when the pipes are full and flows with backwater effects due to a limited throughflow. Combined sewer overflows can be considered as a special case of this last flow type. Conceptual models for these three flow regimes are discussed separately in this section. 5.1.2.1
Flow routing
Conceptual modelling techniques for routing the flow through sewer networks are very similar to the techniques discussed in the previous chapter on channel flow routing. These techniques can be broadly divided into three main classes: techniques based on (the unit hydrograph of) a linear reservoir model, multilinear techniques and non-linear reservoir models. Linear reservoir models The linear reservoir model is based on a mass balance equation and assumes a linear relation between the storage in the reservoir and the outflow from the reservoir (see equations ( 5.11 )). It was shown before that this differential equation can be solved as:
−1 −1 Q(t ) = exp ⋅ Q(t − ∆t ) + 1 − exp ⋅ I (t ) k k
( 5.24 )
The sewer system can now simply be represented by a number of linear reservoirs. The most elementary application of this technique is to replace the complete pipe and channel network by one single reservoir and use the precipitation falling on the impervious surfaces as inflow for the reservoir (e.g. Calabro, 2001). An even more elementary approach can be found in (Coutu et al., 2012), where the complete upstream catchment and the sewer network are represented by one linear reservoir. More extended forms of the technique (e.g. Sartor, 1999; Mannina et al., 2004) use a cascade of identical linear reservoirs, also called Nash-cascades (Nash, 1960). The parameters of these cascade models (the number of linear reservoirs and the value of the reservoir constant) were calibrated based on the results of a hydrodynamic model of the sewer system. The downstream discharge of a Nash-cascade is given by the gamma distribution:
∆t 1 Q(n ⋅ ∆t ) = ⋅ I (n ⋅ ∆t ) ⋅ k ⋅ (n − 1)! k
n −1
⋅ e − n⋅∆t / k
( 5.25 )
Gelormino and Ricker (1994) divide the sewer network into a number of interconnected virtual tanks. The volume in these tanks is explicitly calculated with the mass balance equation of the stored volume, the inflows and the outflow of the tank and the rain intensity input. The tanks are connected with flow paths in which the flow is based on the linear reservoir model approach: they are assumed to be proportional to the volume in the upstream tank. The main difference with the Nash-cascade is that the tanks do not need to be put in series, but can be arranged in any arbitrary way. (Ocampo Martínez, 2007) A similar technique is used in the P8 urban catchment model of the U.S. Environmental Protection Agency (Palstrom and Walker, 1990). The German KOSIM-model uses the Kalinin-Miljukov method (Euler, 1983), that is equivalent to the Nashcascade, but where the number of tanks n and the reservoir constant k are derived from the geometrical characteristics of the pipe. A pipe of length L is divided into a (natural) number of tanks with a characteristic length Lc (Solvi, 2007):
n=
L Lc
( 5.26 )
D Lc = 0.4 ⋅ S
Where D and S are respectively the diameter and the slope of the pipe. The reservoir constant is calculated as follows:
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
52
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
D2 k = 0.64 ⋅ L ⋅ Qmax
( 5.27 )
With Qmax the maximum discharge of a pipe, calculated with the Colebrook-White equation:
k s 2.51 ⋅ υ + ⋅ 2⋅ g ⋅ D⋅S Qmax = A ⋅ − 2 ⋅ log D ⋅ 2 ⋅ g ⋅ D ⋅ S 3.71 ⋅ D
( 5.28 )
In equation ( 5.28 ), A is the cross-sectional area, g the gravitational acceleration, ν the kinematic viscosity and ks the pipe roughness. The starting point of the traditional (or static) linear reservoir is the assumption that the outflow of the reservoir is proportional to the storage in the reservoir. In the case of the dynamic linear reservoir, it is assumed that the rate of change of storage is proportional to that of the outflow and the proportion factor is a function of time (Cluckie et al., 1999):
dQ dS = A1 (t ) ⋅ dt dt
A1 (t ) = k ⋅ t
( 5.29 )
1−α
With α a dimensionless parameter, that can be considered as a stability factor, and A1(t) a function that portraits the dynamic characteristics of the storage decay. Substitution of equations ( 5.29 ) into the mass balance equation ( 5.11 ) yields (Yuan, 1994):
dQ 1 α −1 = ⋅ t ⋅ [I (t ) − Q(t )] dt k
( 5.30 )
The unit pulse response function of equation ( 5.30 ) is given by:
Q(t ) =
− tα t α −1 ⋅ exp k α ⋅k
( 5.31 )
Or for n dynamic linear reservoirs in series:
Q(t ) =
α ⋅ t n⋅α −1
− tα exp ⋅ (α ⋅ k )n ⋅ (n − 1)! α ⋅k
( 5.32 )
The response of the reservoir to a given input hydrograph can be found by superposition of the unit hydrographs due to the inflows occurring in each unit period. The dynamic linear reservoir is applied in the SEWSIM model (Ruan and Wiggers, 1998) and the RHINOS model (Cluckie et al., 1999). Multilinear models Multi linear models use a combination of different linear relationships to breakdown the sewer system into multiple linear subsystems (Saagi, 2013). Most famous example of multi linear models is the Muskingum routing method (equation ( 5.23 )). The Muskingum routing procedure is applied in CITY DRAIN, an open source package for the simulation of integrated urban drainage systems (Achleitner et al., 2007a). Both the river model block, the catchment model block and the sewer model block of the package use the Muskingum method for flow routing. In case of the sewer block the user has the possibility to choose between one Muskingum reservoir or a series of n reservoirs (Achleitner et al., 2007b). The Muskingum routing method for sewer systems is also used in other software packages, for example in the French FLUPOL model (Bujon, 1988) and HORUS model (Zug and Phan, 1999) and the Australian MUSIC-model (MUSIC Development Team, 2012).
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
53
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Attempts have been made to calculate the values of K and X, based on the geometry of the sewer pipes. This is similar to the Muskingum-Cunge approach for rivers (Cunge, 1969). Mehmood (1995) proposes the following equations to calculate the Muskingum parameters:
K=
L
ω
( 5.33 )
1 Q X = ⋅ 1 − 2 B ⋅ S ⋅ L ⋅ ω
Where L and S are respectively the length and gradient of the pipe, Q is the discharge at normal depth (h), B is the surface width and ω is the kinematic wave speed along the pipe. This ω can be calculated as:
ω=
1 dQ ⋅ B dh
( 5.34 )
With Q defined by the normal depth relationship of the Colebrook-White equation:
14.8 ⋅ R Q = A ⋅ 32 ⋅ g ⋅ R ⋅ S ⋅ log ks
( 5.35 )
Where g is the gravitational acceleration, R the hydraulic radius, A the cross section of the pipe and ks the pipe roughness. Vaes (1999) developed a multilinear model that divides the total storage in the system into two (artificial) components: the static storage and the dynamic storage (see Figure 5-5). The static storage can be related to the throughflow by a linear or piecewise linear relationship (see Figure 5-6). The maximum static storage is defined as the storage when the system is filled up to the crest of the overflow and there is no throughflow. Although this is an unrealistic situation, this state is most often a very good approximation of the situation at the end of an overflow event. The instantaneous static storage can be assessed in the same way as the storage during the emptying of the system when no rainfall input is present. The remaining part of the storage is the dynamic storage and causes the hysteresis effect that is visible in Figure 5-6. The maximum dynamic storage is related to the maximum overflow discharge. Therefore a relation between the instantaneous dynamic storage and the global outflow from the system (overflow and throughflow) might exist. It is likely that there is also a relation with the inflow. Hence, the dynamic storage is probably described best by a relationship that combines the inflow and outflow. For practical reasons however, the dynamic storage is mostly only linked to the inflow (cfr. Muskingum method). For reasons of calculation simplicity, it was chosen by the developer to uncouple the dynamic storage from the static storage and put them in series (like a cascade). The dynamic storage component is considered as a linear reservoir model where the storage is proportional to the inflow, via a linear or piecewise linear relation. The outflow from this reservoir is then sent to the static storage reservoir, which is also modeled as a linear reservoir model. In this case the storage is proportional to the total outflow (overflow + throughflow).
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
54
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 5-5. Longitudinal profile of a sewer system on which the static and dynamic storage is indicated. (Vaes, 1999)
Figure 5-6. Piecewise linear ‘dynamic’ approach (with slopes kstat and kdyn) for the storage/throughflow - relationship of a small gravitary sewer system (for a composite storm which will just lead to an overflow event) (Vaes, 1999).
Nonlinear models Although the sewer flow can be described in a simplified way by a linear or multilinear reservoir model, the response in reality is hardly linear. To include various non linearities in the system, the output is modeled as a nonlinear function of storage. (Saagi, 2013) A storage reservoir with a very simple nonlinear storage throughflow relationship is discussed in (Gernaey et al, 2011):
Q = c ⋅ S 1.5
( 5.36 )
Where c is a constant parameter with units of L1.5 / T. The storage S in the system is calculated afterwards with the mass balance equation. Sewer hydraulics can then be simulated by a series of such nonlinear reservoirs. A similar non-linear model is used in the SEWSYS model (Ahlman and Svensson, 2005), but with a different exponent, namely 5/3.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
55
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 5-7. Partially linearized throughflow - storage relationship. After (Kamradt, 2008).
A different approach to a nonlinear reservoir method is described in (Kamradt, 2008). The main idea is that the nonlinear storage throughflow relation is replaced by a partially linearized relation. The nonlinear relation is in fact approximated by several straight lines with different gradients ki, that are only valid between two specific nodes (see Figure 5-7). The outflow of the reservoir is then calculated as:
Q(S (t )) = Qi −1 + k i ⋅ (S (t ) − S i −1 )
( 5.37 )
A third nonlinear model was proposed by Motiee (1997). The base of the model is again the mass balance equation, which is now coupled with a storage equation that links together the storage in the reservoir and the flow rates upstream and downstream. Assuming that the evolution of the wetted cross section is always increasing or always decreasing along the reach, the following storage equation can be written:
S (t ) = [α ⋅ AUS (t ) + (1 − α ) ⋅ ADS (t )] ⋅ Dx
( 5.38 )
Where AUS and ADS are respectively the wetted area upstream and downstream and α is a coefficient whose value depends on the flow characteristics. Notice that equation ( 5.38 ) is similar to the storage equation of the Muskingum method, equation ( 5.22 ). However, in the method of Motiee (1997) a different solution is obtained. In the case of free surface flow in the sewer network the value of α can be set to one, which simplifies equation ( 5.38 ) into:
S (t ) = AUS (t ) ⋅ ∆x
( 5.39 )
Or:
S (t ) =
I (t ) ∆x ⋅ ∆x = ⋅ I (t ) uUS (t ) uUS (t )
( 5.40 )
With uUS(t) the upstream flow velocity, which has to be calculated at every time step, for example by the Manning-Strickler formula. The main hypothesis is that the variations of the hydrograph are slow enough to obtain a quasi steady flow at each time step. Equation ( 5.40 ) can be rewritten with a temporal parameter Tp(t) = Δx / uUS(t) to obtain:
S (t ) = I (t ) ⋅ T p (t )
( 5.41 )
The downstream discharge is subsequently calculated with the mass balance equation:
Q(t ) = I (t ) −
S (t ) − S (t − ∆t ) ∆t
( 5.42 )
Since the relation between flow velocity and flow rate is not linear, this model can be considered as a nonlinear model. Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
56
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
5.1.2.2
Backwater Effects
The simple hydrologic models discussed in the previous section have one major disadvantage. Since they do not (explicitly) consider the influence of the downstream water level on the outflow of the pipe, special situations like pressurized flow, overtopping of street levels or backwater effects cannot be directly taken into account. The term backwater describes the situation that the maximum pipe flow is not sufficient to conduct the incoming flow downstream, or situations where a downstream boundary condition inhibits the maximum pipe flow. The excess water is then stored in the adverse direction of the flow by filling the retention volume in the pipe upstream. Too small pipes or an obstructing structure downstream are examples of situations that can cause backwater effects. The results of backwater can be twofold: inducing overflow events or activating storage volume and hence reducing overflow peaks. (Vanrolleghem et al., 2009) It is clear that these backwater effects have to be considered to achieve an accurate sewer network model. Different researchers have tried to integrate extra components in their conceptual models (all of them are already mentioned above), that can simulate backwater effects. A very common practice (e.g. Sartor, 1999; Solvi, 2007; Achleitner, 2007b) is to insert an additional static reservoir model to detent the flow volume higher than the maximum pipe capacity until the sewer has free capacity again. Figure 5-8 shows an example of the effect of such modifications to the outflow hydrograph: the peak of the hydrograph is limited to a maximum value and is prolonged in the time. Sartor (1999) first calibrates a number of static linear reservoirs with dynamic routing results without backwater conditions, resulting in a value for the maximum flow through the network. The next step is to calibrate the storage coefficient for the extra linear reservoir (Kh) which only affects the part of the hydrograph higher than maximum flow without backwater effects. The approach in the CITY DRAIN package (Achleitner, 2007b) is slightly different: the Muskingum routing procedure is put in series with a (virtual) storage reservoir. The user can set a maximum outflow and a maximum capacity for this (virtual) storage reservoir. All flows exceeding this maximum outflow value are stored in the storage reservoir until this reaches its maximum capacity and overflows at the manholes occur. Solvi (2006) developed a combiner-splitter approach, which is located on top of the reservoir cascade representing (a part of) the network (see Figure 5-9). The splitter, with a given maximum throughflow, sends any excess water back to the upstream combiner, that sums the water coming from the upstream pipe and from the downstream backflow. This sum is then used as input for the Kalinin-Miljukov model to calculate a new downstream discharge. The maximum throughflow can be based on hydrodynamic simulations or calculated with the Colebrook-White equation ( 5.28 ). A choice of two splitter models is given to the user: the first one separates the flow into a set flow to be entered by the user with the remaining flow leaving the splitter through the other side (absolute splitter), whereas the second splitter divides the flow into two fractions according to a given flow fraction parameter. Vanrolleghem et al. (2009) suggest to combine this splitter-combiner approach with a non-linear model for flow routing to obtain an accurate and numerically efficient conceptual model for sewer systems.
Figure 5-8. Comparison of hydrodynamic and conceptual routing methods for an event with extreme backwater conditions. (Sartor, 1999)
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
57
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 5-9. Backflow model withcombiner-splitter approach (Solvi, 2007)
The non-linear model of Motiee (1997) discussed in the previous section can also be extended to account for backwater effects in the sewer network. Consider two consecutive reaches from whom the second one affects the first. The problem to be solved is to calculate the correct values of V1(t) and Q1(t) = I2(t), so that equations ( 5.43 ) and ( 5.44 ) are both verified:
Q1 (t ) = I 1 (t ) −
S1 (t ) − S1 (t − ∆t ) = I 2 (t ) ∆t
S1 (t ) = f [I 2 (t )]
( 5.43 ) ( 5.44 )
Equation ( 5.44 ) represents the fact that the storage volume in the upper reach is dependent on the depth h at the entrance of the second reach. This depth is a function of the flow entering the second reach, according to equation ( 5.43 ). Equation ( 5.44 ) shows a continuously increasing relation between S1 and Q1 or I2: the greater Q1, the greater is h1 and the greater is S1. Equation ( 5.43 ) on the other hand shows a linearly decreasing relation between Q1 and S1. The system of equations ( 5.43 ) and ( 5.44 ) has one equilibrium value and it always exists. Different numerical solution techniques can be used to find this equilibrium. In the approach of Vaes (1999) an extra storage, above the crest level of the overflow structure is assumed. During overflow events the water level will rise above the crest level, which can lead to significant extra backwater effects and an increase of the total storage in the system. The storage above the crest level can be considered as dynamic storage. However, it is difficult to link this extra dynamic storage unambiguously to the inflow as proposed before, since the backwater effects are caused by the downstream boundary and not by the upstream boundary. Therefore an extra (again artificial) component of the storage is defined as the storage between the crest level and the maximum water level above the crest of the overflow (see Figure 5-10). This extra storage can be assumed to be a function of the overflow discharge. The extra storage can also be included in the static storage by using piecewise linear relationships for the storage/throughflow relationships (Figure 5-7). The throughflow is still computed as the outflow of a linear reservoir, as mentioned before.
Figure 5-10. Longitudinal profile of a sewer system on which the static, dynamic and extra storage during the overflow event are indicated (Vaes, 1999).
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
58
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
5.1.2.3
Combined Sewer Overflows
Combined sewer overflows (CSO’s) can be found in combined sewer systems at locations with a fixed maximum throughflow, for example upstream of a wastewater treatment plant with a maximum capacity. When the incoming (stormwater) discharge exceeds this fixed threshold, the water level in the CSO-chamber will rise until the maximum capacity is reached. When the chamber is full, the overflow device will divert the exceeding discharge into the receiving water body (Calabrò and Viviani, 2006). Basically two types of conceptual models for the CSO can be distinguished: models that use a maximum throughflow value and models that explicitly calculate the storage in the CSO chamber. The first type of models (e.g. Solvi, 2007; Coutu et al., 2012) consider a CSO-chamber as a special case of backwater effect, and therefore use a flow delimiter function to model the throughflow. This is a valid assumption because backwater effects will arise once the CSO-chamber starts filling. The most simple approach (Coutu et al., 2012) uses a linear low partitioning function to represent the CSO. The outflow of the CSO is equal to the inflow as long as the inflow doesn’t exceed the threshold value Qcrit. Once the inflow exceeds this threshold value the throughflow is limited to the threshold value and all the excess water is diverted out of the network (see Figure 5-11). A more realistic approach (Solvi, 2007) accounts for an increase of the throughflow with increasing inflow using a linear correction factor δ ≥ 1 (see Figure 5-11):
Qout =
δ −1 4
⋅I +
5−δ ⋅ Qcrit 4
( 5.45 )
The leftover of water (Qin – Qout) is again transferred over the overflow structure into the receiving water. The second type of models are based on the mass balance equation ( 5.46 ) to calculate the storage in the CSO-Chamber (Achleitner, 2007a):
S (t ) − S (t − ∆t ) = I − Qout − Qover ∆t
( 5.46 )
With S the storage, I the flow entering the CSO-chamber, Qout the throughflow and Qover the excess flow that goes over the overflow construction. All flows in ( 5.46 ) represent mean flows occurring during the time step Δt. The throughflow and overflow can be calculated with standard hydraulic formulas, based on the water level in the CSO-chamber (Berlamont, 2002):
Qout (t ) = C d ⋅ A0 ⋅ 2 ⋅ g ⋅ h
( 5.47 )
Qover (t ) = C d ⋅ b ⋅ g ⋅ (h − z c )
3/ 2
( 5.48 )
Where Cd is a discharge coefficient, A0 the cross section of the throughflow pipe, b and zc the width and crest level of the overflow device, g the gravitational acceleration and h the water level in the CSO-chamber.
(a) (b) Figure 5-11. Simulating the flow in a Combined Sewer Overflow: (a) flow delimiter function, modelled versus real flow, (b) mass balance model (Achleitner, 2007b). Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
59
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Examples of this modeling approach can be found in the COSSOM-model (Mehmood, 1995), the COSMOSS-model (Calabrò, 2006) the CITY DRAIN-package (Achleitner, 2007a and 2007b) and the KOSIM-software (Kamradt, 2008). The approach can also be used for other storage tanks in the sewer network like retention tanks, bypass tanks and pass-through tanks (Solvi, 2007). It also allows to model the effects of pumping stations, since the volume pumped out of the CSO chamber can easily be inserted in the mass balance equation ( 5.46 ). Instead of considering and simulating the CSO structure as a part of the sewer network, some models carry out a volume balance over the total network to estimate the overflow volume over the whole duration of the event. An example of this methodology can be found in SOBAL (Sewer Overflow BALancing method) or the method of Kuipers (Vaes and Berlamont, 1995). If the total flow in the sewer exceeds the capacity of a downstream throttle device, static offline storage in the network is activated. When this storage reaches its maximum, the remaining water is spilled through the overflow. The method thus gives a sequence of overflow volumes for a given sequence of block rain events. Dynamic effects during an event (overflow discharges, durations, etc.) cannot be considered (Kroll et al., 2007). 5.1.3.
Waste water treatment plants
Conceptual models for the hydrodynamics in waste water treatment plants (WWTP’s) do not seem to exist. Simplified and conceptual models for WWTP’s are mainly focused on the water quality (e.g. …). Considering the structure of WWTP’s, which consist of a number of tanks (e.g. aeration and sedimentation tanks), they might also be modelled with storage reservoirs with a more or less constant throughflow to the river system. In the case that the storage in the sewer system is significantly larger than the storage in the WWTP, it might be possible to neglect the influence of the WWTP.
5.2. Software packages for integrated urban catchment modelling Over the last two decades a number of models and software packages for simulating the water quantity and quality in urban catchments have been developed. Most of these packages are mainly focused on the water quality part, but all of them need a hydrologic part for the transport of the water and pollutants. This section gives a brief description of the hydrological component of a couple of (academic) integrated urban catchment models. A detailed description is not provided, since the different components of the models have already been discussed in the previous section. Different commercial packages are available as well, but not further discussed, because of the lack of detailed information on them. Examples of commercial packages are SIMBA (Erbe et al., 2002), SMUSI (Muschalla et al., 2006) and SIMPOL3 (Freni et al., 2012). 5.2.1.
The reservoir modelling system Remuli
Based on the static-dynamic approach, discussed before in section 5.1.2, Vaes (1999) developed a physically based, conceptual modelling system for the urban catchment called Remuli. This section is therefore mainly based on his PhD-dissertation. The Remuli system is also implemented in the InfoWorks RS software to model the urban response to rainfall in a simplified way (Wallingford, 2012b). The Remuli system is designed for long term simulations of the urban wastewater system, with rainfall input with a 10 minutes time step, as it is available at Uccle, where the Royal Meteorological Institute of Belgium (RMI) is located. The model consists of four main parts that will all be discussed briefly: a runoff module, a smoothing module and two modules for the dynamic and static storage (see Figure 5-12).
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
60
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 5-12. Schematic representation of the reservoir modelling system Remuli (Vaes, 1999).
The runoff module incorporates the rainfall losses from impervious and pervious surfaces to compute the rainfall input into the system. For the pervious and impervious surfaces a fixed runoff-coefficient and depression storages can be specified, with different parameter values for each surface. An extra possibility is to model infiltration at the pervious areas with the Horton equation. The losses due to evaporation are kept constant, since they have not much of an impact on the prediction of overflows. The smoothing module contains an averaging of the runoff (with a time step of 10 minutes) over the concentration time of the urban catchment. In this way the lag time between inflow and outflow peak are well determined and part of the smoothing due to the flow in the system is already incorporated. This is not possible with a single reservoir, since this has an instantaneous response. Optionally, the concentration time can be a related to the peak discharge. It has been seen that the concentration time is smaller for larger discharges. The routing of the flow through the sewer network is done by using one or two linear reservoirs in series: the first one for the dynamic storage (which is in fact optional) and a second one for the static storage. For the dynamic storage reservoir a (piecewise) linear relationship between inflow and storage is assumed, while two (piecewise) linear relations between the outflow, the throughflow and the storage are used for the static storage reservoir. 5.2.2.
SEWSYS – tool
The SEWSYS – tool was developed at the University of Technology in Gothenburg (Sweden) as a part of the Mistra programme that concerns strategic environmental research (Ahlman & Svensson, 2005). The tool is mainly intended for the simulation of (up to 20) substance flows in urban sewer systems, but also contains a number of components to simulate the flows of water in this system. The input of the model consists of sanitary wastewater and storm water. A vector with 24 values controls the daily variation of sanitary wastewater discharges. The hourly values can be given an additional variation by multiplying them with a random number between 0.7 and 1.2, to achieve a variation within the hour when time steps smaller than one hour are used. Storm water due to rainfall-runoff only includes the runoff from impervious surfaces. The precipitation on these surfaces is first reduced with an initial loss and afterwards with a reduction factor that corresponds to the runoff-coefficient. The sum of sanitary waste water and the runoff is used as the input for the sewer model. Routing of the overland flow is thus not incorporated in the model. Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
61
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Both combined and separated sewer systems can be modelled in the SEWSYS-tool, although the routing technique is the same for both systems: a non-linear reservoir. In a first step the outflow is calculated with equation ( 5.36 ), where the exponent has a value of 5/3 and the parameter c is related to the slope, length and Manning’s roughness number of the considered pipe. However, estimation of its value based on hydrodynamic model results is also possible. The second step is used to update the storage in the reservoir, by using the mass balance equation ( 5.19 ). 5.2.3.
The KOSIM-West Modelbase
The German KOSIM modelling tool (ITWH, 2012) is designed for long-term simulations of dry weather flow generation, rainfall-surface runoff and transport (water and pollutant loads) in the sewer system. The model components were implemented in the WEST® software platform, developed at the university of Ghent, by Solvi (2007). Figure 5-13 gives a general overview of the processes and structures that are contained in the West® platform. Two main components can be distinguished in KOSIM: the catchment model and the sewer model. The catchment model computes the water entering the sewer network, which is composed of rainwater from surface runoff and wastewater from households, commercial and industrial sites (dry weather flow). The rainfall losses due to wetting, depression storage, infiltration and evaporation are modelled as depicted in Figure 5-2. A sine function is used to represent the variation of potential evaporation during the year. Once the net rainfall is determined surface runoff towards the sewer system is modelled with a cascade of 3 linear reservoirs in series. The reservoir constant of these linear reservoirs is estimated as a quarter of the concentration time of the catchment. Summation of the outflow from the cascade of linear reservoirs, the dry weather flow and optionally an infiltration directly into the sewer generates the input for the sewer network (Solvi, 2007; Kamradt, 2008). Routing of the pipe flows in the sewer network is done with the combiner-splitter approach (Figure 5-9). The Kalinin-Miljukov method with a cascade of linear reservoirs is applied to simulate the propagation of the flow. The parameters of this model are derived from geometrical characteristics of the main collectors (equations ( 5.26 ) to ( 5.28 )). If the flow exceeds a threshold value the excess water is sent back to the upstream reservoir, where it is routed again. CSO structures can be modelled in a similar way: a flow delimiter function is used that relates the throughflow with the inflow and transfers the leftover to the overflow weir. Another possible approach is to model them as a storage tank and explicitly calculate the stored volume of water (Solvi, 2007).
Figure 5-13. Elements and processes within the KOSIM-WEST model. (Solvi, 2007).
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
62
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
5.2.4.
Cosimat
The first version of Cosimat (Continuous Simulation Model for sewer systems based on Matlab) was built about 15 years ago. It is largely based on the previously discussed Kosim model and implemented in the Matlab/Simulink platform. Cuppens (2006) amended, refined and simplified the modelling principles and established a strategy for calibration of the different model components. The inflow to the sewer network consists of rainfall-runoff and dry weather flow. The net rainfall is determined by subtracting the evapotranspiration depth and continuous losses, due to surface wetting and depression storage, from the actual rainfall depth. Both evapotranspiration and the continuous losses are kept constant. Runoff to the sewer network is then calculated by applying a fixed runoff coefficient to the net rainfall. This procedure is thus familiar to that of Figure 5-1. A diurnal variation, based on the number of inhabitants and their average waste water production, is used for the dry weather flow. Optionally, industrial waste water contributions can be included. A reservoir cascade consisting of three reservoirs is used to simulate the routing process of both surface and in-sewer flow. The reservoirs can be accommodated with a variable reservoir constant, to allow for the influence of the rainfall intensity on the reservoir constant and hence the differences in routing behavior during events of a high and low intensity. After the water has been routed through the reservoir cascade, it leaves the subbasin via control structures like conduits, sluice gates, orifices, pumps or weirs. The flow through this control structure is calculated based on the water levels at both sides (Kroll, 2014). 5.2.5.
The CITY DRAIN© - package
CITY DRAIN is an open source software for integrated urban drainage modelling, developed at the university of Innsbruck (Achleitner, 2006 and 2007a). The software uses the MATLAB/SIMULINK environment, which enables a block wise modelling of the different parts of the system. Blocks for simulating the hydraulics of catchments, rivers and sewers and for the mass transport of pollutants have been implemented within CITY DRAIN. Rainfall runoff generation is done with a combination of a loss model and a flow model. The loss model accounts for an initial loss and permanent losses. A simple basin methodology is used where runoff occurs when the rain volume exceeds the volume of the basin. The loss model is followed by a flow model that converts the effective rainfall to flows using the area of the catchment and a runoff-coefficient. Muskingum reservoirs are used for routing the different flow types: overland flow in the catchment, river flows and flows in the sewer network. The Muskingum scheme was modified to model a number of reservoirs is series and to allow an extra input to each individual sub reach of this cascade. Backwater effects are covered by inserting an extra storage reservoir at the end of the cascade. This reservoir stores the excess water that is held because of the maximum throughflow value. CSO’s are modelled using the mass balance equation in discrete formulation (equation ( 5.43 )). Next to the just mentioned components the CITY DRAIN package contains some other components from the urban drainage system like source blocks for rain input, pumping stations, a black-box WWTP model and a hydropower station.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
63
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 5-14. Screenshot of the CITY DRAIN 2.0 block library (Achleitner, 2007b).
5.3. Discussion Several conceptual modelling techniques for simulating the urban drainage system were discussed in this chapter. Most of these techniques originate or are adapted from the techniques in the previous chapter on conceptual models for water courses. This section summarizes the similarities between these two systems and proposes a modelling technique for urban drainage systems that is based on the conclusions of the previous chapter. This should allow the simulation of both water courses and urban drainage systems with one generalized modelling technique. 5.3.1.
Link with river courses
The techniques proposed for routing overland flow in the urban catchment and routing sewer flows in the case where no backwater effects are present (unit hydrograph models, linear reservoirs, Muskingum reservoirs, Kalinin-Miljukov cascades, …), are also used for channel routing of river courses without significant backwater effects (section 4.1). It was shown there that a unified coefficient routing model exists that expresses the discussed models in one single formulation ( 4.54 ). This unified coefficient routing model on his turn can be further generalized in the linear transfer function model of equations ( 4.58 ) and ( 4.59 ), which covers a broad range of possible model structures. Furthermore, it was shown by Willems (2000) that the model structure of a linear transfer function can be regarded as an arrangement of linear reservoirs. The first application of transfer function modelling to urban drainage systems was performed by Capodaglio (1994). Transfer functions were used in two test cases for both the overland runoff routing and for predicting flows at a downstream point in the sewer network. Wolfs (2013) proved that the static-dynamic approach of Vaes (1999) can be expressed as a transfer function model, with two linear reservoirs in parallel and a time delay of one time step on the second reservoir. For reasons of calibration and computation this is simplified to two linear reservoirs in series, in the REMULI-system. Two more complex model structures (higher order transfer function models and the use of piecewise linear relations between storage and throughflow) were proposed to achieve more accurate results. Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
64
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
As stated before, the response of sewer networks can be strongly non-linear, for example because of backwater effects, pressurized flows or the overtopping of manholes. This might cause problems when applying the linear transfer function model to sewer networks, since this is considered to be a linear model when the parameter values are held constant. A possible solution is to make (one or more) parameter values state dependent, since it can be assumed that the phenomena causing the nonlinear behaviour (e.g. backwater effects) will arise at a certain discharge and become increasingly important at higher discharges. The non-linear model of Figure 5-7, the flow delimiter function of Figure 5-11 and the piecewise linear relations of Wolfs (2013) can all be considered as an example of a state dependent transfer function. The value of the proportionality factor between storage and through flow will not be constant and will depend on the value of the storage, or the discharge. However, it seems that this state dependent transfer function has not been thoroughly tested and applied on sewer systems so far. Besides channel routing techniques, conceptual models for sewer networks also make use of reservoir routing techniques. For example when dealing with retention tanks, combined sewer overflows or backwater effects. The static-dynamic approach can in fact also be solved with reservoir routing, since in the original solution scheme the storage is calculated at every time step. The same solution techniques as the ones mentioned in the previous chapter can then be used to compute the flow rates in the sewer system. 5.3.2.
Proposed modelling technique
Based on the analysis above, a modelling technique for sewer networks similar to that for river courses is proposed. The technique is mainly focused on simulating the flow propagation in the sewer network, when the inflow hydrograph from the urban catchment is known. The proposed methodology still has to be tested to a number of test cases and might need some modifications. The modelling approach starts with identifying the locations where important backwater or storage effects occur, like retention tanks, locations with overflow structures and sewer pipes with a limited throughflow. As many as possible sewer pipes located upstream of these locations will be modelled with one storage reservoir, so that throughflow and overflow can be related to the storage in the reservoir. When too many upstream pipes are selected, an important hysteresis will show up in the storage-discharge relation. This is caused by the static-dynamic principle explained in section 5.1.2.1. in this case, a transfer function can be inserted to model the transformation of the input hydrograph in the more upstream parts of the sewer network to the point where the backwaters occur. Initially, a simple linear transfer function model can be used to model this transformation. If this is not sufficient, e.g. because of pressurized flows, use can be made of a nonlinear transfer function with state dependent parameters. The objectives (accuracy and speed) of the conceptual model will determine to which extent the urban drainage system will be modeled. If someone is interested in the (order of) magnitude of the overflow volumes to the adjacent river, the model can be faster but less detailed than in the case where the modeler is interested in the flow rates at different locations in the sewer network.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
65
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
6. Estuaries A classical definition of estuaries is that proposed by Pritchard (1976): “An estuary is a semi-enclosed coastal body of water which has a free connection with the open sea and within which sea water is measurably diluted with fresh water derived from land drainage”. The connection of an estuary (or a tidal river) with the open sea implies that the conceptual models discussed in the previous chapters, where the water levels are mainly determined by the discharge, are no longer valid. On the contrary, in estuaries and tidal rivers the main factors influencing the water levels are the tidal amplitude, the storm surge at the inlet and possibly the fresh water inflows. This chapter discusses four simplified models that can be used to determine the tidal propagation in estuaries and tidal rivers. Only the first two (the Muskingum water-stage methods and the CALAM-model) can be classified as conceptual models. The last two (idealized models and black box models) are included since they have previously been applied to the Western Scheldt estuary (and some other estuaries) and they also allow a fast calculation of water levels along tidal rivers and estuaries.
6.1. Muskingum water-stage methods The Muskingum water-stage routing methods (e.g. Franchini and Lamberti, 1994; Si-min et al., 2009) are based on the same principle as the classical Muskingum routing method (section 4.1.6): the total storage in a river reach can be split up in two parts: prism storage and wedge storage (see Figure 4-6). The main difference between the two methods is that the classical Muskingum method describes these two storages in terms of discharges, while the water-stage method calculates the storage in a reach in terms of water levels. This adaptation should make it possible to introduce a water level boundary at the downstream side and hence the application to tidal rivers. The classical Muskingum routing method is based on the water balance and a linear relation between storage and in- and outflow (USACE, 1994):
dS = I −Q dt
( 6.1 )
S = K ⋅ Q + K ⋅ X ⋅ (I − Q ) = K ⋅ [X ⋅ I + (1 − X ) ⋅ Q ]
( 6.2 )
The solution of equations ( 6.1 ) and ( 6.2 ) is a linear function:
Q(t ) = C 0 ⋅ I (t − ∆t ) + C1 ⋅ I (t ) + C 2 ⋅ Q(t − ∆t )
( 6.3 )
With I the inflow, Q the outflow and C0, C1 and C2 constant parameters of the model that can be expressed as a function of the travel time K, a weighting factor X and the time step Δt. This classical formulation ( 6.3 ) can be transformed in a water-stage routing method by transforming the discharges to water levels, with so-called rating curves. This can be done for example by using the Manning formula (Berlamont, 2004):
Q=
1 1/ 2 ⋅ S0 ⋅ B ⋅ H 5 / 3 n
( 6.4 )
Where n is the roughness coefficient and S0 the water surface slope. The inflow at the upstream section and the outflow at the downstream section can then be expressed in a similar, but generalized way as (Franchini and Lamberti, 1994):
I = aU ⋅ h bU
( 6.5 )
Q = a D ⋅ H bD
Where the parameters a and b can be calibrated based on dynamic model results or measurements. For reaches without tidal influence these parameters are assumed to be constant, because of the time-invariant Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
66
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
stage-discharge relation. Combining equations ( 6.3 ) and ( 6.5 ) yields the following non-linear solution for the downstream water level H2:
]
[
1 ⋅ C 0 ⋅ aU ⋅ h2bU + C1 ⋅ aU ⋅ h1bU + C 2 ⋅ a D ⋅ H 1bD H2 = aD
1 / bD
( 6.6 )
Where C0, C1 and C2 are the same constant parameters as in ( 6.3 ) and subscripts 1 and 2 represent the start and end time of the routing interval. The main purpose of this model was to simulate and forecast water levels, when only recorded water level data are available. Since this model is in fact a simple adaptation of the classical Muskingum routing method, it will probably only be valid for river courses without backwater effects or tidal influences. Furthermore it is only valid for relatively slow descending waves, compared with the characteristic time of the channel (Franchini and Lamberti, 1994). Another approach is to express the storage in a reach in terms of the upstream and downstream water level (Si-min et al., 2009):
S = L ⋅ B ⋅ H = L ⋅ B ⋅ [x 0 ⋅ h + (1 − x 0 ) ⋅ H ]
( 6.7 )
With L the length of the reach, B the average width of the channel, h and H the upstream and downstream water depth and x0 a water depth weighted coefficient. To combine equations ( 6.7 ) and ( 6.1 ) the discharges and water depths must be linked with the rating curves of ( 6.5 ). Solving this system of equations yields the following solution for H2:
A0 a ⋅ (1 − x0 ) ⋅ H 2 + D ⋅ H 2bD = Dt 2 aU A ⋅x A a ⋅ h1bU + h2bD + 0 0 ⋅ (h1 − h2 ) + 0 ⋅ (1 − x0 ) ⋅ H 1 − D ⋅ H 1bU Dt Dt 2 2
)
(
( 6.8 )
Where A0 = L . B is the average area of the cross section over the length of the considered reach. Equation ( 6.8 ) is a nonlinear and implicit function for H2, because of the exponent bD. However, under the assumption of linear stage-discharge relations (bU = bD = 1 and aU = aD) the solution can be simplified into:
H 2 = C 0 ⋅ h2 + C1 ⋅ h1 + C 2 ⋅ H 1
( 6.9 )
with:
0.5 ⋅ ∆t − C0 =
A0 ⋅ x0 aU
A0 ⋅ (1 − x0 ) + 0.5 ⋅ ∆t aU A 0.5 ⋅ ∆t + 0 ⋅ x0 aU C1 = A0 ⋅ (1 − x0 ) + 0.5 ⋅ ∆t aU A0 ⋅ (1 − x0 ) − 0.5 ⋅ ∆t aU C2 = A0 ⋅ (1 − x0 ) + 0.5 ⋅ ∆t aU
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
( 6.10 )
WL2014R00_131_1
67
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 6-1. Accuracy of water-stage simulating with bidirectional Muskingum water-stage method, during validation phase (Si-min et al., 2009).
The stage routing method of equations ( 6.9 ) and ( 6.10 ) is also not applicable to tidal rivers, since the stage-discharge relations of these reaches vary in time and are strongly non-linear. Si-min et al. (2009) therefore suggested to make the parameters C0, C1 and C2 time-dependent by introducing two extra parameters KE and XE:
0.5 ⋅ ∆t − KE ⋅ XE KE ⋅ (1 − XE ) + 0.5 ⋅ ∆t 0.5 ⋅ ∆t + KE ⋅ XE C1 = KE ⋅ (1 − XE ) + 0.5 ⋅ ∆t KE ⋅ (1 − XE ) − 0.5 ⋅ ∆t C2 = KE ⋅ (1 − XE ) + 0.5 ⋅ ∆t C0 =
( 6.11 )
With:
K0 =
A0 aU
H max − h XE = x0 ⋅ exp − AX (h − H min ) H − h + AK ⋅ HM KE = K 0 ⋅ max h − H min + AK ⋅ HM
( 6.12 )
The model now contains seven parameters: the travel time K0, the water depth weighted coefficient x0, the maximum and minimum water stage Hmax and Hmin, the average water-stage HM, an influence index AK and a power index AX. Equations ( 6.9 ), ( 6.11 ) and ( 6.12 ) compose the linear, time-dependent Muskingum water-stage routing model that can be used for tidal rivers. However, the model is still a linear model and the results, as presented by Si-min et al. (2009), are rather poor (see Figure 6-1 for an example). Differences between the observed and predicted water levels up to 1.5 m can be noticed.
6.2. CALAM – model The CALAM – model was developed in the Netherlands as part of a master thesis at the technical university of Delft (Urbanus and Vreeburg, 1987). The main intention of the model was a fast simulation of the water movements in the Rhine-Meuse delta, combined with a water quality model that could simulate the propagation of an accidental pollution. The CALAM - model is a reservoir type model and has the advantage that both upstream discharges and downstream water levels can be used as boundary
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
68
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
conditions. As such the model should also be applicable to estuaries, tidal rivers and rivers with a downstream water level boundary.
Figure 6-2. Schematization of the Rhine-Meuse Delta in the CALAM - model, with nodes and segments (Urbanus and Vreeburg, 1987)
The studied estuary or tidal river is replaced by a system of nodes and segments (see Figure 6-2 for an example). Every segment of the system is represented by a reservoir with a constant trapezoidal cross section. The nodes make up the connections between the different segments and can be found at the boundaries of the system, at the intersection of multiple segments and at intermediate points. Water movement is modelled with discretized formulations of the one-dimensional shallow water equations: the continuity equation and the momentum equation. The former is solved in the nodes, while the latter is solved in the segments. The derivation of the discretized form of both equations is described succinctly below. A more detailed derivation can be found in (Urbanus and Vreeburg, 1987). The continuity equation in a vertical cross section can be formulated as:
∂u ∂w + =0 ∂x ∂z
( 6.13 )
With u the flow velocity along the channel (x-direction) and w the vertical velocity (z-direction). Integration of equation ( 6.13 ) over the water depth yields a more familiar form of the continuity equation:
∂Q ∂h + B⋅ =0 ∂x ∂t
( 6.14 )
Where Q is the discharge, B the width of the cross section, h the water level above datum and t the time. The water level h can be replaced by the water depth a when the bottom level of the cross-section is fixed:
∂Q ∂a = −B ⋅ ∂x ∂t
( 6.15 )
The momentum equation for a flow in the x-direction is written as:
∂u ∂u 2 ∂h g ⋅ u ⋅ u + +g⋅ + =0 ∂t ∂x ∂x C2 ⋅a
( 6.16 )
With g the gravitational acceleration and C the Chézy coefficient. The last term in equation ( 6.16 ) represents the friction slope. Equation ( 6.16 ) is also integrated over the water depth and over the width of the cross section to obtain:
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
69
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
∂h g ⋅ Q ⋅ Q ∂Q ∂ Q 2 =0 + + + g ⋅ A⋅ ∂x C 2 ⋅ A ⋅ R ∂x A ∂t
( 6.17 )
Where A and R are respectively the area and hydraulic radius of the cross section. Now, considering that:
∂h ∂z 0 ∂a ∂a = + = −S 0 + ∂x ∂x ∂x ∂x 2 ∂ Q 2 ⋅ Q ∂Q Q 2 ∂A ⋅ − ⋅ = ∂x A A ∂x A ∂x
( 6.18 )
and that the cross section of each segment is constant, equation ( 6.17 ) can be rearranged and simplified into:
Q ⋅Q 1 ∂Q 2 ⋅ Q ∂Q ∂a =0 + 2 − S0 + ⋅ + ⋅ 2 ∂x C ⋅ A 2 ⋅ R g ⋅ A ∂t g ⋅ A ∂x
( 6.19 )
Substitution of equation ( 6.15 ) yields the final formulation of the momentum equation, which will be used in the model:
Q ⋅Q ∂a 1 ∂Q 2 ⋅ Q ⋅ B ∂a ⋅ − ⋅ − S0 + + 2 =0 2 g ⋅ A ∂t ∂t ∂x C ⋅ A 2 ⋅ R g⋅A
( 6.20 )
In the next step, the continuity equation ( 6.14 ) and the momentum equation ( 6.20 ) are discretized to solve them as a system of linear equations and to obtain the discharges and water levels in the system. After integrating equation ( 6.14 ) over the length of the adjacent segments i, the continuity equation in node j can be expressed as:
∑ Qi = F (h j )⋅ i
dh j
( 6.21 )
dt
With Q representing the average incoming and outgoing flow rates, F(hj) the horizontal surface of node j according to the water level in this node, hj. The horizontal surface of a node is calculated as the sum of half the horizontal surfaces of each adjacent segment. Discretization of equation ( 6.21 ) yields:
h j (t ) − h j (t − ∆t ) ∆t
θ ⋅ ∑ Qi (t ) + (1 − θ ) ⋅ ∑ Qi (t − ∆t ) = F (h j ) ⋅ i
i
( 6.22 )
With θ a constant factor between zero and one (usually 0.5 as in equation( 4.67 )). The different terms in equation ( 6.22 ) can be linearized as follows:
h j (t ) = h j (t − ∆t ) + ∆h j
( 6.23 )
Qi (t ) = Qi (t − ∆t ) + ∆Qi Which yields the following equation:
∑ Q (t − ∆t ) + θ ⋅ ∑ ∆Q i
i
i
i
∆h j = F (h j ) ⋅ ∆t
( 6.24 )
Where ΔQi and Δhj are the unknowns, respectively the change of the discharge in segment i and the change of the water level in node j, over the time step Δt.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
70
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
The discretized formulation of the momentum equation is derived in a similar way. After integrating over the length of a segment i, discretization and linearization of h and Q and rearranging, equation ( 6.20 ) becomes:
L g ⋅ A ⋅ ∆t + 2 ⋅ θ ⋅ ε ⋅ Qi (t − ∆t ) ⋅ ∆Qi + Qi (t − ∆t ) ⋅ B ⋅ L 1 dε + ⋅θ ⋅ ⋅ Qi (t − ∆t ) ⋅ Qi (t − ∆t ) − θ ⋅ ∆hi ,1 + − 2 2 dh g ⋅ A ⋅ ∆t
Qi (t − ∆t ) ⋅ B ⋅ L 1 dε + ⋅θ ⋅ ⋅ Qi (t − ∆t ) ⋅ Qi (t − ∆t ) + θ ⋅ ∆hi , 2 = − 2 2 dh g ⋅ A ⋅ ∆t hi ,1 − hi , 2 − ε ⋅ Qi (t − ∆t ) ⋅ Qi (t − ∆t ) + S 0 ⋅ L
( 6.25 )
With hi,1 and hi,2 the water levels at both sides of segment i, Qi the discharge in the segment and ε = L/(C²A²R). The elaborated derivation of equation ( 6.25 ) can be found in (Urbanus and Vreeburg, 1987). Both equations ( 6.24 ) and ( 6.25 ) can be rewritten more distinctly as follows:
Aa ⋅ ∆Qi + Ab ⋅ ∆h j = 0
Ba ⋅ ∆Qi + Bb ⋅ ∆hi ,1 + Bc ⋅ ∆hi , 2 = 0
( 6.26 )
Which yields a system of equations that can be expressed as a matrix:
a11 a1 a 1 a a m1 a m
a1n ∆Q1 b1 a n ∆Q b ⋅ = a mn ∆Qn bn
( 6.27 )
All values of ΔQi can be obtained by solving this system of equations. The values of Δhj follow then from equation ( 6.26 ). As mentioned before the CALAM – model can handle two types of boundary conditions: water levels at the downstream point and discharges at the upstream side. Both can be constant values or time series. The initial conditions of the model consist of the water levels in each node and the assumption of discharges equal to zero in all segments. As a result of this last initial condition the CALAM – model needs a warming up period of approximately 24 hours. Solving the matrix equation of ( 6.27 ) involves the calculation of the inverse of matrix aij. This might become complicated and time consuming, when the size of the matrix and the number of non-zero elements increases. A solution to this problem might be to neglect the upstream discharges and thus only considering the propagation of the tide in the estuary and river. This should allow to achieve a straightforward calculation scheme that proceeds from the downstream side to the upstream side. A number of tests still need to be conducted to confirm this.
6.3. Idealized models Idealized models are based on the cross-sectionally averaged shallow water equations for water motion and assume highly schematized geometries and boundaries. Their relative simplicity makes it possible to analyze the models by standard mathematical techniques, resulting in analytical solutions, without neglecting the physical mechanisms that control the dynamics in estuaries (Schramkowski et al., 2004). The models are mostly used for simulating the propagation of the tide in estuaries and tidal rivers (e.g. Speer and Aubrey, 1985; Friedrichs and Aubrey, 1994; Lanzoni and Seminara, 1998) and for long-term morphodynamic simulations (e.g. Schuttelaars and de Swart, 2000; Hibma et al., 2003).
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
71
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Figure 6-3. Top view of an idealized estuary or embayment. Exponential varying width (left; Todeschini et al., 2008) and rectangular basin (right; Schuttelaars and de Swart, 2000).
Figure 6-4. Diagram of a cross section with tidal flats as used in idealized models.
Idealised models consider estuaries and tidal rivers as straight channels closed at one end and connected with a tidal sea at the other end. The channel is assumed to have a length L, with a cross section as depicted in Figure 6-4. The width of the channel b0 can vary along the longitudinal direction by an exponential decrease or can be held constant, resulting in a rectangular basin (Figure 6-3). Since the boundary on the landward side is assumed to be closed, the effects of fresh water inflows are neglected. This is an acceptable approximation in case of the Western Scheldt, because the river discharge at the Dutch-Belgian border is smaller than one percent of the tidal prism (Schuttelaars & de Swart, 2000). The channel width is considered to be much larger than the reference water depth, which implies that the water motion can be described by the cross sectionally integrated, one-dimensional shallow water equations. For a tidal channel with linearly sloping intertidal flats (like in Figure 6-4) these equations can be expressed as (Friedrichs and Aubrey, 1994):
∂h ∂ + [A ⋅ u ] = 0 ∂t ∂x ∂u ∂u ∂h +u⋅ +g⋅ +F =0 ∂t ∂x ∂x
b⋅
( 6.28 ) ( 6.29 )
Where b is the total estuary width, h is the water surface elevation (as depicted in Figure 6-4), A is the channel cross section area, u is the cross sectionally averaged velocity along the direction of the channel, g is the gravitational acceleration and F represents the bottom friction. Furthermore it is assumed that the flow velocity on the tidal flats can be neglected. The boundary conditions at the entrance x = 0 are restricted to the case of simple harmonic tidal waves that can be represented by trigonometric functions (e.g. Speer and Aubrey, 1985; Lanzoni and Seminara, 1998):
h(0, t ) = a ⋅ cos(ω ⋅ t )
( 6.30 )
Schuttelaars and de Swart (1997, 2000) used a basic tide (M2) and a first overtide (M4) as the boundary conditions for their idealized model of the Western Scheldt estuary: Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
72
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
h(0, t ) = cos(t ) + β ⋅ cos(2t − φ )
( 6.31 )
With β the ratio of the M4 and the M2 amplitudes at the entrance and φ (= φM4 – 2.φM2) the phase difference between the M4 and M2 tidal constituents. Note that this last boundary condition is expressed as a dimensionless function. Elementary idealized models (e.g. Schuttelaars and de Swart, 1997) assume a very simplified form of the momentum equation ( 6.29 ):
g⋅
∂h +F =0 ∂x
( 6.32 )
Combining equations ( 6.30 ) and ( 6.32 ) yields a first order wave equation for tidal elevation. This wave equation has a trigonometric solution with a constant amplitude:
h( x, t ) = a ⋅ cos(ω ⋅ t − k ⋅ x )
( 6.33 )
whereas the amplitude of a classical cooscillating tide in a finite channel undulates along the channel (Friedrichs and Aubrey, 1994). More detailed models (e.g. Lanzoni and Seminara, 1998; Schuttelaars and de Swart, 2000) use the complete momentum equation to account for bottom friction, reflection of the tidal wave, length of the estuary and the convergence and dissipation of the estuary. Reference is made to the mentioned papers for the analytical solutions of ( 6.28 ) and ( 6.29 ) and ( 6.30 ) or ( 6.31 ). All of these solutions are composed of classical trigonometric and exponential functions and are dependent on the type of estuary.
6.4. Response surface models The response surface methodology is an empirical model that consists of a group of mathematical and statistical techniques used in the development of an adequate functional relationship between a response of interest y and a number of associated controls or inputs xi. In general such a relationship is unknown but it can be approximated by a low-degree polynomial of the form (Khuri & Mukhopadhyay, 2010):
y = f (x ) ⋅ β + ε
( 6.34 )
With x = (x1, x2, …, xk) the controls or inputs, f(x) a vector function that consists of powers and crossproducts of powers of elements of x, β a vector with unknown constant coefficients and ε a random experimental error assumed to have a zero mean. The two commonly used models in the response surface methodology are the first degree model: k
y = β 0 + ∑ β i ⋅ xi + ε
( 6.35 )
i =1
And the second-degree model: k
k
y = β 0 + ∑ β i ⋅ xi + ∑∑ β ij ⋅ xi ⋅ x j + ∑ β ii ⋅ xi2 + ε i =1
i< j
i =1
( 6.36 )
IMDC et al. (2005) applied a first-degree surface response model for the tidal part of the river Scheldt to model maximum water levels and flood volumes. It appeared that the maximum water level at two locations of interest could be approximated sufficiently by considering only three control inputs: the downstream water level H, the tidal amplitude A and the wind speed V:
H x = θ 0 + θ H ⋅ H + θ A ⋅ A + θ V ⋅ V + θ AV ⋅ AV + ε Hx
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
( 6.37 )
73
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Where θi are the constant parameters of the model. An interaction term AV was introduced to account for the influence of the wind speed on the tidal amplitude. The magnitude of θAV was small compared to the other coefficients, so this last term could just as well be neglected. A similar first-degree model was obtained for the total flood volume:
g (V x ) = α 0 + α H ⋅ H + α A ⋅ A + α V ⋅ V + α AV ⋅ AV + ε Hx
( 6.38 )
With αi the coefficients of the model and g(Vx) a transformation that can be linear or logarithmic. Again, the interaction term AV is of minor importance compared to the other three controls. Equations ( 6.37 ) and ( 6.38 ) were calibrated with the maximum values of water levels and flood volumes of a number of hydrodynamic simulations for different scenarios. The main drawback of the surface response approach is that this particular model is only capable of predicting these maximum values. This also implies that the model is not suited for continuous simulations, since the model does not have a time component to allow for the travel time of the tidal wave. The accuracy in continuous simulations is therefore rather poor, with deviations between the hydrodynamic model and the surface response model up to one meter for the lower water levels. A similar research was conducted by Bárcena et al. (2012) for the river Saja in northern Spain (region of Cantabria), which has a shallow mesotidal estuary with a length of 5.5 km and a mean width of 150 m. A second order surface response model was calibrated to approximate the water levels and current velocities at different locations in the estuary. The approximating function used in this model takes a quadratic form, because of the non-linearity of the main processes in an estuary and because of the adaptability of this form:
z = a ⋅ Q 2 + b ⋅ A2 + c ⋅ Q ⋅ A + d ⋅ A + e ⋅ Q + f
( 6.39 )
With z either the water level or the flow velocity, Q the upstream river inflow, A the astronomic tidal amplitude and a to f constant coefficients, determined by means of a least square method. The model of Bárcena et al. (2012) is used as a continuous model and was calibrated with the results of a twodimensional hydrodynamic model, which on his turn was calibrated with measurements.
6.5. Conclusion Four simplified models with a limited computation time were discussed in this chapter. Only two of them (the bidirectional Muskingum method and the CALAM – model) can actually be classified as conceptual models, since they are based on the reservoir concept. The other two types of models (idealized models and the surface response methodology) are empirical or black-box models, but were also discussed since they have been applied before to the Western Scheldt estuary and have a rather simple calculation scheme and therefore need limited calculation times. The starting point of the bidirectional Muskingum water-stage method, dividing the total storage in prism storage and wedge storage and linking this with the water levels, looks promising, since it allows to calculate tidal propagation if a downstream (sea) level is known. However, the results indicate that the differences between observed and predicted water levels can become quite large, which probably limits the applicability of this method. Idealized models provide an analytical solution for the propagation of a tidal wave, when a sine or cosine function with constant amplitude is applied as the boundary condition at the downstream side. This should allow a fast calculation of water levels at different locations in an estuary or tidal river. The assumption of constant tidal amplitude implies that the model does not include the influences of neap and spring tides and storm surges on the water levels. Another drawback of idealized models is that they don’t explicitly model water levels in floodplains, tributaries and controlled holding basins. Unlike idealized models, the surface response models can handle tidal oscillations with varying amplitude, since the water level or tidal amplitude can be used as an input for the approximating function. They are, however, black box models with a lack of physical background and the accompanying drawbacks like the uncertainty that arises when the inputs go outside the calibration range or the difficulties in identifying the Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
74
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
model structure, with all its relevant inputs. Furthermore, the surface response models have no time component and are also not able to explicitly model water levels in floodplains. Because of all the drawbacks that come with the other model types, the CALAM – model looks the most suited for simulating water levels in estuaries and tidal rivers, with conceptual models. This model can handle both upstream discharge boundaries and downstream water level boundaries, that may be timedependent. The use of the reservoir concept also allows to model water levels at different locations, link it with the previously discussed reservoir routing technique and extend the model with a water quality component. In the end, the main factor influencing the model choice will be the specific aim of the model. Two extremes can be considered: knowing the variation of a water level at a certain location (for example at the confluence of the rivers Dender and Scheldt) or a detailed study of the water levels at different locations along the tidal river (for example for flood control). Idealized models and black box models like the surface response model seem sufficient for the first application, whereas the second application will most likely need a reservoir-type model.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
75
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
7. References Achleitner, S., Möderl, M., & Rauch, W. (2007a). CITY DRAIN © – An open source approach for simulation of integrated urban drainage systems. Environmental Modelling & Software, 22(8), 1184–1195. doi:10.1016/j.envsoft.2006.06.013 Achleitner, S., & Rauch, W. (2007b). CITY DRAIN 2 . 0 – an open source Matlab / Simulink library for integrated simulation of urban drainage systems. Ahlman, S., & Svensson, G. (2005). SEWSYS - a tool for simulation of substance flows in urban sewer systems (p. 48). Gothenburg, Sweden. Andrews, F. T., Croke, B. F. W., & Jeanes, K. W. (2010). Robust estimation of the total unit hydrograph. International Congress on Environmental Modelling and Software Modelling for Environment’s Sake (p. 8). Bárcena, J. F., García, A., García, J., Álvarez, C., & Revilla, J. a. (2012). Surface analysis of free surface and velocity to changes in river flow and tidal amplitude on a shallow mesotidal estuary: An application in Suances Estuary (Nothern Spain). Journal of Hydrology, 420-421, 301–318. doi:10.1016/j.jhydrol.2011.12.021 Bastin, G., & Moens, L. (2000). HYDROMAX : A Real-Time Application for River Flow Forecasting Hydromax : A Real-Time Application for River Flow Forecasting . In International Symposium on Flood Defence (Vol. 1, p. 9). Bastin, G., Moens, L., & Dierickx, P. (2009). On-line river flow forecasting with “ Hydromax ” : successes and challenges after twelve years of experience. In 15th IFAC Symposium on System Identification (pp. 1774–1779). Bentura, P. L. F., & Michel, C. (1997). Flood routing in a wide channel with a quadratic lag-and-route method. Hydrological Sciences Journal, 42(2), 169–190. Berlamont, J. (2002). Hydraulica (p. 400). Leuven: Acco. Berlamont, J. (2004). Rioleringen (p. 420). Leuven: Uitgeverij Acco. (ISBN 90-334-3891-7) Beven, K. (1989). Changing ideas in hydrology - the case of physically-based models. Journal of Hydrology, 105, 157–172. Beven, K. J. (2001). Rainfall-runoff modelling: The Pimer (pp. 1–360). Chichester: John Wiley & sons, LTD. Beven, K., Young, P., Leedal, D., & Romanowicz, R. (2009). Computationally efficient flood water level prediction (with uncertainty). In Samuels (Ed.), Flood Risk Management: Research and Practice (pp. 281–289). London: Taylor & Francis Group. Bishop, E., & Phelps, R. R. (1961). The support functionals of a convex set. Proceedings of the seventh symposium in pure mathematics of the American mathematical society (pp. 27–36). Bujon, G. (1988). Prévision des débits et des flux polluants transités par les réseaux d ’égouts par temps de pluie. Le modèle FLUPOL. La Houille Blanche, 1(1), 11–23. doi:http://dx.doi.org/10.1051/lhb/1988001 Borah, D. K. (2011). Hydrologic procedures of storm event watershed models: a comprehensive review and comparison. Hydrological Processes, 25(22), 3472–3489. doi:10.1002/hyp.8075 Boughton,W.C., 1993. A hydrograph-based model for estimating the water yield of ungauged catchment. In: Proceedings: Hydrology and Water Resources Symposium, Newcastle, June 30–July 2, 1993. The Institution of Engineers, Australia, pp. 317–324.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
76
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Cabus, P., De Jongh, I., & Cauwenberghs, K. (2007). Hydrologisch modelleren op onbemeten stroomgebieden : Een eenvoudig modelconcept op basis van PDM en waargenomen maximale afvoercoëfficiënten. Tijdschrift Water, 32, 55–59. (in Dutch) Calabro, P. S. (2001). Cosmoss: conceptual simplified model for sewer system simulation. A new model for urban runoff€ quality. Urban Water, 3(1), 33–42. Calabrò, P. S., & Viviani, G. (2006). Simulation of the operation of detention tanks. Water research, 40(1), 83–90. doi:10.1016/j.watres.2005.10.025 Camacho, L. ., & Lees, M. . (1999). Multilinear discrete lag-cascade model for channel routing. Journal of Hydrology, 226(1-2), 30–47. doi:10.1016/S0022-1694(99)00162-6 Capodaglio, A. G. (1994). Transfer function modelling of urban drainage systems, and potential uses in real-time control applications. Water Science & Technology, 29(1-2), 409–417. Chapra, S. C. (2012). Applied Numerical Methods with MATLAB for Engineers and Scientists (3rd ed., p. 672). New York (N.Y.): MacGraw-Hill. Chow, V. T., Maidment, D. R., & Mays, L. W. (1988). Applied Hydrology (p. 572). New York (N.Y.): MacGraw-Hill. Cluckie, I. D., Lane, A., & Yuan, J. (1999). Modelling large urban drainage systems in real time. Water Science and Technology, 39(4), 21–28. Coutu, S., Del Giudice, D., Rossi, L., & Barry, D. a. (2012). Parsimonious hydrological modeling of urban sewer and river catchments. Journal of Hydrology, 464-465, 477–484. doi:10.1016/j.jhydrol. 2012.07.039 Cunge, J. A. (1969). On The Subject Of A Flood Propagation Computation Method (Muskingum Method). Journal of Hydraulic Research, 7(2), 205–230. Cuppens, A. (2006) Implementation of a conceptual sewer modelling software on the Matlab-Simulink platform for the catchment of Herent. MSc-dissertation. Vrije Universiteit Brussel / Katholieke Universiteit Leuven. Decloedt, C., Willems, P. (2013). Methods and experiences in radar based fine scale rainfall estimation, 90 pp. Leuven: KU Leuven, Hydraulics division. De Keyser, W., Gevaert, V., Verdonck, F., De Baets, B., & Benedetti, L. (2010). An emission time series generator for pollutant release modelling in urban areas. Environmental Modelling & Software, 25(4), 554–561. DHI. (2011). Mike 11, A modelling system for rivers and channels, reference manual. Dooge, J. C. I. (1959). A general theory of the unit hydrograph. Journal of Geophysical Research, 64(2), 241–255. Erbe, V., Risholt, L. ., Schilling, W., & Londong, J. (2002). Integrated modelling for analysis and optimisation of wastewater systems – the Odenthal case. Urban Water, 4(1), 63–71. doi:10.1016/S14620758(01)00060-7 Euler, G. (1983). Ein hydrologisches Näherunsverfahren für die Berechnung des Wellenablaufs in teilgefüllten Kreisrohren. Wasser und Boden, 35(2), 85–88. Fenton, J. D. (1992). Reservoir doi:10.1080/02626669209492584
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
routing.
Hydrological
WL2014R00_131_1
Sciences
Journal,
37(3),
233–246.
77
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Fiorentini, M., & Orlandini, S. (2013). Robust numerical solution of the reservoir routing equation. Advances in Water Resources, 59, 123–132. doi:10.1016/j.advwatres.2013.05.013 Franchini, M., & Lamberti, P. (1994). A flood routing Muskingum type simulation and forecasting model based on level data alone. Water resources research, 30, 2183–2196. doi:10.1029/94WR00536 Fread, D. L. (1983). A unified Coefficient Routing Model. Presented at invited seminars at: Pennsylvania State University, February 28 - March 4,1983, 15 p. Fread, D. L., & Hsu, K. S. (1993). Applicability of Two Simplified Flood Routing Methods: Level-Pool and Muskingum-Cunge. In H. W. Shen, S. T. Su, & F. Wen (Eds.), ASCE National Hydraulic Engineering Conference (pp. 1564–1568). San Francisco, CA. Freni, G., Maglionico, M., & Di Federico, V. (2012). State of the art in Urban Drainage Modelling, CARE-S (p. 415). Friedrichs, C. T., & Aubrey, D. G. (1994). Tidal propagation in strongly convergent channels. Journal of Geophysical Research, 99(C2), 3321–3336. doi:10.1029/93JC03219 Gelormino, M., & Ricker, N. (1994). Model-predictive control of a combined sewer system. International Journal of Control, 59(3), 793–816. Ghiaseddin, N. (1986) An environment for development of decision support systems. Decis. Support Syst. 2(3), 195–212. doi:10.1016/0167-9236(86)90028-X Hibma, A., Schuttelaars, H. M., & Wang, Z. B. (2003). Comparison of longitudinal equilibrium profiles of estuaries in idealized and process-based models. Ocean Dynamics, 53(3), 252–269. doi:10.1007/s10236-003-0046-7 HIC (2014). Hydrologisch Informatiecentrum. www.waterstanden.be Hughes, D. A., & Murrell, H. C. (1986). Non-linear runoff routing - a comparison of solution methods. Journal of Hydrology, 85, 339–347. IMDC, HKV, & Nv, P. (2005). Kustverdediging - Sigmaplan. Studie naar het werken met overstromingsrisico’s: deelopdracht3. Overstromingsrisico Schelde. Toepassing van de responsoppervlakmethode. (p. 100). ITWH (2012). KOSIM 7.4 Anwenderhandbuch. Insitut für technisch-wissenschaftliche Hydrologie GmbH, Hannover, Germany. Kamradt, B. (2008). Comparison of sewer modelling approaches within IUWS modelling. Technische Universität Darmstadt. Khuri, A. I., & Mukhopadhyay, S. (2010). Response surface methodology. Wiley Interdisciplinary Reviews: Computational Statistics, 2(2), 128–149. doi:10.1002/wics.73 Kim, D. H. (2011). A new nonlinear hydrologic river routing model. PhD-dissertation. Georgia institute of technology. Knight, D., & Shamseldin, A. (2006). River Basin Modelling for Flood Risk Mitigation. (D. Knight & A. Shamseldin, Eds.) (p. 608). London, UK: Taylor & Francis Group plc. Kožar, I., Lozzi-Kožar, D., & Jeričević, Ž. (2010). A note on the reservoir routing problem. European Journal of Mechanics - B/Fluids, 29(6), 522–533. doi:10.1016/j.euromechflu.2010.06.005 Kroll, S., Blumensaat, F., Dirckx, G., Thoeye, C., Gueldre, D., & Steene, V. De. (2007). Assessment of CSO activity using simple volume balancing volume. In NOVATECH 2007 - 6th International Conference on Sustainable Techniques and Strategies in Urban Water Management (pp. 1639–1646). Lyon, France. Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
78
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Kroll, S., Wambecq, T., Weemaes, M., & Willems, P. (n.d.). Semi-automated model buildup and calibration of conceptual models for backwater influenced sewer systems. [in Press] Langergraber, G., Alex, J., Weissenbacher, N., Woerner, D., Ahnert, M., Frehmann, T., Halft, N., et al. (2008). Generation of diurnal variation for influent data for dynamic simulation. Water Science and Technology, 57(9), 1483–1486. doi:10.2166/wst.2008.228 Lanzoni, S., & Seminara, G. (1998). On tide propagation in convergent estuaries. Journal of Geophysical Research, 103(C13), 30793–30812. Laurenson, E.M., Mein, R.G. (1997). RORB – Version 4, Runoff routing Program – User Manual. Department of Civil Engineering, Monash University/Montech Pty. Ltd. Lees, M. J. (2000). Data-based mechanistic modelling and forecasting of hydrological systems. Journal of Hydroinformatics, 2(1), 15–34. Lekkas, D. F., Imrie, C. E., & Lees, M. J. (2001). Improved non-linear transfer function and neural network methods of flow routing for real-time forecasting. Journal of Hydroinformatics, 3(1995), 153–164. Linsley, R. K., Kohler, M. A., & Paulhus, J. L. H. (1975). Hydrology for engineers (2 nd. ed.). New York: McGraw Hill Book Company. MacArthur, R., & De Vries, J. J. (1993). Introduction and Application of Kinematic Wave Routing Techniques Using HEC-1 (p. 68). Davis, CA. Madsen, H. (2000). Automatic calibration of a conceptual rainfall–runoff model using multiple objectives. Journal of Hydrology, 235(3-4), 276–288. doi:10.1016/S0022-1694(00)00279-1 Mannina, G., Freni, G., & Viviani, G. (2004). Modelling the Integrated Urban Drainage Systems. In L. Bertrand-Krajewski, M. Almeida, J. Matos, & S. Abdul-Talib (Eds.), Sewer Networks and Processes within Urban Water Systems (WEMSno.) (pp. 3–12). London, UK: IWA Publishing. Meert, P., Wolfs, V., Villazon Gomez, M. F., & Willems, P. (2012). Conceptual Model Developer (CMD) for IWRS modeled river systems. Leuven: KU Leuven, Hydraulics division. Meert, P., & Willems, P. (2013). Evaluatie van berekeningsmethoden voor het bepalen van de benodigde buffercapaciteit van kleinschalige opvangsystemen in het kader van erosiebestrijding. Eindrapport, KU Leuven Afdeling Hydraulica, voor Vlaamse Overheid – Departement LNE – Afdeling Land en Bodembescherming, Ondergrond, Natuurlijke Rijkdommen (p. 246). Mehmood, K. (1995). Studies on sewer flow synthesis with special attention to storm overflows. PhDdissertation. University of Liverpool. Mockus, V., & Hjelmfelt, A. T. (2004). Estimation of Direct Runoff from Storm Rainfall. National Engineering Handbook, Part 630 Hydrology (p. 79). Mockus, V., & Styner, W. (1972). Chapter 17. Flood Routing. National engineering handbook - Section 4:Hydrology. Moore, R. J., & Bell, V. a. (2001). Comparison of Rainfall-Runoff models for flood forecasting. Part 1: Literature review of models. (p. 105). Brisol, UK, Environment Agency. Moore, R. J., & Bell, V. a. (2002). Incorporation of groundwater losses and well level data in rainfall-runoff models illustrated using the PDM. Hydrology and Earth System Sciences, 6(1), 25–38. doi:10.5194/hess-6-25-2002 Moore, R. J. (2007). The PDM rainfall-runoff model. Hydrology and Earth System Sciences, 11(1), 483–499. doi:10.5194/hess-11-483-2007
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
79
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Motiee, H., Chocat, B., & Blanpain, O. (1997). A storage model for the simulation of the hydraulic behaviour of drainage networks. Water Science and Technology, 36(8-9), 57–63. doi:10.1016/S02731223(97)00618-5 Muschalla, D., Reussner, F., Schneider, S., & Ostrowski, M. (2006). Dokumentation des Scmutzfrachtmodells SMUSI Version 5.0. Institut für Wasserbau und Wasswewirtschaft, Technische Universität Darmstadt, Germany MUSIC Development Team (2012). MUSIC user guide. Version 5.1. Cooperative Research Centre for Catchment Hydrology, Melbourne, Australia. Nash, J. E. (1960). A unit hydrograph study, with particular reference to British catchments. ICE Proceedings, Volume 17, Issue 3 (pp. 249–282). Nash, J.E., Barsi, B.I., (1983). A hybrid model for flow forecasting on large catchments. Journal of Hydrology 65 (1-3), 125–137. Nathan, R.J., McMahon, T.A., 1990. Evaluation of automated techniques for base flow and recession analyses. Water Resources Research 26 (7), 1465–1473. Nielsen, S.A. and Hansen, E. (1973), Numerical simulation of the rainfall runoff process on a daily basis. Nordic Hydrology, 4, 171-190. Ocampo Martinez, C. A. (2007). Model Predictive Control of Complex Systems including Fault Tolerance Capabilities : Application to Sewer Networks. Technical Univeristy of Catalonia. Omni Environmental. (2007) The Non-Tidal Passaic River Basin Nutrient TMDL Study Phase II Watershed Model and TMDL Calculations, 187. Princeton, NJ. Palmstrom, N., Walker, W.W., Jr., (1990). P8 urban catchment model: User’s guide. Program documentation, and evaluation of existing models, design concepts, and Hunt-Potowomut data inventory. The Narragansett Bay Project Report No. NBP-90-50. Perumal, M. (1992). Multilinear muskingum flood routing method. Journal of Hydrology, 133(3-4), 259–272. doi:10.1016/0022-1694(92)90258-W Perumal, M. (1994). Multilinear discrete cascade model for channel routing. Journal of Hydrology, 158(1-2), 135–150. Pilgrim, D. H. (1987). Australian rainfall and runoff: a guide to flood estimation. (D. H. Pilgrim, Ed.) (pp. 129– 149). Barton, ACT, Australia: The institution of engineers. Ponce, V. M., Lohani, A. K., & Scheyhing, C. (1996). Analytical verification of Muskingum-Cunge routing. Journal of Hydrology, 174, 235–241. Pritchard, D. W. (1967). What is an estuary: physical viewpoint. In G. H. Lauff (Ed.), Estuaries (pp. 3–5). Washington DC: American Association for the Advancement of Science. Rodríguez, J. P., McIntyre, N., Díaz-Granados, M., Achleitner, S., Hochedlinger, M., & Maksimović, Č. (2013). Generating time-series of dry weather loads to sewers. Environmental Modelling & Software, 43, 133–143. doi:10.1016/j.envsoft.2013.02.007 Romanowicz, R. J., Young, P. C., Beven, K. J., & Pappenberger, F. (2008). A data based mechanistic approach to nonlinear flood routing and adaptive flood level forecasting. Advances in Water Resources, 31(8), 1048–1056. doi:10.1016/j.advwatres.2008.04.015 Ruan, M., & Wiggers, J. B. M. (1998). A Conceptual CSO Emission Model: SEWSIM. Water Science and Technology, 37(1), 259–267.
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
80
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Saagi, R. (2013). Sewer System Model - Literature Review (SANITAS - Sustainable And Integrated Urban Water System Management) (p. 27). Sartor, J. (1999). Simulating the influence of backwater effects in sewer systems using hydrological model components. Water Science and Technology, 39(9), 145–152. doi:10.1016/S0273-1223(99)00227-9 Schramkowski, G. P., Schuttelaars, H. M., & Swart, H. E. De. (2004). Non-linear channel-shoal dynamics in long tidal embayments. Ocean Dynamics, 54(3-4), 399–407. doi:10.1007/s10236-003-0063-6 Schuttelaars, H. M., & Swart, H. E. De. (1997). An idealized long term morphodynamic model of a tidal embayment 1 Introduction. Eur. J. Mech., B/Fluids, 15, 55–80. Schuttelaars, H. M., & de Swart, H. E. (2000). Multiple morphodynamic equilibria in tidal embayments. Journal of Geophysical Research, 105(10), 24105–24118. Sherman, L. K., 1932, Streamflow from rainfall by the unit hydrograph method, Eng. News-Record, 108, 501-505. Smith, A. A. (1980). A generalized approach to kinematic flood routing. Journal of Hydrology, 45(1), 71–89. Solomatine, D. P., & Ostfeld, A. (2008). Data-driven modelling: some past experiences and new approaches. Journal of Hydroinformatics, 10(1), 3. doi:10.2166/hydro.2008.015 Solvi, A. (2006). Modelling the sewer-treatment-urban river system in view of the EU Water Framework Directive. PhD-dissertation. Ghent University. Szolgay, J., Danacova, M., & Papankova, Z. (2006). Case study of multilinear flood routing using empirical relationships between the flood wave speed and the discharge. Slovak Journal of Civil Engineering, 1, 1–9. Tan, B.Q., O’Connor, K.M., 1996. Application of an empirical infiltration equation in the SMAR conceptual model. Journal of Hydrology 185 (1–4), 275–295. Theurer, F. D., Comer, G. H., & Richardson, H. H. (1981). The modified Attenuation-Kinematic routing model.pdf. Proceedings of the International Symposium on Rainfall-Runoff Modelling (pp. 553–564). Technische Universität Hamburg – Harburg (2006). The method of Kalinin-Miljukov. Available at: http://daad.wb.tu-harburg.de/fileadmin/BackUsersResources/Hydrology/8HydrologicalModelling/additional_info/Method_of_Kalinin.pdf Todeschini, I., Toffolon, M., & Tubino, M. (2008). Long-term morphological evolution of funnel-shape tidedominated estuaries. Journal of Geophysical Research, 113(C5). doi:10.1029/2007JC004094.To Urban Drainage and Flood Control District. (2008). Colorado Urban Hydrograph Procedure: Excel-based Computer program (User Manual). Urbanus, J. F. X., & Vreeburg, J. H. G. (1987). Calamiteitenmodel Noordelijk Deltabekken. Technische Universiteit Delft. USBR (United States Bureau of Reclamation) (1988). Colorado River Simulation System Documentation, Denver. United States Department of Agriculture, National Resources Conservation Service (1965). TR-20, Computer Program for Project Formulation – Hydrology, May 1965 United States Department of Agriculture, Soil Conservation Service (1986). TR-20 Hydrology Technical Note 2, “Chapter 1: The Modified Att-Kin Routing Model” (p. 46).
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
81
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
U. S. Army Corps of Engineers (1994). Engineering and design - Flood-Runoff Analysis (p. 214). Washington DC. Vaes, G. and Berlamont, J. (1995), Methodologieën voor de bepaling van de overstortfrequentie met behulp van beperkte neerslaginvoer, Leuven: KU Leuven. Vaes, G. (1999). The influence of rainfall and model simplification on combined sewer system design. PhDdissertation. University of Leuven. Vandekerckhove, L., Swerts, M., Leyman, N., Mennens, K., Neven, H., Desmet, J., & De Vrieze, M. (2001). Code van goede praktijk voor het opmaken van een gemeentelijk erosiebestrijdingsplan (p. 83). Vanrolleghem, P., Benedetti, L. & Meirlaen, J. (2005) Modelling and real-time control of the integrated urban wastewater system. Environ. Model. Softw. 20(4), 427–442. doi:10.1016/j.envsoft.2004.02.004 Vanrolleghem, P. A., Kamradt, B., Solvi, A., & Muschalla, D. (2009). Making the best of two Hydrological Flow Routing Models: Nonlinear Outflow-Volume Relationships and Backwater Effects Model. Proceedings 8th international conference on Xenobiotics in the Urban Water Cycle. Vélez, C., Gimeno, P., Sheresta, N., Tolesa, O., Griensven, A. van, Bauwens, W. & Pereira, F. (2013) Surrogate Models in Decision Support Systems for Water Quality Management. El Riesgo en la Gest. del Agua. Conf. interancional AGUA 2013, 10. Cali. Villazon Gomez, M. F. (2011). Modelling and conceptualization of hydrology and river hydraulics in flood conditions, for Belgian and Bolivian basins. PhD-dissertation. KU Leuven. Villazon Gomez, M. F., Meert, P., Wolfs, V., & Willems, P. (2012). Conceptual Model Developer ( CMD ) Manual. Leuven: KU Leuven, Hydraulics division. (p. 76 ) VMM (2014). Overstromingsvoorspeller, http://www.overstromingsvoorspeller.be Voll, J. (2005). Bestimmung des Einflusses von veränderten hydrologischen Wellenablaufberechnungen auf die Simulationsergebnisse des NA-Modells Fränkische Saale. (in German) Universität der Bundeswehr München. Wagener, T., Wheater, H. S., & Gupta, H. V. (2004). Rainfall-Runoff Modelling in Gauged and Ungauged Catchments (p. 332). London, UK: Imperial College Press. Wallingford (2012a). InfoWorksTM Collection Systems (CS) Help Documentation (version 13.0). UK, Wallingford Software. Wallingford (2012b). InfoWorksTM River Systems (RS) Help Documentation (version 13.0). UK, Wallingford Software. Wery, B. (1990). Identification des systèmes hydrologiques - Application à la prévision des crues. PhD-dissertation. Université de Louvain. (in French) Wikipedia. (2013) Decision support system. in Wikipedia. http://en.wikipedia.org/wiki/Decision_support_system
Retrieved
October
2,
2013,
from
Willems, P. (2000). Probabilistic immission modelling of receiving surface waters. PhD-dissertation. KU Leuven. Willems, P. (2009). A time series tool to support the multi-criteria performance evaluation of rainfall-runoff models. Environmental Modelling & Software, 24(3), 311–321. doi:10.1016/j.envsoft.2008.09.005 Willems, P. (2013). Waterloopmodellering (p. 268). Leuven & Den Haag: Acco. ISBN 978-90-334-9296-9, D/2013/0543, NUR 955. (in Dutch) Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
82
Development of conceptual models for an integrated catchment management: Subreport 1 – Literature review of conceptual modelstructures
Willems, P. (2014). Parsimonious rainfall-runoff model construction supported by time series processing and validation of hydrological extremes - part 1: Step-wise model-structure identification and calibration approach. Journal of Hydrology. [in press] Willems, P., Mora Serrano, D., Vansteenkiste, T., Teferi Taye, M., & Van Steenbergen, N. (2014). Parsimonious rainfall-runoff model construction supported by time series processing and validation of hydrological extremes - part 2: Intercomparison of models and calibration approaches. Journal of Hydrology. [in press] Wolfs, V., Villazon, M. F., & Willems, P. (2013). Development of a semi-automated model identification and calibration tool for conceptual modelling of sewer systems. Water science and technology, 68(1), 167– 75. doi:10.2166/wst.2013.237 Yevdjevich, V. M. (1959). Analytical Integration of the Differential Equation for Water Storage. Journal of Research of the National Bureau of Standards - B. Mathematics and Mathematical Physics, 63(1), 43–52. Young, P. C. (1998). Data-based mechanistic modelling of environmental, ecological, economic and engineering systems. Environmental Modelling & Software, 13(1), 105–122. Young, P. (2001). The Identification and Estimation of Nonlinear Stochastic Systems. In A. I. Mees (Ed.), Nonlinear Dynamics and Statistics (pp. 127–168). Boston: Birkhauser. Yuan, J. (1994). Hydrological modelling with weather radar data in urban drainage systems. University of Salford. Zoppou, C. (2001). Review of urban storm water models. Environmental Modelling & Software, 16(3), 195– 231. doi:10.1016/S1364-8152(00)00084-0 Zug, M., & Phan, L. (1999). Horus, un modèle conceptuel de simulation de la pollution en réseau d’assainissement - structure et validation. Revue des sciences de l’eau / Journal of Water Science, 12(4), 643–660. doi:10.7202/705370ar
Final version F-WL-PP10-1 Version 04 RELEASED AS FROM: 12/11/2012
WL2014R00_131_1
83
Waterbouwkundig Laboratorium Flanders Hydraulics Research
Berchemlei 115 B-2140 Antwerp Tel. +32 (0)3 224 60 35 Fax +32 (0)3 224 60 36 E-mail:
[email protected] www.waterbouwkundiglaboratorium.be
