Česká zemědělská univerzita v Praze Fakulta životního prostředí

Česká zemědělská univerzita v Praze Fakulta životního prostředí Katedra vodního hospodářství a environmentálního modelování

Aplikace srážko-odtokového modelu Boussmo Diplomová práce

Vedoucí diplomové práce: Ing. Michal Kuráž Diplomant: Michal Steinhart

2010

Prohlášení Prohlašuji, že jsem tuto diplomovou práci vypracoval samostatně pod vedením Ing. Michala Kuráže. Uvedl jsem všechny literární prameny, ze kterých jsem čerpal.

V Praze 30.4.2010

............................. Michal Steinhart

Poděkování Mé poděkování patří Ing. Michalu Kuráži za trpělivé konzultování a kvalitní vedení při vypracování této práce. Dále děkuji Ing. Jiřímu Pavláskovi za poskytnutí dat a cenných rad. Děkuji všem na katedře vodního hospodářství a environmentálního modelování za rady a ochotu. Mé hlavní poděkovaní patří mým rodičům a prarodičům, bez jejichž všestranné podpory by tato práce nikdy nevznikla.

Diploma thesis: Application of rainfall-runoff model Boussmo Abstract: Boussmo is a conceptual hydrological model based on a numerical solution of the Boussinesque equation for the subsurface flow and the kinematic wave equation for the surface flow. The model neglects evapotranspiration and an unsaturated zone and thus is useful for rainy seassons in tropical areas. The aim of my work is an estimation of a saturation on the experimental catchment Modrava 2 and an application of the Boussmo model for this catchment. For the estimation of the unsaturated zone an the empirical concept of Antecedent Precipitation Index was selected. The soil water content is expressed from a relation of a rainfall volume before an upswing limb of a hydrogram and API. A part of this work is a model calibration and its validation.

Key words: Boussinesque equation, Antecedent precipitation index (API), Calibration, Validation.

Abstrakt: Boussmo je konceptuální hydrologický model založený na numerickém řešení Boussinesquovi rovnice pro podpovrchový odtok a rovnice kinematické vlny pro povrchový odtok. Model zanedbává evapotranspiraci a nenasycenou zónu a hodí se tak pro období dešťů v tropických oblastech. Cílem této práce je odhad nasycenosti experimentálního povodí Modrava 2 a aplikace Boussma na toto povodí. Pro odhad nenasycené zóny byl vybrán empirický koncept předchozího srážkového indexu. Půdní vlhkost je vyjádřena ze vztahu objemu srážky před vzestupnou větví hydrogramu a API. Část této práce je kalibrace a validace Boussma.

Klíčová slova: Boussinesqueova rovnice, Předchozí srážkový index (API), Kalibrace, Validace.

1. 2. 3.

ÚVOD...................................................................................................................- 9 CÍLE PRÁCE......................................................................................................- 10 REŠERŠE ...........................................................................................................- 11 3.1 Úvod do srážko-odtokových modelů a jejich rozdělení .............................- 11 3.1.1 Příprava matematického modelu ........................................................- 11 3.1.2 Rozdělení srážko-odtokových modelů .................................................- 13 3.1.3 Kybernetické modely (“black box“) ...................................................- 16 3.1.4 Koncepční modely (“grey box“).........................................................- 17 3.1.5 Hydrodynamické modely (“white box“) .............................................- 19 3.2 Popis Boussma............................................................................................- 22 3.2.1. Matematický model .............................................................................- 23 3.2.1.1. Boussinesqueova rovnice(BR) ....................................................- 23 3.2.1.1.1 BR ustáleného proudění na nakloněné nepropustné rovině ....- 26 3.2.2.1.2 BR neustáleného proudění na nakloněné nepropustné rovině- 28 3.2.1.2 Rovnice kinematické vlny...........................................................- 29 3.2.2 Procedurální model ............................................................................- 30 3.3 Předchozí srážkový index. ..........................................................................- 32 3.3.1 Porovnání modelů APIc a SAC-SMA..................................................- 40 3.4 Kalibrace a optimalizace parametrů ...........................................................- 40 3.4.1 Kalibrace a validace...........................................................................- 40 3.4.2 Automatická optimalizace parametrů .................................................- 44 3.4.3 Objektivní funkce ................................................................................- 47 4. METODIKA .......................................................................................................- 50 4.1 Charakteristika povodí Modravy 2 a použitých dat....................................- 50 4.1.1 Popis povodí .......................................................................................- 50 4.1.2 Základní geometrické charakteristiky.................................................- 52 4.1.3 Geologické poměry povodí .................................................................- 54 4.1.4 Půdní a vegetační kryt povodí ............................................................- 56 4.1.5 Posouzení srážkových dat ...................................................................- 56 4.2 Výběr srážko-odtokových událostí .............................................................- 60 4.3 Stanovení API .............................................................................................- 61 4.4 Závislost počátečních ztrát na API ............................................................- 62 4.5 Příprava dat pro kalibraci modelu Boussmo...............................................- 64 4.6 Automatická optimalizace ..........................................................................- 65 5. VÝSLEDKY.......................................................................................................- 65 5.1 API a odhad počátečních ztrát ....................................................................- 65 5.2 Kalibrace.....................................................................................................- 70 5.3 Validace ......................................................................................................- 72 6. DISKUZE ...........................................................................................................- 79 7. ZÁVĚR ...............................................................................................................- 80 8. PŘEHLED LITERATURY A POUŽITÝCH ZDROJŮ ....................................- 81 9. PŘÍLOHY ...........................................................................................................- 84 -

Seznam příloh: Příloha č.1: Konfigurační soubor modelu Boussmo Příloha č.2: Porovnání výsledků modelu Sacramento s modelem APIc (Smith 2000) Příloha č.3: Thomsonův přeliv (Experimentální povodí Modrava) Příloha č.4: Lokalizace povodí Modrava 2 (Experimentální povodí Modrava) Příloha č.5: Průzkumné profily (Levý 2008) Příloha č.6: Skripty v programu R Příloha č.7: Regresní analýza pro různé typy API Příloha č.8: Průběh konvergence chyby Příloha č.9: Validace nevykazující dobrou shodu

1.

ÚVOD

Tato diplomová práce nepřímo navazuje na mou bakalářskou práci, v které jsem se zabýval hydrologickou studií vybraného povodí. Převážná část práce představovala odvození geometrických charakteristik z digitálních mapových podkladů a výpočet odtokových čísel CN pro odvození n-letých maximálních průtoků pomocí hydrologického modelu DesQ-MaxQ. V podstatě se jednalo o popis vnějších vstupů do modelu než o samotný model. V této práci je mou snahou se zaměřit na vnitřní strukturu a principy samotného modelu. Úvodní část práce je věnována rozdělení srážko-odtokových modelů do nejběžnějších skupin. Pro každou skupinu modelů je uveden nejpoužívanější či nejtypičtější zástupce. Celá tato část je směřována pro vymezení a charakterizování srážko-odtokového modelu Boussmo, kterému je věnována další část práce, týkající se zejména odvození rovnic, na kterých je model založen. Před aplikací modelu na povodí Modravy 2 je nutné nejprve provést odhad nenasycené zóny, zvolit pro tento účel vhodnou metodu a zapracovat ji jako komponentu pro odhad počátečních ztrát do modelu. Zbývající část práce je zaměřena na kalibraci a validaci takto upraveného modelu Boussmo a na posouzení jeho použitelnosti na experimentálním povodí Modrava 2.

-9-

2.

CÍLE PRÁCE Cílem mé práce je kalibrace a validace srážko-odtokového modelu Boussmo na

experimentálním povodí Modrava 2. Před samotnou kalibrací modelu je nezbytné učinit odhad stavu nenasycené zóny a stanovit pro model počáteční ztráty. Nutnou součástí práce je výběr vhodných srážkoodtokových epizod z poskytnutých dat.

- 10 -

3. 3.1

REŠERŠE Úvod do srážko-odtokových modelů a jejich rozdělení

3.1.1 Příprava matematického modelu Ačkoliv je každý model pouhým zjednodušením hydrologického či jiného procesu a jeho výsledné simulace budou vždy zatíženy nějakou chybou, je nepostradatelným nástrojem pro získání rámcové představy o chování sledovaného systému. Pod pojmem chování systému v hydrologii se rozumí komplexní reakce systému na vstup neboli transformační funkce povodí, podle níž se transformuje efektivní déšť na povrchový odtok. Pro účel vymezení transformační funkce je nutné pomocí řady parametrů popsat fyzikálně geometrické vlastnosti systému jako jsou plocha, hydraulická drsnost, sklon, půdní vlastnosti, retence a mnoho dalších (Hrádek, Kuřík 2008; Kovář 1990). Obecně platí, že čím komplikovanější model, tím více parametrů používá a je třeba si dát pozor na přeparametrizování modelu, které vede k větším možným neurčitostem ve výsledku (Beven 2002). Matematické modely svou variabilitou a flexibilitou zcela zastoupily fyzikální modely, které se ve vodohospodářské praxi uplatňují už jen okrajově v některých speciálních odvětvích hydrauliky. Rozvoj matematických modelů nastal zcela přirozeně s rozvojem výpočetní techniky, která přináší kvalitnější měření dat a potřebný prostor pro časově náročné výpočty numerické matematiky, jejímž prostřednictvím se výrazně zpřesnili výsledky matematických modelů a mohlo tak dojít k jejich implementaci. Při sestavování matematického modelu je nutné mít v první fázi co nejpodrobnější představu o chování budoucího modelu tj. sestavit si percepční model, který je vždy do jisté míry ovlivněn subjektivní představou osoby, která daný model navrhuje. Druhou fází je přepis perceptuálního modelu do vhodných matematických rovnic. Tento přepis je vždy pouhým zjednodušením perceptuálního modelu. Nejčastěji tomuto účelu slouží rovnice kontinuity a hybnosti, které jsou často doplněny konstitučními vztahy nebo dalšími doplňkovými rovnicemi podle modelovaných procesů. Důležitou součásti je určení okrajových podmínek. Závěrečnou fází je poté model procedurální, který představuje naprogramování konceptuálního modelu. Výsledný algoritmus je také třeba přizpůsobit hardwarovému a softwarovému vybavení.

- 11 -

U takto připraveného simulačního prostředku je nutné kontrolovat chování jeho jednotlivých elementů. Tato průběžná kontrola se provádí numerickými testy a experimenty, kdy se porovnává známý průběh funkce nebo výsledné chování elementu s analyticky získanými nebo měřenými hodnotami (Ředinová 2004). Poté se model kalibruje a validuje. Jednotlivé hlavní kroky při sestavování modelu (obr. 1)

revize rovnic

úprava kódu

revize parametrů

Perceptuální Perceptuální model model

Konceptuální Konceptuální model model

Procedurální model

zvyšující se aproximace

revize předpokladů

Kalibrace modelu

Validace modelu

ne

Shoda? ano

Obrázek č.1. Základní schéma při sestavování modelu (Beven 2002)

- 12 -

3.1.2 Rozdělení srážko-odtokových modelů Rozdělení matematických hydrologických modelů slouží k základní orientaci v jejich uplatnění a smyslu použitelnosti. Ne vždy je snadné model stručně popsat a přesně vymezit interpretovatelnost jeho výstupu, obzvláště pokud se jedná o model, který pro účel svého řešení používá kombinovaných přístupů či rozdílných metod založených na rozdílných předpokladech, ať již matematických, fyzikálních či empirických. Z hlediska aplikace se modely rozdělují do dvou základních skupin. A to na modely: •

predikční (návrhové)

•

předpovědní (operativní)

Predikční slouží pro základy civilního inženýrství, pro stavbu vodních děl apod. (Dingman 2002) či jak uvádí Daňhelka (2003) jejich oblast využití je pro návrh, plánování případně pro konzultaci v oblasti vodního hospodářství. Předpovědní modely slouží pro oblast operativní hydrologie, kde model slouží k předpovědi odtoku z povodí na srážku. Dále je vhodné dělit modely dle způsobu schematizace topografie a to z hlediska rozčlenění vstupních a stavových veličin na: •

celistvé

•

distribuované

•

semi-distribuované

•

modely 1D, 2D, 3D

U celistvých jsou parametry v celém výpočetním procesu neměnné. Oproti tomu parametry distribuovaných modelů se mění v čase i v prostoru. Fyzikálně založené modely obsahují pouze data zjištěná měřením nebo odvozením a stává se u nich při dostatečném množství dat, že se nemusí kalibrovat (Daňhelka 2003). Zejména se jedná o hydrodynamické modely viz dále. U semi-distribuovaných se mění pouze některé parametry. Tvoří přechod mezi první a druhou skupinou modelů. Modely 1D interpolují hladinu mezi příčnými profily. Vyžívají se k modelování říčních systému a předpovídání povodní. 2D modely poskytují informace a vodním a rychlostním stavu.

- 13 -

používají se např. pro výpočty týkajících se plavebních kanálů. 3D modely se používají relativně málo, oproti 1D a 2D jsou schopné podat informaci o tom, zda je vtoku větší rychlost u dna či u hladiny (DHI 2000). Pro rozdělení modelů se používá i časové a prostorové hledisko a to na modely: •

kontinuální

•

diskrétní (epizodní)

•

regionální

•

lokální

Kontinuální modely využívají dlouhé časové řady srážek a dalších potřebných hydrometeorologických dat. Diskrétní modely využívají krátkou časovou řadu a nepřívalové srážky do nich nevstupují (Kovář 1990). Samozřejmě řada modelů může sloužit pro oba účely. Regionální se uplatňují na povodích o rozloze v řádu sto až tisíců km2. Za lokální modely se označují takové, které se používají na povodích o rozloze řádu desítek km2. Dle rozsahu výpočtů častí hydrologického cyklu můžeme dále dle Kováře (1990) dělit modely na: •

komplexní (snaha o popis celého hydrologického cyklu)

•

komponentní (popis pouze vybraných částí hydrologického cyklu)

Podle Dingmana je účelné dělit modely i podle způsobu výpočtu a to na:

Výpočty

•

0-dimenzionální

•

Analytické

•

Numerické

•

Hybridní 0-dimenzionálních

modelů

nejsou

založeny

na

formálním

souřadnicovém systému. Obvykle se používají u celistvých modelů. Analytické řešení počítá v souřadnicovém systému s diferenciálními rovnicemi, které se dají vyřešit analyticky. Při numerickém řešení jsou diferenciální rovnice řešeny metodou konečných diferencí, konečných prvků, konečných objemů a řadou dalších. Za hybridní řešení se považuje spojení 0-dimenzionálního řešení a formálního řešení pro model.

- 14 -

Za jedno ze základních rozdělení hydrologických modelů považuje Zeman (1994) hledisko kauzality: •

stochastické

•

deterministické

Jestliže kterákoliv proměnná vystupující v modelu je nahodilá, tj. má nějaké pravděpodobnostní rozdělení, jedná se potom o model stochastický (Kovář 1990). Parametry jsou tedy náhodně generovány a dva shodné soubory vstupních dat mohou dát rozdílné výsledky. Stochastické modely dále dělíme na pravděpodobnostní modely a na modely pro generování časových řad (Zeman 1994). Používají se při extrapolaci časových řad nebo hydrologických parametrů při zachování základních statistických charakteristik. Klasickým příkladem je model ARIMA (Daňhelka 2003) U deterministických modelů je každá proměnná reprezentována jednou hodnotou a jejich vztahy mezi sebou i k parametrům jsou pouze příčinné tedy deterministické. Účelem deterministických modelů v hydrologické aplikaci je popsat co možná nejpřesněji matematickými rovnicemi vztahy určité fyzikální představy, které jsou předmětem našeho zájmu. Čím je popis fyzikálních vztahů lepší tím je pochopitelně přesnější. V praxi však vyšší stupeň přesnosti matematického popisu klade náročnější požadavky na vstupní data. S ohledem na omezenou kvalitu i kvantitu pozorovaných proměnných a tím i odvozených parametrů se vyvinuli dvě hlavní větve deterministických modelů (Kovář 1990): •

hydrodynamické modely (“white box”)

•

hydrologické (parametrické) modely

Pro hydrologické modely jsou typické dva přístupy: •

kybernetický (“black box“)

•

koncepční (“grey box“)

Podrobné rozdělení deterministických modelů dle Kováře (1990) (obr. 2)

- 15 -

Obrázek č.2. Rozdělení deterministických modelů (Kovář 1990)

Ve vodohospodářské praxi se uplatňují ve všech oblastech. Obecně se dá říci, že jsou uživatelům přístupnější, neboť se svým vnitřním uspořádáním snaží přiblížit jednotlivým procesům hydrologického cyklu a jsou tak fyzikálně i matematicky srozumitelnější. Nemají požadavky na existenci extrémně dlouhých řad (Daňhelka 2003). 3.1.3

Kybernetické modely (“black box“) Tento typ modelu je zaměřen převážně na transformační funkci systému.

Ignoruje změny stavových veličin a není podstatná struktura systému. Využívá metod systémové analýzy z oboru kybernetiky ke zkoumání systému (Kovář 1990). Vhodný pro systémy s jednotným chováním a jednoduchou strukturou. Ke správnému fungování potřebuje i výstupní data, kvůli identifikaci funkce, která vystihuje chování systému. Jelikož jsou vztahy mezi vstupními a výstupními daty zpravidla pouze empirické vyžadují tyto modely častou rekalibraci. Princip black box modelu se používá jako komponenta větších modelů (Daňhelka 2003). Mezi modely typu black box můžeme zařadit například Nashův model, triangle (Beven 2002) oba založené na koncepci jednotkového hydrogramu či model akumulačního typu Doogův model nebo modely různých kombinací kaskád a fiktivního systému nádrží jakým je např. Tank model, který je bilančním modelem simulujícím hydrologický cyklus konečným počtem hydrologických nádrží. Schematická struktura povodí je vyjádřena uspořádáním těchto

- 16 -

nádrží. Toto pojetí zanedbává hybnostní a energetické vztahy a počítá pouze se vztahy kinematickými. Optimalizovanými veličinami jsou parametry přepadů (vstupů) a výpustí (výstupů). 3.1.4 Koncepční modely (“grey box“) Pro tyto modely je typické formulovat jednotlivé části hydrologického cyklu nebo cyklus jako celek matematickými vztahy. Jedná se o modely konceptuální, odrážející základní zákonitosti ve zjednodušené (koncepční) formě (Daňhelka 2003). U těchto modelů je snaha o co největší analogii mezi strukturou modelu a strukturou zkoumaného jevu. Tento přístup se vyhýbá prostorovým vztahům a omezuje se na předpoklad, že k prostorovým změnám veličin dochází pouze na reprezentativních bodech objektu. Díky této diskretizaci vede řešení často na obyčejné diferenciální rovnice, kde jedinou proměnnou je čas. Pro většinu koncepčních modelů je nutno v identifikační fázi jejich použití počítat s upřesňováním jejich parametrů některou optimalizačních technik. Dobrým příkladem koncepčního modelu je Stanfordský model, který jako první na světě aplikoval lineární kumulativní rozdělení hodnot některých parametrů kolem jejich průměrných hodnot na povodí. Model má 34 parametrů z nichž nejméně čtyři je nutno optimalizovat. Zbývající parametry je možno vyhodnotit z map, průzkumů a hlavně z měřených dat srážek, průtoků, a některých meteorologických dat (Kovář 1990) dalším příkladem je model BROOK90, který má mnoho parametrů s předurčenými hodnotami a pro obdržení rozumných výsledků je není třeba optimalizovat. Používá se pro všechny typy povrchu. Model stanovuje intercepci a transpiraci z jedné rostlinné vrstvy; půdní a sněhovou evaporaci, akumulaci a tání sněhu a samozřejmě povrchový a podpovrchový odtok (Dingman 2002). V současné době nejpoužívanějším koncepčním modelem je Sacramento, jehož schéma (obr.3). Sacramento Soil Moisture Accounting model (SAC-SMA) je srážkoodtokový model vyvíjený od poloviny 70. let národní meteorologickou službou (NWS) v USA jmenovitě Robertem Burnash s Larrym Ferralem. Je to konceptuální hydrologický model založený na parametrizaci charakteristik půdní vlhkosti. V modelu je aktivní vrstva půdy reprezentována dvěma zónami a to dolní (dlouhodobá zásoba jako např. půdní vlhkost a podzemní voda) a horní (časově krátká zásoba), obě zóny mají vodu vázanou, ovlivněnou adhezí a kohezí a jednu nebo více nádrží s vodou volnou, která není vázáná půdními částicemi a volně se pohybuje ve směru gravitace. V horní zóně srážka nejprve naplní nádrž s vodou vázanou, v které je zadržena a může být - 17 -

odstraněna pouze evaporací, kapacita této nádrže vyjadřuje množství srážek, které jsou nutné k vyplnění všech pórů v horní části půdního profilu. Jakmile se v horní zóně naplní zásobník s vodou vázanou dojde k postupnému plnění zásobníku pro vodu volnou. Poté co se obě horní nádrže naplní dochází k perkolaci do dolní vrstvy nebo se voda dále chová jako podpovrchový odtok. Každá srážka přesahující v horní zóně kapacitu vody vázané a volné vystupuje jako rychlá odezva povodí v podobě přímého odtoku. I dolní zóna obsahuje nádrže pro vodu vázanou a volnou, přičemž jakmile dojde k naplnění nádrže s vodou vázanou perkolací se zahájí plnění dvou nádrží s vodou volnou. Odtok z těchto dvou nádrží generuje krátkodobý a dlouhodobý základní průtok (Burnash, Ferral 1996). Jak uvádí Daňhelka (2003) předpokládá se, že odvodnění spodní zóny probíhá podle Darcyho zákona. Velikost toku lze stanovit jako součin hydraulické vodivosti a gradientu mechanické energie. V modelu Sacramento je součinitel vodivosti násoben zbytkovou volnou vodou. Tento předpoklad bohužel neumožňuje simulaci většího počtu typů odtoků, jak je možno pozorovat ve skutečné přírodě. Za předpokladu, že model obsahuje dva typy spodních zón s volnou vodou (primární zónu, která se velmi pomalu vyprazdňuje a poskytuje základní odtok pro dlouhodobé období a druhý typ, který podporuje odtok pro období s velice řídkými srážkami) je možná kombinací těchto zón, primární a sekundární, které se odvodňují nezávisle na Darcyho zákonu a umožňují aproximovat různé typy odtoků vyskytujících se v reálné přírodě (Burnash, Ferral 1996). V každém časovém kroku jsou vstupy a výstupy z různých zásobníků sčítány k určení celkového objemu. Tento model je vhodný i pro povodí o rozloze větší než 1000 km2.

- 18 -

Obrázek č. 3. Schéma modelu SAC-SMA (Daňhelka 2003).

3.1.5 Hydrodynamické modely (“white box“) Tyto modely jsou svým pojetím opakem kybernetických modelů. Jsou založeny na fyzikálním základě a více méně respektují principy zachování hmoty, hybnosti a zachování energie. Skutečná podstata systému je vyjádřena pomocí diferenciálních rovnic. Praktickou stránkou hydrodynamického modelu je algoritmus řešení těchto rovnic, převedených do algebraických lineárních rovnic. Zatímco struktura systému je u konceptuálních hydrologických modelů součástí modelu, u hydrodynamických je vložena přímo do základních rovnic. Pro sestavení a implementaci je nutné mít dle Kováře (1990) následující informace: dobře vymezené přírodní zákony podle kterých daný přírodní proces probíhá a je popsán formou parciálních diferenciálních rovnic (např. rovnice kontinuity a pohybové). Geometrický systém potřebný k diskretizaci diferenciálních rovnic do rovnic diferenčních. Numerické schéma, které umožní převedení výchozí rovnice do diferenčního tvaru.s využitím geometrického systému. Dále hodnoty hydrologických a hydraulických proměnných a parametrů. Nakonec je důležité správné určení počátečních a okrajových podmínek (Kovář 1990). Obecný odtokový model obvykle zahrnuje sub-modely tří dominantních procesů. První procesem je PRODUKCE efektivního deště z příčinného deště včetně vyčíslení příslušných ztrát. Druhým je TRANSFORMACE efektivního deště do povrchové a podpovrchového odtoku a posledním je proces PROPAGACE charakteristik odtoku v oblasti řešení času a prostoru. Procesy produkce a transformace jsou zpravidla modelovány pomocí konceptuálních modelů. Pro proces propagace se - 19 -

mnohem lépe uplatní hydrodynamický model. Odvození základních hydrodynamických rovnic provedl St. Venant (Kovář 1990). Jako zástupce těchto modelů jmenujme TOPMODEL, který jak uvádí Beven (2001) je založen na principu hydrologické podobnosti, která spočívá v podobnosti různých bodů na povodí pomocí jednoduché teorie o topografii a půdách. Základní myšlenkou je předpoklad, že body se stejným topografickým indexem budou mít v systému stejné chování. Model tedy počítá hodnoty distribuční funkce pouze pro reprezentativní body se stejnými hodnotami indexu. Tím snižuje délku výpočtu. Pro TOPMODEL je dynamika saturované zóny aproximována po sobě jdoucími ustálenými stavy na ploše α a hydraulický gradient saturované zóny je aproximován lokálním topografickým sklonem tan β . Pro výpočet topografického indexu slouží rovnice:  α  ln    tan β  a pokud se hodnota transmisivity T0 mění v prostoru, je index roven:

 α  ln    T0 tan β  Tento koncept se však nedá použít na povodí v oblastech se silným sezónním suchem. Jako další zástupce je model MIKE SHE vycházející z modelu SHE (Systéme Hydrologique Européen) z roku 1977. který je v součastné době považován za nejúplnější hydrodynamický přístup. Představuje vysoce integrovaný, fyzikálně založený distribuovaný modelovací systém. Tento systém, jehož schéma (obr. 4) umožňuje simulovat všechny fáze pevninského cyklu. Základním modulem je modul pro pohyb vody (dále WM–water movement module). K tomuto modulu se dají připojit další moduly, které simulují přídavné procesy jako jsou například advekce nebo disperze, biodegradaci, půdní eroze atd. Výhodou těchto modulů je, že mohou vystupovat jako samostatné jednotky či v interakci s ostatními. Charakteristiky povodí a vstupní data jsou zobrazena na vodorovném plánu ve výpočetní síti. Uvnitř každého elementu jsou popsány vertikální změny půdy a další hydrogeologických vlastností a to ve všech vodorovných vrstvách majících proměnlivou hloubku.

- 20 -

Povrchový odtok je schematizován dvojdimenzionální aproximací Saint Venantových rovnic, soustředný odtok pouze jednodimenzionálním aproximací těchto rovnic, pro pohyb vody nenasycenou zónou je použita Richardsonova rovnice, podzemní odtok řeší Boussinesqueho vztah. Síť je též možno použít v 3D nebo v kvazi3D nebo 2D (DHI 2000).

Obrázek č.4. Ilustrační schéma modelu Mike She (DHI 2000).

- 21 -

3.2

Popis Boussma Srážko-odtokový předpovědní model Boussmo [buzmo] navrhl Michal Kuráž

společně s Jiřím Pavláskem a naprogramoval Michal Kuráž. Název Boussmo vznikl spojením prvních pěti písmen ze slova 'bouss'inequova rovnice, na které je model převážně založen a z prvních dvou písmen ze slova 'Mo'drava, což je experimentální povodí, pro které byl navržen. Součastná verze modelu je určena pro epizodní modelování. Základ modelu tvoří numerické řešení boussinesqueovy rovnice pro podpovrchový odtok a rovnice kinematické vlny pro odtok povrchový. Rozlišení těchto dvou odtoků je řešeno pomocí dvou nádrží a půdních podmínek. Boussmo zatím nepočítá s evapotranspirací a zanedbává nenasycenou zónu. Po skončení výpočtu podá model informace o odtokovém koeficientu, o celkovém objemu srážek a o celkovém objemu podpovrchového a povrchového odtoku. Podpovrchový odtok je dále dělen na odtok ústící do koryta a na odtok, který se dostane pod úroveň koryta (obr. 5). Toto rozdělení je závislé na šířce koryta zadávané uživatelem (Kuráž 2009).

Obrázek č.5. (Kuráž 2009)

Kde d představuje šířku koryta a vyjadřuje odtok do koryta a h vyjadřuje odtok pod korytem. K je Darcyho hydraulická vodivost a úhel α je sklon povodí. Dále můžeme dle kapitoly 3.1 charakterizovat Boussmo jako deterministický, celistvý,

1D

model

s koncepčními

prvky,

který

počítá

hydrologického cyklu. Vhodný je pro povodí o lokálním měřítku.

- 22 -

určité

komponenty

3.2.1. Matematický model V této kapitole je popsán konceptuální model Boussma, pro který bude podrobně popsána a odvozena Boussineqova rovnice a popsána rovnice kinematické vlny. 3.2.1.1.

Boussinesqueova rovnice (BR)

Obecné rovnice pro nestacionární trojrozměrné proudění podzemní vody mají velice obtížné řešení, proto se přistupuje k řadě zjednodušení. Velmi často se uplatňuje hydraulický přístup a zavedení některých dalších předpokladů. Hydraulický přístup je založen na způsobech řešení, která předpokládají, že většina zvodní má malou výšku ve srovnání s horizontálními rozměry, to vede k předpokladu, že proudění má převážně vodorovný charakter a jeho vertikální složky se zanedbávají. Při uvažování horizontálního proudění se ekvipotenciály berou jako vertikální přímky. Převaha horizontálního proudění ve zvodní je základem Dupuitových postulátů (Valentová 2007). Dupuit své řešení proudění ve zvodni s vlnou hladinou publikoval v roce 1863. Řešení je založené na zjednodušujících postulátech., které vycházejí z předpokladu, že sklony hladiny podzemní vody jsou většinou velmi malé 1/1000 až 1/10000 a proto je možné považovat horizontální směr proudění (Valentová 2007). Dupuitovy postuláty mají následující znění: • hydraulická výška H (x,y,z) je rovna výšce hladiny podzemní vody h (x,y), proudnice jsou vodorovné přímky a ekvipotenciály svislice •

gradient potenciálu je dán sklonem volné hladiny a je po svislici konstantní

- 23 -

Obrázek č.6. Dupuitovy postuláty (Valentová 2007)

Dupuit vyšel z předpokladu, že pokud je úhel θ velmi malý, přichází v úvahu nahrazení úhlu sin θ = dh/ds sklonem hladiny tg θ = dh/dx . Ekvipotenciály jsou brány jako svislice a hydraulická výška není funkcí vertikální souřadnice z (tzn. H = h(x) místo H = h(x,z)). Hustota toku se pak pomocí Dupuitových axiomů dá psát ve tvaru (Valentová 2007):

vx = − K

dh , h = h(x) dx

(1)

Pro řešení proudění podzemní vody na nakloněném nepropustném podloží se převážně používá Boussinesqových aproximací, které byly odvozeny pro řešení drenážní soustavy na svahu. Tyto aproximace (obr. 7,8) vycházejí z dvou různých verzí Dupuitova postulátu aplikovaného na nakloněné nepropustné podloží: 1) Ve své první publikaci v roce 1877 vycházel Boussinesq z předpokladu, že hladina podzemní vody a proudnice jsou skoro rovnoběžné s nakloněným nepropustným podložím a proto je hydraulický potenciál konstantní v rovině kolmé na nepropustné podloží. 2) V druhé publikaci v roce 1904 uvedl Boussinesq teorii, že proudnice jsou horizontální, což je základní Dupuitův předpoklad. Tento postup je jednodušší a je určen pro mírnější svahy (Pavlásek 2005).

- 24 -

z M ϕ(x) srov. rovina

h x . sinθ

h.cos θ x(M) x

θ x .Obrázek č.7. Schéma Boussinesqueovy první aproximace (BPA) (Pavlásek 2005). z N h

ϕ(x)

srov. rovina

x . tanθ x (N)

θ x

Obrázek č.8. Schéma Boussinesqueovy druhé aproximace (BDA) (Pavlásek 2005).

Kde: θ

- sklon nepropustného svahu

φ(x)

- hydraulický potenciál [-]

h

- výška hladiny

x,z

- označení osy koordinačního systému

V praxi se dle Valentové (2007) zpravidla pro řešení nestacionárních úloh proudění podzemní vody používá Boussinesqueova rovnice ve tvaru:

∂  ∂h  ∂  ∂h  N S ∂h h + h + = ∂x  ∂x  ∂y  ∂y  K K ∂t

- 25 -

(2)

Tato nelineární diferenciální rovnice je pro výpočet proudění v homogenním izotropním prostředí, které je dotováno vertikálním přítokem N. V rovnici (2): K …………….Hydraulická vodivost [m/s] S…………….. Storativita [-] h………………Výška hladiny [m] 3.2.1.1.1

BR ustáleného proudění na nakloněné nepropustné rovině

Pro odvození Boussinesqueovy rovnice pro stacionární proudění se uvažuje homogenní prostředí, hydraulická vodivost je proto reprezentována konstantní hodnotou. Nepropustné podloží je nakloněné a hladina podzemní vody je volná. Díky BPA lze pomocí Darcyho zákona psát vzorec pro rychlost proudění podzemní vody jako:

dϕ ( x) dx

(3)

ϕ ( x ) = h cos θ − x sin θ

(4)

v = −K Z obrázku č.7. lze stanovit:

Specifický průtok na jednotku šířky : h

qx = ∫ − K 0

dϕ ( x) dz dx

(5)

Hydraulický potenciál (φ(x)) je konstantní podél osy z a jeho hodnota se mění pouze s osou x, proto lze celý zlomek, kde se vyskytuje, vytknout před integrál: dϕ ( x) qx = − dx

qx = − K

h

∫ K dz

≈

0

dϕ ( x ) h dx

(6)

Dalším krokem je rozepsání φ(x) pomocí rovnice (4) rovnice (6) pro specifický průtok tak získá tvar:

- 26 -

 dh  qx = − Kh  cos θ − sin θ   dx 

(7)

Rovnici (7) vynásobíme 1/cos θ:

qx  dh  = − Kh  − tan θ  cos θ  dx 

(8)

qx dh = − Kh + K h tan θ cos θ dx

(9)

Po vynásobení závorky jsou na pravé straně rovnice dva členy. První je průtok způsobený sklonem hladiny vzhledem k nakloněnému nepropustnému podloží a druhý je průtok způsobený sklonem nepropustného podloží. S narůstajícím sklonem vzrůstá význam druhého členu rovnice (Pavlásek 2005). V této fázi odvození jsou v rovnici (9) dvě neznámé. Specifický průtok (qx) je v podstatě funkcí přítoku, kterou si můžeme vyjádřit z rovnice kontinuity.

z

∆xr

R

q1

∆h q2 ∆x θ

x

Obrázek č.9. Odvození rovnice kontinuity pro BPA za ustáleného proudění (Pavlásek 2005)

Rovnici kontinuity pro konstantní přítok na hladinu podzemní vody R můžeme psát q2-q1 = R ∆xr , kde ∆xr je odvozena pomocí (obr. 9) jako: ∆x.cos θ + ∆h.sin θ

(10)

Rovnici 10 dosadíme zpět do rovnice kontinuity a vydělíme dx obdržíme tak tvar:

dqx dh = R (cos θ + sin θ ) dx dx Rovnici (11) je vydělena cos θ a upravena na tvar:

- 27 -

(11)

1 dqx  dh  = R 1 + tan θ  cos θ dx  dx 

(12)

Spojením rovnic (8) a (12) získáme rovnici:

dh  d    dh  R 1 + tan θ − Kh  − tan θ   =  dx  dx    dx  R + R tan θ

(13)

dh d  dh  dh = − K  h  + K tan θ dx dx  dx  dx

(14)

Vydělíme celou rovnici (14) K:

R R dh d  dh  dh + tan θ = −  h  + tan θ K K dx dx  dx  dx

(15)

Dostáváme výsledný tvar diferenciální rovnice pro ustálené proudění pro BPA:

d  dh  dh  R  R  h  − tan θ 1 −  + = 0 dx  dx  dx  K  K 3.2.2.1.2

(16)

BR neustáleného proudění na nakloněné nepropustné rovině

V tomto případě se mění výška hladiny v závislosti na změně vertikálního přítoku. Rovnici kontinuity pro měnící se přítok na hladinu podzemní vody R můžeme vyjádřit:

∆q = q2-q1 = R ∆xr - ∆h/∆t . ∆x. µ

(17)

∆xr je odvozena pomocí. (Obr. 10) jako: ∆xr = dx.cos θ + dh.sin θ

∆xr

(18)

R

∆h/∆t

q1

q2

∆x θ

x

Obrázek č.10. Odvození rovnice kontinuity pro BPA za neustáleného proudění (Pavlásek 2005).

- 28 -

Člen µ představuje aktivní efektivní pórovitost, díky které rovnice vychází v objemových jednotkách. Kombinací rovnic (17) a (18) a vydělením výsledné rovnice cos θ.dx je získána rovnice kontinuity pro neustálené proudění: 1 cos θ

dh  µ dh  dq     = R 1 + tan θ − dx  cos θ dt  dx  

(19)

Spojením rovnic (19) a (9) obdržíme Boussinesqovu rovnici pro neustálené proudění v homogenním prostředí, její tvar je následující: d  dh  dh  RR µ dh  h  − tan θ 1 −  = dx  dx  dx  K  K K cos θ dt

(20)

Po drobných úpravách obdržíme tvar rovnice, který je použit pro výpočet průtoku podzemní vody v modelu Boussmo:

µ

dh R  dh  d  dh  dh = 1 + tan θ  +  h  − tan θ K cos θ dt K  dx  dx  dx  dx

(21)

3.2.1.2 Rovnice kinematické vlny Jelikož se na povodí Modravy 2 povrchový odtok vyskytuje pouze při výjimečných událostech jakou popsal Pavlásek (2008), nebude tato rovnice při simulacích modelu použita. Rovnice vychází ze schématu (obr. 11)

Obrázek č.11. Schéma pro odvození kinematické vlny (Kineros 2).

- 29 -

Z pohledu velmi malého měřítka představuje povrchový odtok extrémně komplexní 3D proces. Avšak ve větším měřítku na něj může být nahlíženo jako na 1D odtokový proces. Tento proces je vyjádřen vztahem:

Q = α ⋅ hm

(22)

kde Q je průtok na jednotku šířky[L3.T-1], α ,m jsou parametry vztažené ke sklonu, drsnosti a odtokovému režimu. Rovnice (22) se spojí s rovnicí kontinuity: ∂h ∂Q + = q ( x, t ) ∂t ∂x

(23)

kde t je čas, x je vzdálenost podél svahu, q je hodnota bočního přítoku [L2.T-1]. Pro povrchový odtok spojením rovnic (22) a (23) získáme rovnici kinematické vlny ve tvaru: ∂h ∂h + α mh m −1 = q ( x, t ) ∂t ∂x

(24)

Rovnice kinematické vlny jsou zjednodušením Saint Venantových rovnic (Kineros2). Tato diferenciální rovnice je v modelu řešena metodou konečných diferencí (Kuráž 2009).

3.2.3

Procedurální model Algoritmus je založena na přístupu schematizace pomocí nádrží, kde odtok

z povodí je popsán Boussinesqueovou rovnicí (BR) a povrchový odtok rovnicí kinematické vlny (KV). Schéma modelu (obr. 12)

- 30 -

srážka přetečení nádrže

max.přítok =K (zbytek jde do nádrže pro povrchový odtok)

přetečení nádrže

Nádrž podpovrchového odtoku (konečný objem)

přetečení nádrže se dostává do nádrže pro povrchový odtok

BR

KV

Nádrž povrchového odtoku (nekonečný objem)

obě vyústění plní tok v místě uzavěrového profilu

Obrázek č.12 Schéma systému nádrží v modelu Boussmo (Kuráž 2009)

Procedurální model je napsán v programovacím jazyce F (podskupina standardu Fortran 95-2008). Inicializační procedura načte data ze vstupní složky uložené in/boussmo.conf. Konfigurační soubor je umístěn v příloze č.1, kde uživatel definuje parametry, kterých je celkem 16. Počáteční podmínka je vyjádřena jako jakási průměrná srážka z období před srážkovou událostí. Použitím této hodnoty je počáteční hladina v podpovrchové nádrži počítána za použití stacionární verze Boussinequeovi rovnice (dále BR). Počáteční plnění nádrže je počítáno integrací počáteční vodní hladiny. Poté jsou srážková data načtena a uložena do přiděleného pole. Je spuštěna procedura bouss, která zavádí inicializační test půdních vlastností a srážek. Pokud je objem srážka menší než Darcyho nasycená hydraulická vodivost, pak je celý objem této srážky v nádrži pro podpovrchový odtok. Pokud je objem srážky vyšší, hodnota objemu hydraulické vodivosti je ponechána v nádrži podpovrchového odtoku a zbytek připadá pro nádrž, která je na řešena rovnicí kinematické vlny. Procedura volume je volána, zjišťuje objem nádrže, pokud je předchozí dopadnuvší objem srážky pro odhadovaný časový krok vyšší než zbytek objemu podpovrchové nádrže, pak je přebytek zprůměrován podle časového kroku pro výpočet srážky pro kinematickou vlnu. Dále je BR řešena iterativně. Řešič kinematické vlny je volán na konci této procedury se srážkovými daty vyčíslenými jak bylo popsáno v předchozím odstavci. Kvůli vyšší nelinearitě rovnice (24) v porovnání s nelinearitou v BR (21) je potřebný časový krok pro konvergenci - 31 -

menší než potřebný čas pro BR při nevhodně zvoleném časovém kroku má rovnice oscilační chování. Bouss procedura volá proceduru kinematix, která určuje pouze časový krok pro řešič kinematické vlny. Tato procedura volá dále privátní proceduru solver, která zkouší řešit iterativně rovnici kinematické vlny. V každé iteraci je ověřována konvergence – chyba má být menší než v předchozí iteraci, pokud ne procedura je ukončena s definicí kódové chyby a procedura kinematix zkouší volat řešič pro případ sníženého časového kroku. Procedura kinematix v momentě, když kumulativní čas řešené kinematické vlny dosáhne časové periody definované jako časový krok BR rovnice (Kuráž 2009). 3.3

Předchozí srážkový index. Srážko–odtokový

konceptuální

model

BOUSSMO,

řešící

numericky

boussinesquovu rovnici pro podpovrchový odtok pomocí metody konečných diferencí a rovnici kinematické vlny pro odtok povrchový, zanedbává evapotranspiraci a stav nenasycené zóny a je tak vhodným modelem pro tropické oblasti zejména při období dešťů, kdy se dá předpokládat úplné nasycení celého půdního profilu. Proto při aplikaci tohoto modelu na experimentální povodí Modravy 2, je nezbytné pro použití v těchto podmínkách učinit odhad vlhkosti půdy, respektive odhad nasycenosti povodí jakožto podstatnou složku počátečních podmínek vstupujících do modelu. Model Boussmo můžeme zařadit dle Daňhelky (2003) do takzvaných modelů výzkumných, pro které je charakteristické přesnější popis problému, jejichž hlavním cílem je studium problému S-O vztahů. Tyto modely je schopen provozovat pouze úzký okruh uživatelů, často zainteresovaných na vývoji modelu. Modely jsou často aplikovány na experimentální povodí s nadstandardní pozorovací sítí velkého počtu charakteristik pro srážkoodtokový proces Počáteční stav nasycenosti půdy ovlivňuje hodnoty potenciální retence, která představuje největší možnou retenci daného povodí. Retenční kapacita půdy je dále ovlivněna: tloušťkou půdní vrstvy, průměrnou pórovitostí půdy (Diermanse 2001). Retence je například jednou ze složek pro odvození čísla CN – křivek. Metoda CN (Curve Number Method) byla vyvinuta v USA Službou pro ochranu půd (US–Soil Conservation Service – US SCS). Metoda umožňuje odvození objemu „přímého odtoku“ a kulminačního průtoku na zemědělsky a lesnicky využívaných povodích i na povodích urbanizovaných do velikosti plochy povodí cca. 5 km2 (Hrádek, Kuřík 2002).

- 32 -

H o Ra = Hd Rp

Metoda vychází ze vztahu:

Ho …výška přímého odtoku [mm] Hd …výška výpočtového deště [mm] Ra …aktuální retence povodí [mm] Rp …potenciální retence povodí [mm] a je charakterizována třemi skupinami předchozích vláhových podmínek (PVP) podle úhrnu předchozích dešťů za 5 dnů. Metoda PVP (předchozí vláhové podmínky) se v modelech pro určení obsahu vody v půdě používá zřídka. Přehled skupin PVP (tab. 1).

Skupina PVP

Celkový úhrn předchozích srážek v [mm] za 5 dnů v období mimovegetačním

vegetačním

I

<13

<36

II

13-18

36-53

III

>28

>53

Tabulka č. 1. Přehled skupin předchozích vláhových podmínek.

Pro model BOUSSMO byla pro odhad nasycenosti povodí vybrána empiricky založená metoda předchozího srážkového indexu API (z ang. zkratky Antecedent Precipitation Index) označovaného často i jako úhrn předchozích srážek UPS. Proceduru API poprvé definoval A.M. Kohler (1951). Obecný tvar vypadá takto:

APIn =

n

∑C n =1

i

. Pi

[mm]

(25)

kde: n – celkový počet dní před výskytem příčinných srážek, obvykle se n volí 5 nebo 30 i – pořadí dne počítané nazpět ode dne, ke kterému je API určován C – evapotranspirační konstanta, pro naše podmínky obvykle C = 0.93 P – denní úhrn srážek v milimetrech v i-tém dni před výskytem příčinných srážek

- 33 -

Za konceptuální model efektivního vodního vstupu je považován vztah: Weff = W-ztráty, kde W je celkový srážkový vstup během události a ztráty = ET + ∆Sc + ∆D + ∆Θ , kde ET je část vody evapotranspirovaná během události, ∆Sc je hodnota přiřazená hodnota zásobě na vegetačním pokryvu, ∆D je hodnota určená depresní zásobě tj. voda přidaná do jezer, rybníků, kaluži a podobně, ∆Θ je hodnota určená pro zásobu půdní vody během události. Jelikož srážkové události bývají obvykle krátkého trvání a jsou doprovázeny vysokou vlhkostí a malým slunečním zářením, proto je ET obvykle malé. Kapacita zásobnosti vegetačního krytu je v řádu 1 mm × index listové plochy a takto definován se poměrně rychle zaplní. ∆Sc je taktéž obvykle zanedbatelné pro srážky, které generují výraznou odpověď povodí v podobě přímého odtoku v literatuře též často označovaného jako bouřkový nebo rychlý odtok. Depresní zásoba je prostorově variabilní a tudíž obtížně odhadnutelná, proto se obvykle kombinuje v konceptuálně založených modelech se zásobou půdní vody. Tyto skombinované zásobní komponenty jsou typicky modelované stejně jako proces infiltrace. Takto je koeficient Weff / W určen velkou měrou stupněm kapacity možné zásoby, která již obsahuje nějakou půdní vlhkost. Tedy kolik vody se ještě může infiltrovat. Jednou z operativních metod k nalezení vztahu mezi efektivní srážkou a předchozími podmínkami na povodí je právě API (Dingman 2002). Po předchozí atmosférické srážce se vyjmenované zásobní komponenty uvedené ve vztahu ztráty = ET + ∆Sc + ∆D + ∆Θ postupně prázdní. Tento proces je aproximován přes empirický srážkový index Ia(d), který je počítán každý den na základě: I a (d ) = I a (0).k d

(26)

kde Ia(0) je hodnota pro den se srážkou, k je konstanta (obvykle 0.80< k <0.93), a d je počet dní od poslední dešťové srážky. Hodnoty Ia(0) a k jsou empiricky určené pro jednotlivá povodí. Významově Ia(0) reprezentuje celkovou zásobu povodí při povrchu (obvykla vyjádřená jako hloubka vody), a Ia(d)je množství vody z předchozí srážky, která zůstává v zásobě do dne d. Poté je nutné najít empirický vztah mezi Weff a Ia(d) pro minulé srážky. Tento empirický vztah byl do detailu probrán v článku (Linsley, Kohler 1951). Ve 40.letech 20. století byla snaha většiny hydrologů upřena na zjednodušení vztahu mezi odtokem a srážkami. Fundamentální problém představovalo odhadnutí

- 34 -

hodnoty odtoku z daného množství srážkového úhrnu. V roce 1951 Kohler a Linsley popsali vztah mezi srážkami a odtokem pomocí grafické metody koaxiálních vztahů, která vzala v potaz sezónní vlivy, předchozí podmínky na povodí, trvání srážky, srážkový úhrn a její pomocí vymezili část srážky, která se podílí na odtoku. Před nimi byla snaha o uplatnění infiltrační teorie, která však narážela na velkou variabilitu přírodních povodí a řešení bylo prakticky nemožné i s hustou sítí srážkoměrných stanic. K tomu navíc se přímá aplikace infiltrační teorie dala použít pouze pro stanovení části povodňového hydrogramu a to povrchového odtoku. Předpovědi na řeku, ale vyžadují celkový tok, tedy včetně hypodermického a podzemního odtoku. Právě tyto dvě komponenty představují hlavní složku v povodňovém hydrogramu pro některá povodí. Výběr správných parametrů byl určujícím faktore při navrhování techniky na odhad odtoku. Po intercepci, infiltraci a naplnění depresní zásoby nastává odtok. V tomto se logicky jevilo nějaké použití rozdílu mezi srážkou a odtokem jako závislé proměnné. Tato diference se nazývala jako ztráta nebo také zdržení respektive retence povodím. Po zjištění tohoto zdržení a srážky může být odtok spočítám přímým odečtením. Velikost doby zdržení od dané srážky závisí na vlhkostním nasycení půdy na začátku srážky a na charakteristikách samotné srážky jako jsou množství, intenzita atd. Zatímco charakteristiky příslušného deště mohou být stanoveny pomocí adekvátní srážkoměrné sítě, přímé stanovení vlhkostních podmínek pro celé povodí je značně složité. Půdní typy, povrchové charakteristiky, rozdílný vegetační pokryv a využití půdy to vše přispívá ke komplexitě problému stanovení půdní vlhkosti. Číselně měřitelné faktory, které byli využity pro vyjádření vlhkostních podmínek byly dny od posledního deště, odtok na začátku srážky a předchozí srážka. První je necitlivé, druhé je dobrý ukazatel ve vlhkých regionech, ovlivněn sezónností a nepostihuje změny způsobené deštěm během týdne. Předchozí srážka je univerzálně přijatelná a poskytuje dobré výsledky pokud je dobře odvozen koeficient, navíc zahrnuje vliv sezónnosti a teploty. Obecně dle Kohler a Linsley (1951) má předchozí srážkový index tvar : I=b1.P1 + b2.P2 + b3.P3 + …+ bi.Pi

(27)

…Pi představuje množství srážek, které spadnou v i dni před srážkou …bi je konstanta, která je funkcí času, pokud je vyžadována hodnota indexu ze dne na den většinou bi s časem logaritmicky klesá. Jinými slovy, během doby bez srážky : It = I0 kt kde t je počet dní mezi It a počátečním indexem I0, pro t =1

- 35 -

(28)

I1 = k I0

(29)

Takto je index určitého dne roven indexu předchozího dne vynásobeného koeficientem k. Jestliže déšť nastane v kterýkoliv den, množství naměřené srážky se přidá k indexu. Názorněji (obr.13)

Obrázek.č 13. Průběh API po dvacet dní.

K výpočtu předchozího srážkového indexu postačují údaje o denních srážkových úhrnech z daného povodí pro sledované období. Teoreticky hodnota koeficientu zdržení k by mohla představovat fyziografické charakteristiky povodí, ale ze zkušeností s tímto faktorem se ukázalo, že není rozhodující a jeho rozmezí se např. pro východní a centrální část Spojených států amerických pohybuje mezi 0.85 do 0.9. Pro naše zeměpisné šířky se používá hodnota k = 0.93 (Kovář 1990). Tato hodnota k byla také na příklad použita v roce 2002 pro výpočet UPS na povodí Jizera pro hydrometeorologické vyhodnocení katastrofální povodně v srpnu 2002 (Řičicová 2003). Důležité je zvolit vhodnou hodnotu počátečního srážkového indexu API (0). Využívá se možnosti zahájit výpočet indexu na počátku suchého období s nízkou hodnotou indexu nebo zahájit výpočet 2-3 týdny před první srážkou s předpokládanou hodnotou indexu rovnou normální 10 denní srážce za roční období, která aproximuje průměr hodnoty indexu pro zájmové území (Linsley, Kohler 1951). API (0) [mm] reprezentuje průměrný stav nasycenosti povodí před samotným výpočtem předchozího srážkového indexu. API bylo původně koncipované pro reprezentování aktuálního stavu půdní vlhkosti v modelech předpovídající srážko-odtokovou událost. Základ metody API vyšel z potřeby nalezení vztahu mezi snadno měřitelnou veličinou, zde tuto veličinu

- 36 -

představují srážky a obtížněji stanovitelnou veličinou, půdní vlhkostí. Zatímco srážkové charakteristiky mohou být určeny z adekvátní sítě srážkoměrných stanic tak přímé stanovení vlhkosti celého povodí je značně obtížné kvůli velké heterogenitě zájmového území. Půdy inklinují k tomu, že jsou lokálně různorodé v jejich charakteristikách, tak že infiltrační kapacita a rychlost generování povrchového odtoku se může měnit od místa k místu. Na mnoha místech, zejména porostlých vegetacích, srážky velmi zřídka přesáhnou infiltrační kapacitu půdy dokud se půda zcela nenasytí (Beven 2002). Lepší odhady půdní vlhkosti mohou přijít s nárůstem dostupnosti prostorových dat z dálkového průzkumu země jako jsou srážky a satelitní snímky o různých vlnových délkách zahrnující aktivní a pasivní mikrovlnné senzory použité pro odhad povrchové půdní vlhkosti. Vztah mezi srážkami a půdní vlhkostí spočívá ve faktu, že půdní vlhkost klesá logaritmicky nebo asymptoticky s časem, ve kterém nedošlo ke srážkové události (Linsley, Kohler 1951). API je možné popsat jako váženou sumaci denních srážek. Váha je dána každému srážkovému dni obvykle exponenciální nebo reciproční funkcí času. Nejstarší srážková událost obdrží nejmenší váhu. API se uplatňuje v mnoha srážkoodtokových modelech jako komponenta pro odhad půdní vlhkosti. Četné experimenty a početné studie provedené po celém světě zabývající se problematikou formování svahového odtoku jasně ukazují, že počáteční obsah půdní vody má přímí vliv na infiltrační kapacitu a následkem toho na povrchový odtok (Descroix, Nouvelot,Vauclin 2002). API se používá jako klíčová proměnná pro povrchový odtok v prostředích s malou datovou dostupností (Fedora, Beschta 1988) např. v příbřežních oblastech státu Oregon. V subtropických horách severního Mexika, kde se půda i snahy chovají dle Hortonových předpokladů, jak je ostatně často pozorováno i v dalších tropických nebo subtropických oblastech, byl vyvíjen a popsán Descroixem, Nouvelotem a Vauclin (2002) jednoduchý deterministický model NAZASM, který je založen právě na API. Tito autoři ve svém článku upozorňují na nutnost, aby data určená ke kalibraci obsahovala srážkové události jak ze suchých tak i na srážky bohatších let. API není vhodným podkladem pro nevýrazné odtokové vlny a také zimní a jarní epizody spojené s táním sněhu, protože v těchto případech přestává reprezentovat aktuální vláhové poměry v povodí. Proto je výběr srážkoodtokových událostí omezen

- 37 -

především jen na vegetační období. Tento fakt může v sušších a rovinných povodích znamenat výrazné zmenšení podkladového datového souboru (Daňhelka 2003). Model, využívající předchozí srážkový index je ku příkladu APIc (Antecedent Precipitation Index Continuous). Model je vhodný pro kontinuální operativní předpovědi. Koncept modelu je vystaven na koaxiální grafické korelaci odvozené ze souboru srážkových epizod a příslušných odtokových vln. Výhoda této metody spočívá v tom, že je logicky koncipována s využitím empirických a intuitivních vztahů, které mohou být velmi názorně prezentovány ve formě grafu. APIc byl používán v rámci hydrologického předpovědního systému AquaLog na povodí Jizery a Sázavy (Daňhelka 2003). Aqualog je víceúčelový vodohospodářský modelovací systém určený pro podporu rozhodovacích procesů v oblasti vodního hospodářství. Aqualog byl vyvinut skupinou odborníků zaměřených na hydrologii, hydrauliku, kvalitu vody a výpočetní techniku. Jedná se o integrovaný programový prostředek určený pro simulace, předpovědi a řízení odtokového procesu a kvality vody v historickém i v reálném čase. Kromě oblasti operativní hydrologie může tedy pracovat i jako prostředek pro odvození návrhových charakteristik vodohospodářských opatření a při sledování různých scénářů kritických situací (Aqualog). V roce 2001 byl však nahrazen APIc modelem Sacramento (SAC-SMA), zároveň však zůstala zachována možnost provádění výpočtu modelu APIc, který je dále využíván pro srovnávací studie. Lze jej charakterizovat jako metodiku vycházející z tradičních koaxiálních funkcí a používající Sittner-SchaussMonro algoritmus jako základ. Používání modelu má v našich podmínkách tradici a s manuálním postupem tvorby předpovědi je jeho používání podepřeno zkušenosti z hydrologických pracovišť jako jsou již zmiňované povodí Sázavy a Jizery, ale také povodí Vltavské kaskády po vodní dílo Orlík (Buchtele 1972). Snadná konstrukce a řešení modelu však sebou nesou i nedostatky. V modelu se uvažuje pouze odtok povrchový a podzemní. Další složky odtoku se zanedbávají, jsou to např. podpovrchový, přímý z nepropustných ploch, podzemní s krátkou a dlouhou dobou zdržení. Problémy může činit skutečnost, že koaxiální vztah je odvozen z celých odtokových epizod přitom pro samotné modelování je nejdůležitější nástup povodně a odezva dešťových srážek. Samozřejmě značnou nejistotu do modelu vnáší vykreslování čar v okrajových polohách, kde se v reálných aplikacích dle Daňhelky (2003) uživatel pohybuje nejčastěji. Nezávislé proměnné v modelu APIc jsou uvedeny (tab. 2):

- 38 -

Proměnná

Jednotka

Popis

P API TS

[mm] [mm] [hod]

úhrn příčinné srážky index předchozích srážek trvání příčinné srážky

WK

[-]

kalendářní číslo týdne (slouží jako charakteristika ročního evapotranspiračního cyklu)

RI

[-]

retenční index

Tabulka č.2. Nezávislé proměnné modelu APIc

Předchozí stav vlhkosti v povodí je vyjádřen indexem předchozích srážek pro 5 a pro 30 dní. Pro používání modelu APIc vidí Daňhelka (2003) budoucnost v možnosti testování radarových vstupů do modelu na vybraných povodích, kde již byl v minulosti využíván. Tam kde by ho bylo třeba nově zavádět a kalibrovat se zdá být neproduktivní. Na mnoha povodích je tento model zpravidla nahrazován modelem Sacramento. Ukázka koaxiálně grafického korelačního vztahu srážka – odtok v podobě grafu (obr. 14)

Obrázek.č.14. Koaxiálně grafický korelační vztah srážka-odtok.

Pro potřeby modelů, které v sobě obsahují komponentu API pro odhad stupně nasycenosti povodí je vhodné stanovit maximální hodnotu API (APImax). Při dostatečně vydatné srážce, která svou intenzitou překročí infiltrační rychlost půdy může

- 39 -

samozřejmě dojít k úplnému nasycení nenasycené zóny a k tvorbě povrchového odtoku. Dosažení APImax se předpokládá pouze u velmi významných povodňových událostí. Například model CLS používaný k modelování celkového odtoku, při překročení prahové hodnoty APImax přechází na vyšší transformační funkci s větším odtokovým koeficientem a strmějším nástupem (Blažková 1993). Anderson (1984) uvádí hodnotu APImax v rozsahu mezi 20-25 cm. 3.3.1 Porovnání modelů APIc a SAC-SMA Jeden z klíčových částí AHPS – Advanced Hydrologic Prediction Service navrhnutý NWS – National Weather Service je Sacramento Soil Moisture Accounting (SAC-SMA) model. Sacramento nahrazuje epizodní model předchozího srážkového indexu (APIc). Nakalibrované SAC-SMA pomůže zlepšit přesnost předpovědí RFC (River Forecast Centers), protože aplikace modelu APIc byla založena na většině územích na regionálních empirických vztazích, než na specifických charakteristikách povodí. Takže některá RFC nakalibrovala parametry SAC-SMA pro specifické charakteristiky daných povodí. I u nás došlo v roce 2001 k nahrazení modelu APIc za SAC-SMA v modelovacím systému Aqualog, který je používán ČHMÚ na celém povodí Labe. V práci M.B Smitha (2000), která hodnotila přínosy v zavedení Sacramenta je porovnání výsledků modelu Sacramento s modelem APIc, jelikož oba modely jsou uživatelsky velmi modifikovatelné, byla při porovnání snaha vyloučit tento lidský faktor tím, že simulace predikce se prováděly na mnoho let dopředu. Při porovnání byla použita kritéria RMSE (root mean square error) a chyba kulminace mezi pozorovaným a kulminačním průtokem. Výsledky v příloze č. 2 převzaté z práce Smitha (2000). 3.4

Kalibrace a optimalizace parametrů

3.4.1 Kalibrace a validace Před použitím hydrologického modelu na vybrané povodí je nezbytné učinit kalibraci a validaci. Kalibrace je proces hledání optimálních hodnot vstupních parametrů a koeficientů, které model obsahuje. V podstatě se jedná o iterační postup, kdy při známých cílových hodnotách jsou měřeny hodnoty vstupních parametrů tak, aby se výstupy s určitou přesností, kterou určuje zvolená objektivní funkce nebo soubor objektivních funkcí, blížily naměřeným hodnotám. Kalibrace slouží k odhadnutí hodnot parametrů, které nemohou být přímo změřeny z terénu či mapových podkladů ani nijak

- 40 -

v dostatečné míře odvozeny. Kalibrace je tak nutná pro většinu celistvě koncepčních a empirických modelů, které obsahují parametry, jejichž hodnoty musí být odhadnuty pomocí tohoto procesu. Parametry není možné pro jeden model stanovit obecně, pro každé povodí je třeba kalibrovat zvlášť. Kalibrace by se měla opakovat i při jakémkoliv použití modelu, kdy jsou použita data, která jsou naměřena za jiných podmínek na povodí, než jaké byly na povodí při měření dat, ze kterých byla kalibrace provedena. Data pro kalibraci by měla odrážet pokud možno všechny jevy, cykly, fenomény na povodí. Toho se často dosahuje tím, že se použije co možná největší soubor dat, která jsou vhodná. I zde však platí, že kvalita informací obsažená v datech je podstatnější než délka časové řady (Yapo et kol. 1996). Součástí kalibrace je citlivostí analýza. Ta slouží pro vymezení vlivu změny parametru na výsledek modelu. Každou kalibraci je nutno ověřit validací. Smysl validace spočívá v ověření, že pro nakalibrované hodnoty daných parametrů jsou výstupy blízké měřeným veličinám. Pro validaci se používá taková část z datového souboru, která nebyla použita při kalibraci. Jak uvádí Daňhelka (2003) je třeba si uvědomit, že získání schopnosti dobře kalibrovat koncepční model není krátkodobá záležitost. Zatímco matematické vztahy koncepčního modelu jsou poměrně jednoduché, kombinace jednotlivých procesů, komponent, chyb modelu a dat vytváří z procesu kalibrace komplexní systém, který není deterministický a kdy pracujeme v oblasti numerické neurčitosti. Významnou roli zde hraje profesionální zkušenost uživatele, který model kalibruje a jeho intuice případně zvláštní nadání. Beven (2002) vymezuje tři druhy kalibrace: •

Metody kalibrace, které předpokládají existenci optimálního souboru parametrů, ale které ignorují odhad neurčitosti simulovaných dat. Tyto metody mají rozsah od jednoduchých jako je “pokus a omyl” až po sofistikované metody automatické optimalizace.

•

Metody předpokládající existenci optimálních parametrů, ale dělající jistý odhad objektivní funkce kolem tohoto optima k odhadnutí intervalu neurčitosti výsledku. Tyto metody se označují jako analýzy spolehlivosti.

- 41 -

•

Metody neuvažující existenci optimální sady parametrů a jsou založeny na principu ekvifinality. Ekvifinality předpokládá existenci více přijatelných modelů, které jsou ve shodě s naměřenými daty.

Důvody obtížnosti kalibrace a optimalizace parametrů (Beven 2002, Daňhelka 2003): •

rozdílnost v měřítku dostupných měřících technik a měřítku, ve kterém je hodnota parametru požadována. Model počítá s většími prostorovými úseky, než jaké jsme schopni změřit.

•

při optimalizaci se vypočtené výsledky porovnávají s naměřenými, která ovšem mohou obsahovat chyby.

•

náhodné nebo systematické chyby při měření vstupních dat (srážky, teplota) používané ke stanovení vstupních podmínek v povodí (čas, prostor)

•

chyby v důsledku nekalibrovaných parametrů

•

chyby v důsledku nekompletní struktury modelu Úspěšné použití matematických modelů předpokládá nejen vyhodnocení shody

mezi měřenými a modelově vypočtenými výstupními údaji, ale zejména stanovení hodnot parametrů modelu. Změny hodnot ovlivňují citlivost modelu. Je třeba vymezit rozsah změn tak, aby parametrům zůstal fyzikální význam a ještě poskytovaly dobré výsledky (Kovář 1990). Nejjednodušší analýza se dá provést metodou pokusu a omylu, ale ta se při větším počtu parametrů stává neúnosnou. Jinou alternativou je výběr všech sad parametrů z celého parametrického prostoru tj. z intervalu všech možných sad parametrů a vyhodnotit výsledky kriteriem míry dobré shody. Grafem těchto výsledků je průběh daného kritéria (objektivní funkce) v prostoru parametrů. (obr. 15) je znázorněn průběh objektivní funkce v podobě hory, kde vrchol je tvořen hodnotami optimálních parametrů. Toto optimum často není lehké najít (Beven 2002).

- 42 -

Obrázek č.15. Průběh objektivní funkce v parametrickém prostoru (Beven 2002)

Důvody slabé identifikovatelnosti parametrů (Beven 2002, Yapo et kol. 1996): 1.

Interakce parametrů: vede k více lokálním optimům (obr.16) nebo tzv. hřebenům (obr.17) v průběhu objektivní funkce. Různé dva parametry dávají stejnou nebo podobnou hodnotu míry shody simulovaných a měřených dat. Tyto interakce mají jasné fyzikální opodstatnění: jestliže

se model skládá

z hortonovského odtoku, podpovrchového odtoku a saturovaného odtoku, pak pravděpodobně existují soubory parametrů s dobrou shodou, které využívají vždy především jeden z těchto odtoků. Lokální optima lze pak nalézt v různých částech parametrického prostoru.

Obrázek č.16 (Beven 2002)

2.

Obrázek č. 17 (Beven 2002)

Nestacionarita parametru: Parametr je ovlivněn řadou různých charakteristik povodí nebo daty, která nejsou použita během kalibrace.

3.

Neinformovanost dat: Data dostatečně nepopisují hydrologické podmínky tak přesně aby se dali identifikovat parametrem.

- 43 -

4.

Neadekvátnost kritéria: Objektivní funkce nedokáže získat informaci obsaženou v datech

5.

Necitlivost: Změny parametru mají nesignifikantní vliv na výstup modelu (obr. 18) parametr je “plochý“. V takovém případě kdy je parametr nedůležitý, může být fixován číselnou hodnotou (Yapo et kol 1996).

Obrázek č. 18. (Beven 2002)

Přestože je použití numerických kriterií významné, hodnota grafického porovnání simulovaných a měřených hydrogramů by neměla být přehlížena. Subjektivní analýza grafu poskytuje názorný celkový dojem na spolehlivosti modelu (Daňhelka 2003). 3.4.2 Automatická optimalizace parametrů Jelikož je proces kalibrace obtížný a komplexní je potřeba robustní a spolehlivé automatické kalibrace. Dle Yapo et kol. (1996) zavedení automatické kalibrace vyžaduje správný výběr kalibračních komponent jako je: 1. Soubor kalibračních dat. 2. Objektivní funkce. 3. Optimalizační algoritmus. 4. Rozmezí parametrů, které bude v prostoru parametrů hledáno. 5. Vhodnou validační proceduru. V současnosti stále probíhá diskuze, zda automatická kalibrace by mohla zcela nahradit tradiční časově náročné a subjektivní manuální kalibrační procedury.

- 44 -

Procedury automatických kalibrací, používané v současnosti se většinou testují v rámci jednotlivých modelů. Zobecnit získané výsledky i na ostatní modely je však obtížné. (Daňhelka 2003). V posledních letech se výzkum zaměřuje na vývoj globální optimalizační metody pro automatickou kalibraci nejen hydrologických koncepčních modelů, ale rovněž hydrodynamických modelů. Nalezení globálního extrému u nekonvexních funkcí s lokálními extrémy dosud není uspokojivě vyřešeno. Objevují se nové přístupy obsahující prvky umělé inteligence jako např. genetické algoritmy (GA) či shuffled complex evolution (SCE), které prokázaly být účinný nástrojem především v lokalizaci globálního optima modelu při kalibraci s použitím jednotlivé kriteriální funkce (Daňhelka 2003, Yapo et kol. 1996). Další metody, které naleznou globální extrém funkce jsou např. gradietní metody , které nemají umělou inteligenci či simulované žíhání (Beven 2002). Genetické algoritmy (nepřímá metoda stochastická) : Vychází z analogie evoluce v živé přírodě. GA pracuje se sadou možných řešení označovaných jako “chromozomy“ dále nazývanými jako populace. V základním schématu jsou chromozómy reprezentovány jako binární řetězce. Schéma genetického algoritmu a genetických operátorů popisuje Kučerová et Hrstka (2004) následujícím způsobem: Jako první krok je nezbytné vygenerovat (ve většině případů náhodně) startovní populaci (soubory parametrů) možných řešení, kterým bude přidělena hodnota objektivní funkce. Poté se opakuje sekvenční smyčka tak dlouho dokud nedojde k dosažení určeného kritéria: 1. tvorbě předepsaného počtu nových jedinců (“chromozomů“) použitím genetických operátorů crossing-over a mutace 2. hodnoty optimalizační funkce jsou přiřazeny novým individuím 3. populační velikost je snížena na původní hodnotu operátorem selekce Shuffled complex evolution (nepřímá metoda) : Jedná se o formu genetického algoritmu, která byla vyvinuta pro potřeby srážkoodtokových simulací. Je kombinací gradientní metody a teorií genetických algoritmů. V tomto algoritmu, jsou různá jednosměrná hledání uskutečněna paralelně z každého náhodného startovního bodu. Po každé iteraci vícenásobného hledání jsou hodnoty

- 45 -

aktuálních parametrů smíchány a tvoří simplexy, z kterých se vytváří nový startovací bod pro další iteraci (Beven 2002). Gradientní metoda : Tato metoda byla vyvíjena od šedesátých let kdy vstoupily do modelování počítače. Je založena na postupném zlepšování objektivní funkce. Algoritmus začíná z jednoho bodu objektivní funkce a podle známého gradientu průběhu funkce postupuje směrem kam roste (v případě hledání maxima objektivní funkce) či klesá (v případě hledání minima). Dostupné techniky např. Newtonova metoda nebo LevenbergovaMarquardt mohou být rozděleny do dvou základních typů: První jsou algoritmy vyžadující gradient objektivní funkce definován analyticky pro každý z parametrového prostoru. Tyto metody se v hydrologii moc nepoužívají, protože je často nemožné analyticky definovat změny na pro komplexní strukturu modelu. Mnohem používanější jsou algoritmy přímého hledání, které hledají podle zkušebních pravidel cestu ze součastného bodu s cílem najít zlepšenou hodnotu objektivní funkce. Pokud je metoda použita pro kalibraci parametrů, je nutné začít z více náhodně vybraných bodů v prostoru parametrů. Když jsou si výsledné soubory parametrů podobné, je možné se domnívat, že funkce má jen jedno optimum (Beven 2002). Metoda simulovaného žíhání (nepřímá, stochastická) : Také tato metoda používá náhodně vybrané počáteční body v prostoru parametrů. Název algoritmu vychází z analogie mezi optimalizovanými parametry modelu a chováním částic v chladnoucí tekutině. Na počátku jsou částice náhodně rozděleny v prostoru. Jak tekutina ztrácí teplotu, minimalizuje se energie systému. Pokud je proces příliš rychlý minimální energie se projeví lokálně. Pokud tekutina chladne pomalu, může být výsledkem globální energetické minimum. Společným pravidlem všech variant této metody je přijetí nového souboru parametrů. Jestliže má tento soubor lepší hodnotu objektivní funkce než soubor minulý, je přijat. V případě, že lepší není může být ještě přijat, a to s určitou pravděpodobností, která se s klesající teplotou, tedy v průběhu iterací, snižuje. Přijímání souborů parametrů s horší hodnotou objektivní funkce zabezpečuje, že algoritmus nebude zastaven na lokálním optimu (Ředinová 2004).

- 46 -

3.4.3 Objektivní funkce Cílem optimalizace je nalézt globální nebo lokální extrém vybrané objektivní funkce. Hledaným extrémem je buď minimum nebo maximum funkce. Automatická optimalizace parametrů spočívá především v jejich funkční specifikaci. Všechny parametry musí být popsány funkčními vztahy společně s proměnnými. Soustava těchto vztahů je potom optimalizována podle vybrané objektivní (kriteriální) funkce. Výběr objektivní funkce je dán účelem modelu (Kovář 1990). Mezi nejpoužívanější patří objektivní funkce založené na sumě chyb nejmenších čtverců nebo chyb variance. Objektivní funkce můžeme rozdělit na statistické a hydrologické (YU et kol. 1994, Dawson 2007):

Statistické objektivní funkce: 1) Absolutní maximální chyba (Absolut Maximum Error - AME)

(

A M E = m ax Q i − Qˆ i kde

),

Qi…měřená hodnota i-tého pořadí

Qˆ i …simulovaná hodnota i-tého pořadí 2) Vrcholový rozdíl (Peak Difference- PDIFF)

( )

PDIFF = max ( Qi ) − max Qˆ i , kde

Qi…měřená hodnota i-tého pořadí

Qˆ i …simulovaná hodnota i-tého pořadí 3) Střední absolutní chyba (Mean Absolute Error - MAE) 1 n MAE = ∑ Qi - Qˆ i , n i =1 kde

n…počet měření Qi…měřená hodnota i-tého pořadí

Qˆ i …simulovaná hodnota i-tého pořadí 4) Střední chyba (Mean Error - ME)

- 47 -

ME = kde

(

1 n Qi − Qˆ i ∑ n i =1

)

n…počet měření Qi…měřená hodnota i-tého pořadí

Qˆ i …simulovaná hodnota i-tého pořadí 5) Směrodatná odchylka (Mean Square Error – MSE)

(

1 n ∑ Qi − Qˆi n i =1

MSE = kde

)

2

,

n…počet měření Qi…měřená hodnota i-tého pořadí

Qˆ i …simulovaná hodnota i-tého pořadí 6) Relativní absolutní chyba (Relative Absolute Error - RAE) n

RAE =

∑ Q − Qˆ i

i =1 n

∑ Q −Q

,

i

i =1

kde

i

n…počet měření Qi…měřená hodnota i-tého pořadí

Qˆ i …simulovaná hodnota i-tého pořadí Q …průměr měřených hodnot

7) Korelační koeficient (Correlation Coefficient - CC)

∑ ( Q − Q ) ( Q − Qɶ ) ⌢

n

CC =

i =1

i

∑ (Q − Q ) ∑ ( n

i =1

kde

i

2

i

n

i =1

⌢ Qi − Qɶ

)

, 2

n…počet měření Qi…měřená hodnota i-tého pořadí

- 48 -

Qˆ i …simulovaná hodnota i-tého pořadí Q …průměr měřených hodnot Qɶ …průměr simulovaných hodnot

Hydrologické objektivní funkce: 8) Výkonnostní koeficient (Coefficient of Efficiency – CE)

∑ (Q − Q ) ⌢

n

CE = 1 −

i =1 n

∑ (Q − Q )

, 2

i

i =1

kde

2

i

n…počet měření Qi…měřená hodnota i-tého pořadí

Qˆ i …simulovaná hodnota i-tého pořadí Q …průměr měřených hodnot 9) Chyba odhadu kulminačního průtoku (Error of Peak Discharge - PQp) EQ p =

kde

⌢ Qp − Qp

Qp

×100% ,

⌢ Q p …simulovaný kulminační průtok

Q p …měřený kulminační průtok 10) Chyba celkového objemu (Error of Total Volume - EV) n

EV =

⌢

∑ V − V  i

i =1

n

∑V i =1

kde

i

,

i

n…počet měření ⌢ Vi …simulovaný objem i-tého časového kroku

- 49 -

Vi …měřený objem i-tého časového kroku 11) Chyba v době kulminačního průtoku (Error of Time Peak - ETp)

⌢ ETp = Tp - Tp , kde

⌢ Tp …čas simulovaného kulminačního průtoku Tp…čas měřeného kulminačního průtoku

4.

METODIKA

4.1

Charakteristika povodí Modravy 2 a použitých dat.

4.1.1 Popis povodí V Národním parku (NP) Šumava se v souvislosti se vznikem kůrovcové kalamity dostal do popředí požadavek na detailní popis vlivu různých lesních ekosystémů na hydrologický cyklus v měřítku malých lesních povodí. V roce 1997 byl proveden terénní průzkum, který vedli Křovák a Kuřík (2001) a následný výběr vhodných lokalit pro výstavbu experimentálních povodí na území NP Šumava zaměřených na posouzení vlivu rozdílných typů lesního ekosystému na odtokové poměry pokud možno s co nejpodobnějšími základními fyzicko-geografickými charakteristikami. Proto byla vybrána tři velmi malá povodí nacházející se v blízkosti Modravy, tyto povodí se od sebe odlišují vegetačním pokryvem, jejich vzájemná vzdálenost není delší než 14 km a všechna tři povodí mají severní expozici . Na povodí Modrava 1 se jedná o vegetační pokryv, který je tvořen odumřelým lesním porostem v důsledku silné kůrovcové kalamity, dále se jedná o paseku s 10-15 letými smrky na Modravě 2 a posledním vegetačním pokryvem je zdravý 150 let starý smrkový les s příměsí buku, kde dochází k přirozenému zmlazení. Tento les na Modravě 3 ještě nebyl zasažen kůrovcovou kalamitou (Kuna, Kuřík 1998). Hlavní funkcí nově vybudovaných experimentálních povodí byl popis hydrologické funkce lesních ekosystémů a vliv lesních porostů na vybrané komponenty hydrologického cyklu. Tyto komponenty srážkoodtokového procesu byly měřeny v období bez sněhové pokrývky na povodí. Detailně je studován srážkoodtokový proces zaznamenaný v různém časovém měřítku. Předmětem výzkumu bylo stanovení vzájemných závislostí vybraných komponent srážkoodtokových událostí. Modravská povodí byla vybudována Katedrou vodního hospodářství a Katedrou biotechnických

- 50 -

úprav krajiny FLE ČZU v roce 1998 v rámci výzkumných aktivit grantového projektu VaV 620/6/97 „Obnova biodiverzity a stability lesních ekosystémů v pásmu přirozeného rozšíření smrku na území NP Šumava“. Modrava 2 představuje území, kde jsou prováděna hydrologická měření srážek, odtoku, teploty a vodivosti vody od roku 1998. V místě uzávěrového profilu o souřadnicích 48º 58´27.3´´ severní šířky a 13º 30´42.5´´ východní délky v souřadném systému WGS84 byl vybudován měrný trojúhelníkový Thomsonův přeliv (příloha č.3), z jehož přepadové výšky je vypočítáván průtok. Přepadovou výšku snímá tlakové čidlo každé 2 minuty. Srážky jsou měřeny překlopným srážkoměrem se záchytnou plochou 200 cm2 ve výšce 1 m nad zemským povrchem v časovém kroku 2 min. Srážkoměr je umístěn v bezprostřední blízkosti měrného profilu. Teplota je měřena ve výšce 2 m nad zemským povrchem a konduktivita vody na měrném přelivu v časovém kroku 1 hod. Naměřené hodnoty jsou zaznamenány sběrnou jednotkou od firmy NOEL. Naměřená data jsou dále zpracována do vhodnějších časových kroků, obvykle po 1 hodině a to pro celý hydrologický rok. Dle Bevena (2002) srážková data naměřená v časovém kroku po 1 hodině a podrobnější jsou nutná pokud odtok z povodí rapidně reaguje na srážky. Naměřené hodnoty teploty a konduktivity vody se dále zpracovávají bez dalších úprav rovnou z měření. Na povodí Modravy 2 se zahájilo měření uvedených meteorologických charakteristik od roku 1998. V letech 1998–2005 se měřilo ve vegetačním období po roztaní sněhu, tedy v rozmezí květen/červen–říjen. V roce 2006 se přistoupilo ke kontinuálnímu měření po celý rok s tím, že srážkoměr je v zimních měsících odpojen. Měření se také rozšířilo o základní charakteristiky sněhové pokrývky jako jsou výška sněhu nebo objemová hmotnost sněhu., Následující rok 2007 bylo měření rozšířeno i o kvalitativních charakteristiky sněhové pokrývky. Pro účely zimního měření je povodí vybaveno 19 sněhoměrnými latěmi. V následujících letech se počítá s rozšířením měření o kontinuální záznam hydropedologických charakteristik. Výzkum Katedry vodního hospodářství je v poslední době zaměřen na analýzu souboru srážkoodtokových událostí s cílem navzájem

porovnat

odtokové

odezvy

povodí,

aplikovat

vybrané

epizodní

srážkoodtokové modely při simulaci odtokového procesu na lesních povodích a odvodit obecnější

závislosti

aplikovatelné

při

popisu

vybraných

charakteristik

srážkoodtokového procesu experimentálních povodí při povodňových stavech.

- 51 -

Charakter výsledné povodňové vlny je ovlivněn řadou okolností, jako jsou srážkový úhrn, jeho časové a prostorové rozložení na povodí, geometrické a orografické charakteristiky povodí, charakter povrchu povodí, předchozí vláhové poměry, hydropedologické a hydrogeologické charakteristiky povodí, antropogenní zásahy v povodí a mnoho dalších, proto je podrobné měření všech dostupných charakteristik nezbytné k detailnějšímu porozumění a popsání srážkoodtokového procesu. 4.1.2

Základní geometrické charakteristiky. Experimentální povodí Modrava 2 leží při hranicích České republiky a Německa

ve vrcholové části Šumavy (příloha č.4) konkrétně na severním svahu lokality Malá Mokrůvka v pramenné oblasti Ptačího potoka, jehož hydrologické pořadí je 1-08-01002 a má identifikační číslo 10113444. Tato lokalita se nachází 5 km jižně od Filipovy Huti na hranici s Bavorskem. Místní název lokality je „Medvědí doupě“. Povodí má nejvyšší bod s nadmořskou výšku 1330 m n.m. na hoře Malá Mokrůvka a v místě uzávěrového profilu 1197 m n.m. Průměrná nadmořská výška povodí dosahuje 1267 m n.m. Celková rozloha povodí zaujímá 0.17 km2. Plocha F = 0.17 km2

FL = 0.1 km2 a FP = 0.07 km2

(FL,P….. plocha levého resp. pravého svahu povodí) Modrava 2 se svou velikostí řadí do kategorie velmi malých povodí. Za velmi malá povodí se považují taková povodí drobných vodních toků, na nichž jsou maximální průtoky převážně vyvolány přívalovými dešti a vyznačují se malou rozvinutostí říční sítě. Předpokládá se rovněž, že přívalový déšť zasáhne celou plochu daného povodí. Tok bývá výrazně vyvinut pouze v údolnici (případně údolnice není úplně rozvinutá), přítoky se nepovažují za významné z hlediska vytváření maximálního odtoku v uzávěrového profilu povodí. Na tvorbě maximálního odtoku se rozhodující měrou podílí svahový odtok (Hrádek,Kuřík 2002). Míra asymetrie vyjádřená pomocí součinitele a [-]:

a =

0.1 − 0.07 FL − FP F −F = L P = = 0.18 FL + FP F 0.17

Orografická rozvodnice dosahuje délky 2.2 km. Tvar povodí společně se sklonovými poměry povodí ovlivňuje dobu soustřeďování povrchového odtoku z povodí do uzávěrového profilu. Nejpoužívanější

- 52 -

charakteristikou tvaru povodí je součinitel tvaru povodí α , který vyjadřuje poměr mezi střední šířkou povodí B a délkou údolnice LU, kde se délkou údolnice rozumí údolí hlavního toku prodloužené až na rozvodnici.

α=

B Lu

Toto vyjádření vychází z přiblížení tvaru povodí na obdélník, jehož plocha je rovna ploše F, strany tohoto obdélníka jsou střední šířka povodí B a údolnice Lu. Vzhledem k tvaru povodí byl k idealizování tvaru povodí použit trojúhelník. B=

α=

2F 0.17 = = 0.453 km Lu 0.75

F B 0.453 = = = 0.302 2 L u 2Lu 2* 0.75

Schematizace tvaru povodí náhradními geometrickými obrazci se používá při popisu povodí a při analýze povrchového odtoku, kdy se zohledňuje vliv tvaru povodí na soustřeďování vody v uzávěrovém profilu (Hrádek, Kuřík 2002). Vypočtená hodnota součinitele α řadí povodí mezi povodí vějířovitého tvaru. Mezi podstatné geometrické charakteristiky patří sklonové poměry povodí. Díky nejednotným sklonovým poměrům je snaha o zjištění středního sklonu svahů v povodí. K základním vypovídajícím charakteristikám sklonových poměrů patří absolutní spád povodí ∆H, který představuje rozdíl mezi minimální a maximální nadmořskou výškou v povodí.

∆H = Hmax - Hmin [m]

∆H = 1330 – 1197

∆H = 133 m

Sklonové poměry v povodí jsou vždy velmi různorodé a je tedy účelné je nějakým způsobem zobecnit. Pro tento účel je vhodný zjednodušený vztah pro střední sklon svahů má tuto podobu: Isv =

H max − H min F

.100 [%]

Isv =

1330 − 1197 . 100 170000

Isv = 32.2 %

Vzhledem k tomu, že průměrný sklon svahů v povodí značně ovlivňuje výpočet rychlosti soustřeďování vody v povodí, doporučuje se stanovit střední sklon svahů v povodí nejlépe vztahem vyjádřeným dle Herbsta: Isv =

△h.∑ lsi F

- 53 -

.100 [%]

∆h……v tomto tvaru je konst., zvolený výškový interval mezi vrstevnicemi [m] lsi…….průměrná délka vrstevnic v i-tém intervalu [m] F…….. plocha povodí [m2] Průměrný sklon svahu dle Herbsta dosahuje 0.21, tedy 21 %. Další důležitou charakteristikou sklonových poměrů, obzvláště na velmi malých a malých povodí je průměrný sklon údolnice Iu: Isv =

H max,u − H min,u Lu

[%] =

1270 − 1197 .100 = 9.7 % 750

H max,u …maximální nadmořská výška údolnice (na rozvodnici) Hmin,u …minimální nadmořská výška údolnice (uzávěrový profil povodí)

4.1.3 Geologické poměry povodí Povodí Mokrůvky spadá dle regionálně geologického dělení Českého masívu do moldanubické oblasti. Moldanubikum je tvořeno metamorfity ve vysokém stupni metamorfózy, které jsou prostupovány plutonickými horninami, jež intrudovali v závěru variské orogeneze. Moldanubické horniny jsou překryty kvartérními uloženinami fluvio-deluviálním a periglaciálního původu. V povodí Mokrůvky jsou moldanubické horniny jednotvárné skupiny proterozoického až paleozoického stáří. Zastoupené horniny jsou silimanit a magmatické horniny. Granity jsou eisgarnského typu. Kvartérní uloženiny lze rozdělit do dvou skupin: fluviální usazeniny Mokrůvky a periglaciální usazeniny (sutě). Fluviální sedimenty jsou písčitošterkovité uloženiny s malým rozsahem vázané na vodoteče Mokrůvky. Periglaciální uloženiny se zde vyskytují jako blokovité sutě horninového podloží s malým transportem v předpolí kontinentálního ledovce. Na blocích jsou patrné mrazové praskliny. Sněžná čára v době posledního kontinentálního zalednění byla v úrovni 1000 - 1100 m n.m. Strukturně spadá sledovaná oblast do antiklinálního pásma Kvildy, která začíná na státní hranici u Mokrůvky a vede přes Černou horu na Františkov a Borovou Ladu. Významnější tektonické projevy nejsou v zájmovém území pozorovány, v ose vodoteče Mokrůvky lze očekávat lokální tektonickou zónu (Levý 2008). Vzhledem, ke skutečnosti, že se jedná o experimentální povodí a je něm snaha zkoumat a vymezit co nejpodrobněji jednotlivé složky odtoku, bylo ve vymezeném

- 54 -

prostoru na terénně vytipovaných lineárních profilech (příloha č.5) provedeno firmou Inset s.r.o. geofyzikální měření za účelem vymezení rámcové hydrogeologické rozvodnice se zaměřením zejména na okolí uzávěrového profilu. Cílem tohoto měření v okolí uzávěrového profilu bylo vysledování průběhu skalního podloží s velmi malou puklinovou propustností, vyhledání poruchových zón v podloží, které jsou syceny povrchovými srážkami a vymezení lokálních depresních struktur s mocnější akumulací pokryvných útvarů, jako zón možné lokální akumulace srážkových vod, případně úseků, kterými může být srážková voda odváděna mimo měřený profil (Levý 2008). Pro geologický průzkum k určení reliéfu skalního podloží byly zvoleny dvě metody a to metoda mělkou refrakční seismikou (MRS) spolu s metodou elektrického odporu v multielektrodovém uspořádání (MOS). Měrný profil je umístěn v erozním žlabu potoka Mokrůvka. Pokryvné vrstvy se na obou březích pohybují od 0.5–5 m. U pokryvu v celém rozsahu měřeného profilu je nutno předpokládat vysokou propustnost prostředí. Na západním břehu měrného profilu se vyskytuje výrazná elevace horninového masivu, kde se vyskytuje výrazná porucha v podloží. V místech profilu P4 naznačuje metoda (MOS) značné porušení příčnými poruchami. Na základě současného geofyzikálního průzkumu povodí tj. bez zjištění plošného rozsahu struktur uvedených výše, nelze říci zda se orografická rozvodnice shoduje s rozvodnicí hydrogeologickou. Je tudíž možné, že do měrného profilu přitéká voda poruchovými strukturami nebo doprovodnými depresemi skalního podloží z širšího okolí západního svahu nebo naopak část srážek je odvedena mimo měrný profil. Podmínky proudění podzemní vody jsou rozdílné v různých geologických jednotkách a hydrogeologických strukturách. V horninách krystalinika je pohyb podzemní vody komplikovaný. Nejvýraznější pohyb je v povrchové části masivu postiženého zvětráním a rozpukáním. V důsledku postižení masivu tektonikou jsou značné rozdíly v proudění podzemní vody i ve vlastním masivu. Pro moldanubické horniny v zóně zvětrávání, tj. cca 20-30 m od povrchu, odpovídá koeficient filtrace řádu 10-5 m.s-1 s transmisivitou 10-4-10-5 m2.s-1. V poruchových zónách masivu pak předpokládáme hodnoty koeficientu filtrace řádu 10-210-3 m.s-1. Povrchovou zónu kamenitých a balvanitých sutí je možné označit jako zónu intenzivního proudění podzemní vody s koeficientem filtrace řádu 10-1 až 10-3 m.s-1. (Levý 2008).

- 55 -

4.1.4 Půdní a vegetační kryt povodí Půdní horizont experimentálního povodí je tvořen mělkým skeletem, ve kterém jsou promíchány větve a kořeny v různé fázi rozkladu. Převládá půdní typ podzol a kryptopodzol (Hudečková 2008). Pod ním je zvětralá granitová matečná hornina, která místy vystupuje na povrch. Hloubka půdního profilu je 0,6–0,8 (1) m. Hodnota ustálené infiltrační rychlosti se pohybuje od 0.073 mm/s, kterou svými měřeními stanovil Jačka (2009) do hodnoty od Kudrnové (2007) 0.299 mm/s. Nasycená hydraulická vodivost byla naměřena 3.46x10-6, což je hodnota, kterou se dá označit málo propustná půda (Jačka 2009). Oproti tomu Hudečková (2008) uvádí hodnotu K = 1.43x10-4. Tedy výrazně vyšší. Jačka (2009) se domnívá, že to může být způsobeno rozdílným způsobem odběru, skladováním, sycením vzorku, odběrem v jiných hloubkách a celkově jinými hydraulickými vlastnostmi vzhledem k heterogenitě povodí. Po kůrovcové kalamitě byla v této lokalitě povolena těžba napadeného smrkového porostu. Původní smrkový porost byl starý přibližně 160 let a na části plochy se vyskytoval porost starý 26 let. Porost rovnoměrně pokrýval celou plochu povodí. Po těžbě byla paseka zalesněna smrkem, částečně jeřábem a javorem klenem. V okolí měrného přelivu se nachází oplocenka. V současné době tvoří povrch terénu vysázené a náletové dřeviny a travní porost. V bylinném patru převažuje třtina chloupkatá (Calamagrostis villosa), bika lesní (Luzula sylvatica), metlička křivolaká (Avenella flexuosa), častá je borůvka (Vaccinium myrtillus), ale i dřípatka horská (Soldanella montana) nebo maliník (Rubus idaeus). Výška bylinného patra činí asi 30 cm (Tesař 2003). V současné době se dá lesní pokryv povodí označit jako paseka s mladým 10-15 letým porostem s převládajícím zastoupením smrku. 4.1.5 Posouzení srážkových dat Srážky vznikají kondenzací nebo desublimací vodní páry v ovzduší, na povrchu území, rostlin a předmětů. Svým objemem určují odtok v povodí v místě uzávěrového profilu a ovlivňují proces povrchového odtoku nejen v bezprostředně následujícím období po dopadnutí srážek na povrch, ale mohou ovlivnit odtok v následujícím roce nebo i v delším časovém období (podpovrchový odtok). V našich podmínkách jsou podstatné srážky atmosférické, pouze na povodích v horských oblastech připadají v úvahu i srážky horizontální, vznikající z mlhy a z mraků pohybujících se po terénu. Podrobné rozdělení atmosférických srážek uvádějí Hrádek s Kuříkem (2002). - 56 -

Na území naší republiky se projevuje kontinentální typ ročního režimu srážek. Tento typ se vyznačuje vysokými srážkovými úhrny v letním období a malými srážkovými úhrny v zimním období. Tento typ ročního režimu srážek je typický převážně pro klima v monzunových oblastech a v pásmech mírného klimatu na pevninách (Hrádek, Kuřík 2002).V našich podmínkách se vyskytuje tento kontinentální typ s jedním maximem (převážně v červenci) a jedním minimem (zpravidla v lednu a v únoru). V horských oblastech na severu Čech vyniká druhé maximum v lednu. Absolutní maximum denního úhrnu srážek (Hs,d ) bylo v České republice zaznamenáno 29.07.1897 na stanici Bedřichov-Nová louka a dosáhlo hodnoty 345 mm. Nejvyšší hodnota dlouhodobého průměrného ročního srážkového úhrnu (Hs,ra) je 1705 mm na povodí Bílého potoka v Jizerských horách. Naopak nejsušší oblast je na Žatecku, kde Hs,ra dosahuje hodnoty 410 mm. V závislosti na Hs,ra a dlouhodobé průměrné roční teplotě vzduchu T byly v České republice vymezeny oblasti sucha a vlhké oblasti podle Langova dešťového faktoru f : f =

Hs, ra T

Hs,ra [mm]

T[Cº]

Pro vyhodnocení Langova dešťového faktoru f je definováno 5 oblastí (tab. 3). půdotvorný Oblast

f

klasifikace oblastí

proces

1

60

velmi suchá

černozemní

2

61-70

suchá

hnědozemné

3

71-80

normální

-

4

81-100

vlhká

podzolový

5

nad 100

velmi vlhká

-

Tabulka č. 3. Vymezení klimatických oblastí dle Langova faktoru

Pro určení Langova dešťového faktoru byly použity hodnoty dlouhodobých ročních průměrů z období 1961–1990 z meteorologické stanice Churáňov, která je spravována ČHMU ležící v těsné blízkosti povodí Modravy 2 a dlouhodobé půlroční průměry srážkových úhrnů a teplot z období 1998–2008 naměřených na povodí Modravy 2. Z těchto 10 let ovšem nejsou k dispozici datové záznamy po celý rok, proto byly pro výpočet zvoleny pouze měsíce, které jsou ve většině sledovaných let k dispozici, konkrétně se jedná o měsíce květen až říjen. V těchto měsících byly v jednotlivých letech vypočteny měsíční srážkové úhrny a měsíční teplotní průměry.

- 57 -

Pro období od května do října byly vypočteny pomocí aritmetického průměru z těchto hodnot dlouhodobé průměrné půlroční srážkové úhrny a dlouhodobé průměrné půlroční teploty. Rok 1998 nebyl zařazen do výpočtu dlouhodobého půlročního průměru teploty a srážkového úhrnu, z důvodu velmi krátké časové řady, kdy byly potřebné veličiny měřeny. V letech 1999–2001 nejsou k dispozici data pro určení srážkového úhrnu a teplot z měsíce května. Přehled časového výskytu událostí se (tab. 4).

1998 1999

Rozsah měření 1.8-9.11 9.6-24.10

2000

8.6-9.11

2001 2002

7.6-15.11 15.5-31.10

2003

24.5-4.11

2004

7.5-21.10

Rok

Výskyt události 27.10 - 28.10 8.7 - 10.7 14.7 - 15.7 28.7 - 29.7 10.8 - 11.8 17.10 - 19.10 9.10 - 12.10 7.10 - 9.10 23.9 - 25.9

2005

8.5-20.10

9.7 - 11.7 16.8 - 18.8 22.8 - 23.8.

2006 2007 2008

10.4-31.12 1.1-31.12 1.1.-31.12

6.8 - 8.8 28.9 - 30.9 3.10 - 4.10

Tabulka č.4. Časový výskyt srážko-odtokových událostí

Řada měřených ročních srážkových úhrnů na určitém místě umožňuje vyhodnotit časový průběh srážek pouze v daném místě. Pro okolí meteorologické stanice Churáňov pro Hs,ra má dešťový faktor f hodnotu:1090.7/4.2 = 260 a jedná se tedy o velmi vlhkou oblast. Pro okolí meteorologické stanice Churáňov pro Hs,květen - říjen má faktor f hodnotu: 600/9.85 = 61. .

Vzhledem k použití pouze dlouhodobého půlročního srážkového úhrnu a

dlouhodobé půlroční teploty je výsledná hodnota faktoru f vynásobena dvěma, aby se dala zařadit do Langovy stupnice klimatických oblastí, která je odvozena pro dané roční charakteristiky. Dlouhodobý půlroční srážkový úhrn z dostupných dat se pro dané měsíce zdá být dostatečně reprezentativním zastoupením srážkových úhrnů v průběhu roku a dá se tak zhruba považovat za poloviční hodnotu dlouhodobého srážkového úhrnu. V České republice byl použit Langův dešťový faktor získaný z dat vegetačního - 58 -

období (duben- říjen) k charakteristice klimatického sucha pro jednotlivé roky (Sobíšek 1993). V tomto případě spadá hodnota 61 také do velmi vlhké oblasti. Pro oblast povodí Modravy 2 pro Hs,květen - říjen byl vypočten faktor o hodnotě:698.5/10.9 = 64.1. Opět s ohledem na fakt, že byl stanoven půlroční srážkový úhrn je dešťový faktor dvojnásobně zvětšen a řadí se tak opět do velmi vlhké oblasti. Toto porovnání Langova faktoru neprokázalo výraznou odlišnost datového souboru naměřeného za 10 let na Modravě 2 od meteorologických dat sbíraných 30 let na stanici Churáňov a napomohlo do jisté míry odpovědět na otázku, zda tento časový úsek dat nepatří do klimatického výkyvu projevujícího se buď výrazně sušším či vlhčím charakterem. Pro posouzení naměřených dat z povodí Modravy 2 z jednotlivých let z hlediska vydatnosti srážkových úhrnů slouží odchylka ročních srážkových úhrnů (Hs,r) od Hs,ra na daném místě a je významná nejen pro dané povodí, ale charakterizuje i tendenci možných proměn časového rozdělení srážek na velkých územních oblastech, pro které je daná meteorologická stanice reprezentativní. Tato teorie je obhajována tím, že odchylka ∆ = Hs,ra - Hs,r je závislá především na všeobecné cirkulaci atmosféry a vliv ostatních orografických faktorů je eliminován, jelikož jsou v dané situaci stejné. Na základě porovnání součtové čáry odchylek

∑∆

v jednotlivých letech z dané

srážkoměrné stanice můžeme posoudit rozdílnost srážkových úhrnů v dané oblasti (Hrádek, Kuřík 2002). Součtová čára odchylek pro povodí Modrava 2 za období 1999 2008, která je sestrojena z rozdílu Hs,ba květen – říjen a Hs,r květen – říjen (obr. 19)

- 59 -

Obrázek č.19.

Ze součtové čáry je dobře patrné, že mezi srážkově vydatnější roky patří 2002, 2004 - 2007. Vzhledem ke skutečnosti, že ve srážkových úhrnech pro roky 1999 – 2001 není započítán měsíc květen, dají se považovat srážkové úhrny jednotlivých let za vcelku vyrovnané, tedy bez patrných tendencí ke změně srážkových úhrnů v dané oblasti. 4.2

Výběr srážko-odtokových událostí Výběr srážko-odtokových událostí z desetileté řady dat z období 1998-2008, je

podřízen metodě API pro odhadnutí počátečního stavu vlhkosti. Použitelnost této metody platí pro s-o události pro vegetační období, kdy není půda zmrzlá a kdy nedochází k tání sněhu a je vhodná pro výrazné odtokové vlny. Vybrané úseky srážkoodtokových epizod musely dále splňovat, aby povodňová vlna měla jednoduchý průběh tj. aby měla jeden vrchol a v optimálním případě byl počáteční a koncový průtok vlny stejný. K dalšímu kritériu výběru patřila z hydrogramu zřetelně identifikovatelná efektivní srážka, která danou vlnu vyvolala (obr. 20), na kterém je ukázka jednoho z hydrogramů. Modrá spojitá čára představuje průtok Q [dm3/s] a svislé červené čáry reprezentují srážky P [mm] obě charakteristiky mají společnou časovou řadu po hodinovém kroku.

- 60 -

Obrázek č.20. Ukázka hydrogramu

Z hlediska datové dostupnosti bylo nutné, aby pro danou událost byla kompletní časová řada srážek a odtoků. Navíc kvůli metodě API bylo nezbytné mít k dispozici i úplnou časovou řadu denních úhrnů srážek a denních teplot pro 30 předchozích dní. Snahou též bylo zastoupení každého roku nejméně jednou s-o událostí. Celkový počet událostí byl stanoven na 15. Horáček (2006) použil pro potřeby kalibrace a verifikace Nashova a Diskinova modelu 14 s-o epizod a zejména Beven (2002) uvádí, že pro jednoduchý model, pro který se odhaduje 4 až 5 parametrů vyžaduje nejméně 15 až 20 hydrogramů pro dostatečně robustní kalibraci. 4.3

Stanovení API Pro každou vybranou s-o událost byl vypočten index předchozích srážek pro 5

(API 5) a 30 dní (API 30) s evapotranspirační konstantou, která je s hodnotou C = 0.93 běžná pro naše podmínky. Výpočet indexu byl zvolen tak, že poslední den, pro který se hodnota určuje, je předchozí kalendářní den před výskytem vybrané události. Jelikož jsou srážková data zaznamenávána v hodinovém kroku byly srážkové úhrny agregovány do kroku hodinového. Stanovení API (0) [mm], kterým se výpočet indexu zahajuje, předcházelo rozdělení celé časové řady srážkových dat do 10 denních úseků. Pro každý úsek byl vypočten průměrný srážkový úhrn a konečně hodnota API (0) byla získána aritmetickým zprůměrováním všech těchto desetidenních úseků.

- 61 -

Počet dní bez srážek mezi API (0) a API (n) představuje t (28) a je uvažováno, až do dne, kdy se vyskytne srážka. Daná srážka navýší hodnotu API v den, kdy ke srážce došlo. Srážka ≤ 0.2 mm je při výpočtu indexu brána jako den bez srážkového úhrnu. 4.4

Závislost počátečních ztrát na API Odhad stavu vlhkosti půdy povodí Modravy 2 je založen na nalezení vztahu

mezi efektivní srážkou a stavem půdní vlhkostí na počátku srážkoodtokové události. Objem srážky, který je ekvivalentní vytvářenému odtoku se označuje jako efektivní srážka a bývá silně nelineární v závislosti na předchozích podmínkách na povodí (Beven 2002). Počáteční obsah půdní vody má přímí vliv na infiltrační kapacitu a následkem toho na povrchový odtok. Počáteční ztráta je definována jako objem srážky [mm], který je zadržen povodím, než se začne tvořit povodňová vlna. Tyto počáteční ztráty může též chápat také jako kapacitu půdního profilu, kdy stupeň nasycení nenasycené zóny přesáhne své reziduální hodnoty. Tento objem je určen ze souvislé srážky před začátkem průtokové vlny. Počáteční ztráty jsou vyjádřeny jako závislost objemu srážky před vzestupnou větví hydrogramu na API. Dále byla použita pro odhad počátečních ztrát závislost API v kombinaci s teplotou na objemu srážky. V (obr. 21) je v černém rámečku vymezena srážka představující počáteční ztrátu.

Obrázek č.21. Hydrogram s vymezenou srážkou, která představuje počáteční ztráty

Za účelem zjištění možné závislosti mezi počátečními ztrátami a vlivu API v kombinaci s teplotou byly k hodnotám API (30) a API (5) přičteny průměrné hodnoty - 62 -

teplot pro příslušná období. Tedy teplotní průměr za 5 (T 5) a 30 dní (T 30) pro každou vybranou srážkovou epizodu. Průměrná denní teplota byla počítána z hodinových měření vztahem dle Klabzuby (2001): (t 7 + t14 + t 21× 2) 4 Kde t představuje změřenou teplotu pro danou hodinu v místním slunečním čas, který přibližně odpovídá středoevropskému času. Byla tedy zjišťována závislost mezi počátečními ztrátami (PZ) a API (5). Též závislost PZ na API (5) společně s T (5). Dále byl sledován funkční vztah mezi PZ a API (30) a též možný vztah PZ na API (30) + T (30). Zjištěné hodnoty samotných API a i v kombinaci s teplotou pro jednotlivé s-o události byly vyneseny do grafu, kde osu x představují počáteční ztráty a na ose y jsou sledované charakteristiky předchozích podmínek na povodí. Takto získané body byly pomocí lineární, mocninné a exponenciální regrese proloženy danými trendy a získány tak rovnice popisující vztah mezi počátečními ztrátami a API. Vybraným ukazatelem, toho jak dobře rovnice získané regresní analýzou popisují vztahy mezi sledovanými veličinami, byl zvolen koeficient determinace (R2). R2 poskytuje informaci o tom, jakou část z celkového součtu čtverců odchylek se nám podařilo proložením regresní funkce vysvětlit. Představuje podíl regresního součtu čtverců (SR) a celkový součet čtverců (ST).

R2 = SR / ST SR =

∑ ( yˆ

i

− y)2

i

yˆi je y-souřadnice i-tého bodu na přímce obdržené metodou nejmenších čtverců.

ST = ∑ ( yi − y )

2

i

Pro získání modelu počátečních ztrát byly z rovnic, získaných z regresní analýzy vypočteny inverzní funkce. Z těchto modelových rovnic se pro konkrétní hodnoty API a teploty získaly objemy počátečních ztrát.

- 63 -

4.5

Příprava dat pro kalibraci modelu Boussmo Z důvodů kalibrace a validace modelu se vybrané s-o epizody rozdělily na dva

soubory a to na soubor kalibrační o 7 událostech a validační soubor o 8 událostech. Každá s-o epizoda byla rozdělena na dva textové soubory. Jeden obsahoval srážky a druhý odtoky. Obě veličiny byly seřazeny dle referenčního času, jehož časový krok byl po jedné hodině. Kvůli přizpůsobení se modelu byl referenční hodinový čas nahrazen 5 minutovým časovým krokem. Interpolace časového kroku byla provedena z důvodu výpočtu a na výsledek kalibrace by neměla mít vliv, neboť došlo pouze ke zmenšení časového kroku. Počáteční ztráta je do modelu zakomponována tak, že srážky jednotlivých událostí byly zmenšeny o objem počátečních ztrát tzn. jednotlivé srážky v počátečních časových krocích byly nulovány do doby, než takto došlo k vyčerpání hodnoty počáteční ztráty. Následně bylo nutné takto již upravené srážky převést z hodinových, respektive 300 sekundových úhrnů, na intenzitu srážek [mm/s] a společně s odtoky [dm3/s] převést do jednotek SI. Dalším krokem bylo stanovení celkového odtokového součinitele, který dále slouží jako objektivní kritérium pro automatickou optimalizaci. Odtokový součinitel (OS) je definován jako podíl objemu odtoku (Hq) a objemu srážek (Is). Výška odtoku představuje podíl součtu všech objemů průtoků (Qi) za danou událost a plochy povodí (F) [m2]. n

Hq OS = = Is

∑Q

F

i

i =1

n

∑ Is i =1

(30)

i

Závěrečným bodem příprav dat pro kalibraci představovalo spojení všech srážek i odtoků s-o epizod do jedné časové události. Bylo tak učiněno z důvodů výpočetních potřeb algoritmu. Ze stejných důvodů byl pro kalibraci modelu připraven textový soubor obsahující časové konce srážek. V průběhu přípravy dat byla snaha o co největší automatizaci jednotlivých kroků, aby se tak předešlo možným hrubým chybám při přílišné manipulaci s daty. Datové přípravy probíhaly v tabulkovém procesoru, v kterém byla data poskytnuta a ve statistickém programu R. Ukázka skriptů usnadňujících přípravu dat v příloze č. 6.

- 64 -

4.6

Automatická optimalizace Pro nakalibrování 5 z 16 parametrů modelu byla použita automatická

multikriteriální optimalizace založená na genetickém algoritmu SADE (Simplified Atavistic Differential Evolution) vylepšeném dle (Kučerová 2004). Tato multikriteriální automatická optimalizace hledá pomocí algoritmu maximum z (funkce - ε ). Kde ε (epsilon) představuje více objektivní funkci. Tato funkce je definována jako :

ε =a+b+c a = 100 * (hodnota sumy čtverců vzdáleností měřených a kalibrovaných dat) b = 50 * (rozdíl sumy odtokových koeficientů kalibrovaných a měřených dat) c= 50 * (rozdíl maximální hodnoty největší kalibrované a měřené srážky) Soubor parametrů, jehož výsledná hodnota ε je nejmenší, je považován za optimální soubor a takový to soubor parametrů je dále použit pro účely validace. Touto automatickou optimalizací byly získány tři soubory parametrů pro jednotlivé odhady počátečních ztrát.

5.

VÝSLEDKY

5.1

API a odhad počátečních ztrát Zahajovací hodnota API (0) byla výpočtem stanovena na 4.14 mm a stala se

základem pro určení předchozích srážkových indexů API (5) a API (30) pro všechny s-o události. Přehled výsledných indexů, průměrných teplot a objemů srážek před povodňovou vlnou pro s-o epizody podává (tab. 5). S-o události jsou značeny dle klíče M2026950, kde M2 je označení povodí Modrava 2, dále 02 představuje rok výskytu a poslední čtyřčíslí je referenční čas, odkdy daná událost začala být sledována.

- 65 -

S-O epizoda M2987197 M2994535 M2004685 M2005017 M2015305 M2026950 M2036702 M2036757 M2046380 M2054558 M2055432 M2055597 M2065222 M2076490 M2086627

API (5) [mm] 34.39 13.68 32.33 27.5 24.29 9.1 66.93 97.92 8.68 62.48 17.88 22.57 40.5 15.12 48.28

API (30) [mm] 56.39 44.44 63.85 60.75 44.31 49.29 71.84 108.21 25.09 96.99 62.65 85.19 65.76 66.24 67.64

T (5) [ºC] 4.0 17.0 8.0 14.0 12.9 5.9 7.0 3.0 12.0 11.0 10.0 13.0 11.0 8.3 5.0

T (30) [ºC] 5.5 12.0 12.0 11.0 14.0 4.5 10.0 7.0 9.0 13.0 13.0 12.0 17.0 7.0 7.0

V [mm] 31 20.8 20.6 16.8 23.6 15.6 12.8 37.2 59.8 9.8 14.2 8.4 13.4 29.6 17.8

Tabulka č.5. Přehled předchozích srážkových indexů, teplot a objemů před povodňovou vlnou.

Po vynesení API a API + teplota (hodnoty API(n) + T(n) na počátečních ztrátách V [mm] do grafu, byla ze souboru s-o událostí vyloučena M2036757 pro své hodnoty, které se velmi odlišovali od ostatních a výrazně tak ovlivňovala výslednou závislost posuzovanou koeficientem determinace. Největší hodnoty R2 bylo získáno při modelování funkčního předpisu za použití API (30) + T (30) pro všechny typy regresní analýzy (obr. 22, 23, 24). Pro další typy API je výsledek regresní analýzy v příloze č. 7.

- 66 -

Pro lineární regresní analýzu:

API 30 + T30 140.0 hodnota

120.0

Lineární (hodnota)

API [mm]

100.0 80.0 y = -1.1071x + 95.222

60.0

2

R = 0.5822

40.0 20.0 0.0 0

10

20

30

40

50

60

70

V [mm]

Obrázek č.22

Pro mocninnou regresní analýzu:

API 30 + T30 140.0 hodnota

120.0

Mocninný (hodnota)

API [mm]

100.0 80.0 60.0 -0.4771

y = 278.85x 40.0

2

R = 0.7104

20.0 0.0 0

10

20

30

40

50

60

70

V [mm]

Obrázek č.23

- 67 -

Pro exponenciální regresní analýzu:

API 30 + T30 140.0 hodnota

120.0

Exponenciální (hodnota)

API [mm]

100.0 80.0 60.0 -0.0184x

y = 102.17e

40.0

2

R = 0.6963 20.0 0.0 0

10

20

30

40

50

60

70

V [mm]

Obrázek č.24

Přehled výsledných hodnot R 2 pro jednotlivé typy regresní analýzy zpřehledňuje (tab. 6) Typ API

R2 lineární

API (5) API (30) API (5) +T (5) API (30) +T (30)

0.2062 0.5549 0.2258 0.5822

R2 mocninná R2 exponenciální 0.2527 0.6743 0.2359 0.7104

0.2739 0.6855 0.2231 0.6963

Tabulka č.6. Hodnoty koeficientu determinace.

Pro typ regresní rovnice, která dle R2 nejlépe popisuje průběh závislosti, tedy pro mocninnou rovnici s určením předchozího stavu povodí pomocí API (30) + T (30) se odvodila inverzní funkce:

 x  y = −0.4771    278.85 

- 68 -

Pro posouzení zahrnutého vlivu teploty pro odhad počátečních podmínek na model Boussmo se dále určila inverzní funkce u mocninného modelu pro API (30). Rovnice má podobu:

x   y = −0.5334    277.96  A třetí typ, ze kterého byla určena inverzní funkce byl lineární pro API (30) + T (30) s následujícím tvarem rovnice: y = x - 95.222 / -1.1071 Z těchto rovnic byly vypočteny počáteční ztráty pro každou s-o událost s daným typem API. Přehled konkrétních hodnot počátečních ztrát je uveden (tab. 7)

S-O epizoda

Poč. ztráty přímka API (30)+T(30)

Poč. ztráty mocninná API (30)

Poč. ztráty mocninná API (30)+T( 30)

M2987197 M2994535 M2004685 M2005017 M2015305 M2026950 M2036702 M2036757 M2046380 M2054558 M2055432 M2055597 M2065222 M2076490 M2086627

30.1 35.0 17.5 21.2 33.3 37.5 12.1 -18.1 55.2 -13.3 17.7 -1.8 11.3 19.9 18.6

19.9 31.1 15.8 17.3 31.3 25.6 12.64 5.86 90.8 7.2 16.3 9.2 14.9 14.7 14.2

23.5 28.5 15.3 17.2 26.6 31.5 13.1 6.4 81.9 7.0 15.4 9.1 12.8 16.5 15.8

Tabulka č. 7 . Přehled počátečních ztrát pro jednotlivé typy API Lineární závislost byla zvolena i přes nízký koeficient determinace z důvodu dobré identifikovatelnosti krajních hodnot. Výskyt záporných hodnot vypovídá o skutečnosti, že pro tyto případy není odhad počáteční vlhkosti povodí pomocí lineárního vztahu mezi PZ a API vhodný. Nicméně je tento odhad dále v práci používán pro ostatní nenulové hodnoty. Záporné hodnoty jsou brány jako nulové.

- 69 -

5.2

Kalibrace Rozdělení s-o epizod na soubor kalibrační a validační udává (tab. 8): Validační soubor M2994535 M2005017 M2015305 M2026950 M2036757 M2054558 M2055597 M2076490

Kalibrační soubor M2987197 M2004685 M2036702 M2046380 M2055432 M2065222 M2086627 ---

Tabulka č.8. Rozděleni s-o epizod na kalibrační a validační soubory.

Hodnoty pěti kalibrovaných parametrů získaných pomocí automatické optimalizace metodou SADE pro jednotlivé způsoby určení počátečních ztrát (tab. 9). Seznam pěti kalibrovaných parametrů je následující: K……...Hydraulická vodivost [m/s] H………Hloubka nepropustného podloží pod korytem [m] h……….Průtočná hloubka koryta [m] E………Elasticita [-]

α ……...Sklon podloží [-] API(30) API(30)+T(30) Ozn.parametru API(30)+T(30) lineární vztah mocninný vztah mocninný vztah K

3.20E-05

3.48E-05

3.02E-05

H

40

30.6667

25.4614

h

1.68345

1.66772

1.69794

E

1.99E-05

1.91E-05

2.13E-05

0.1

0.100046

0.1

α

Tabulka č.9. Výsledné hodnoty nakalibrovaných parametrů.

- 70 -

Výpočet optimální sady kalibrovaných parametrů probíhal iterativně. Průběh snižování chyby epsilon pro jednotlivé odhady počátečních ztrát (obr. 25, 26, 27). Pro lineární model API(30)+T(30) počátečních ztrát:

Obrázek č. 25

Pro mocninný model API (30) počátečních ztrát:

Obrázek č. 26

- 71 -

Pro mocninný model API (30)+T(30) počátečních ztrát:

Obrázek č.27

V příloze č. 8 se nachází tabulka s hodnotami průběhů konvergence chyby a údaji o počtu proběhlých iterací pro jednotlivé modely odhadů počátečních ztrát.. 5.3

Validace Následná validace pro určení vhodnosti nakalibrovaných parametrů je

vyhodnocena podle grafických výstupů na základě subjektivního posouzení a dle odtokového součinitele, který byl zvolen jako hledisko objektivní. V grafických výstupech (obr. 28-36) jsou vybrány pro každý způsob odhadu PZ takové epizody, jejichž validace přinesla nejbližší shodu mezi průtoky simulovanými modelem Boussmo a s průtoky naměřenými. U s-o události M2015305 došlo při odečtu počátečních ztrát v případě všech odhadů nasycenosti povodí k úplnému vynulování srážek. U s-o události M2026950 se srážky vynulovali v případech, kdy byla použita pro odhad PZ závislost lineární na API (30) + T(30) a mocninná na API (30) + T (30). V takovém případě nedojde samozřejmě žádný srážkový vstup do nasycené zóny a nedojde tak k žádnému odtoku, proto byly tyto události z validace vynechány. Pro každou variantu odhadu nasycenosti povodí bylo provedeno šest simulací. Uspokojivé shody mezi simulovanými a naměřenými průtoky bylo dosaženo zhruba v polovině případů.

- 72 -

Obrázek č.28. Nespojité body jsou naměřený průtok, plná čára představuje simulovaný průtok.

Obrázek č.29 Nespojité body jsou naměřený průtok, plná čára představuje simulovaný průtok.

- 73 -

Obrázek č. 30. Nespojité body jsou naměřený průtok, plná čára představuje simulovaný průtok.

Obrázek č. 31 Nespojité body jsou naměřený průtok, plná čára představuje simulovaný průtok.

- 74 -

Obrázek č.32. Nespojité body jsou naměřený průtok, plná čára představuje simulovaný průtok.

Obrázek č.33 Nespojité body jsou naměřený průtok, plná čára představuje simulovaný průtok.

- 75 -

Obrázek č. 34. Nespojité body jsou naměřený průtok, plná čára představuje simulovaný průtok.

Obrázek č. 35 Nespojité body jsou naměřený průtok, plná čára představuje simulovaný průtok.

- 76 -

Obrázek č.36. Nespojité body jsou naměřený průtok, plná čára představuje simulovaný průtok.

K objektivnímu zhodnocení validace byla použita absolutní hodnota z rozdílu odtokového koeficientu získaného na základě průtoků simulovaných modelem a odtokového součinitele z měřených dat pro danou s-o. Podrobný přehled těchto rozdílů je uveden (tab. 10-12).

Validace 1 2 3 4 5 6

Rozdíl odtokových koeficientů při lineárním API (30)+T(30) Odtokový Odtokový koef. koef. pozorovaný S-o epizoda |OKP - OKM| modelový (OKP) (OKM) M2005017 M2036757 M2054558 M2055597 M2076490 M2994535

0.02521 0.02197 0.02295 0.02337 0.0226 0.022699

0.036545 0.013044 0.00706 0.008982 0.019703 0.016101 arit. průměr =

0.011335 0.008926 0.01589 0.014388 0.002897 0.006598

relativní chyba [%] 44.96 40.63 69.24 61.57 12.82 29.07

0.010005667 Tabulka č. 10

- 77 -

Validace 1 2 3 4 5 6

Rozdíl odtokových koeficientů při mocninném API(30) Odtokový koef. Odtokový koef. S-o epizoda |OKP - OKM| modelový pozorovaný (OKM) (OKP) M2005017 M2036757 M2054558 M2055597 M2076490 M2994535

0.032162 0.028445 0.029837 0.0322 0.029005 0.029941

0.0267095 0.015265 0.0085 0.018211 0.01662 0.013407 arit.průměr =

0.0054525 0.01318 0.021337 0.013989 0.012385 0.016534

relativní chyba [%] 16.95 46.34 71.51 43.44 42.70 55.22

0.013812917 Tabulka č. 11

Validace

Rozdíl odtokových koeficientů při mocninném API(30) + T(30) Odtokový Odtokový koef. koef. S-o epizoda pozorovaný |OKP - OKM| modelový (OKP) (OKM)

1 2 3 4 5 6

M2005017 M2036757 M2054558 M2055597 M2076490 M2994535

0.039162 0.034345 0.036114 0.039468 0.035019 0.036479

0.026544 0.01548 0.008457 0.018429 0.016531 0.120578

0.012618 0.018865 0.027657 0.021039 0.018488 0.084099

arit.průměr =

0.030461

relativní chyba [%] 32.22 54.93 76.58 53.31 52.79 230.54

Tabulka č. 12 Validace, pro které se neprojevila dobrá shoda mezi naměřenými a simulovanými průtoky pro všechny způsoby odhadu půdní vlhkosti jsou umístěny v příloze č. 9.

- 78 -

6. DISKUZE Z hlediska průměrného nejmenšího absolutního rozdílu odtokových součinitelů dopadla validace nejlépe pro lineární závislost PZ na API(30)+T(30), pro kterou je hodnota průměrného absolutního rozdílu až třikrát menší než u API(30)+T(30) a je nižší než u mocninného odhadu API(30). Tento výsledek je překvapující s ohledem na to, že odhad PZ na lineárním vztahu API(30)+T(30) má nejmenší R2 a ne vždy dostatečně popisuje průběh počátečních ztrát. To může být ovšem ovlivněno tím, že při odhadu vlhkosti lineárním API(30)+T(30) se v případě záporných hodnot počátečních ztrát tyto ztráty neuvažují a celková suma srážek je tedy vyšší než v ostatních dvou případech a tím je významně ovlivněna hodnota odtokového součinitele. Proto by pro posouzení vhodnosti validace měla být zvolena jiná objektivní funkce, než rozdíl odtokových součinitelů. Vhodné by bylo například znormovat jak simulované tak změřené odtoky dle nejvyšší hodnoty a vyčíslit a porovnat získaná rezidua. Průměrný absolutní rozdíl odtokových součinitelů je pro mocninný model API(30) a API(30)+ T(30) znatelně odlišný ve prospěch API(30). Otázkou je, do jaké míry toto porovnání ovlivňuje skutečnost, že simulovaná data byla modelována pro různé sady optimálních parametrů. Zejména by mohl být patrný vliv parametru H, tedy hloubky nepropustného podloží pod korytem, který pokud si uchoval svůj fyzikální rozměr je pro každý způsob odhadu vlhkosti výrazně jiný (tab. 9). Ostatní kalibrované parametry se na první pohled od sebe příliš neliší a parametry jako hydraulická vodivost, hloubka pod korytem a sklon svahu se svými hodnotami příliš neodchylují od reálných charakteristik povodí. Bez provedení citlivostní analýzy nemůžeme určit do jaké míry je výsledek validací ovlivněn kalibrovanými parametry. Z provedené validace lze konstatovat, že pro s-o události, které mají výraznou odezvu na srážku a odtok větší než 0.015 [m3/s] se dosahuje uspokojivé shody mezi naměřenými a simulovanými průtoky. Pro události M2054558 (validace č.3) a M2055597 (validace č.4), jejichž simulované průtoky byly výrazně nadhodnocené je společné krom malého průtoku, že mají velmi vysokou hodnotu předchozího srážkového indexu a jejich počáteční ztráty jsou tedy malé tj. leží v oblasti odhadu počátečních ztrát v krajních bodech, kde je určení závislosti nejvíce nepřesné a problematické. Pro tyto krajní body při odhadu vlhkosti půdy by pravděpodobně stálo za úvahu počáteční ztráty v modelu Boussmo vyjádřit pomocí spíše exponenciálního vztahu, před kterým byla upřednostněna lineární funkční vzdálenost. Pro zpřesnění - 79 -

odhadu počátečních ztrát pomocí API by bylo vhodné rozšířit počet pozorovaných s-o událostí a vyzkoušet API pro menší počet předchozích dní, např. API (15).

7.

ZÁVĚR Prokázala se závislost objemu srážky před vzestupnou větví hydrogramu na API.

Počáteční ztráty byly vyjádřeny jako závislost objemu těchto srážek na API. Jako nejvhodnější odhad stavu vlhkosti půdy se ukázal mocninný vztah počátečních ztrát na předchozím srážkovém indexu z předchozích třiceti dní. Pro takto definovanou komponentu odhadu nenasycené zóny je použitelný model Boussmo pro modelování srážko-odtokové události na povodí Modrava 2, která se projeví rychlým nástupem odezvy a s průtokem povodňové vlny větším než 0.015 [m3/s]. Pro zpřesnění by bylo vhodné parametry modelu podrobit citlivostní analýze a znovu kalibrovat pomocí automatické optimalizace při větším počtu iterací pro nalezení globálního optima ze souboru všech parametrů a jako komponentu pro odhad půdní vlhkosti použít mocninný vztah objemu srážky před vzestupnou větví hydrogramu a API (30).

- 80 -

8.

PŘEHLED LITERATURY A POUŽITÝCH ZDROJŮ

AQUALOG, http://www.aqualogic.cz/Slu_by/Produkty/Software/software_0.html ANDERSON, E.A., NEUMAN, P.J., 1984: Inclusion of Frozen Ground Effects in a Flood Forecasting Model. Fifth Northern Research Basins Symposium and Workshop, March 19-23, Vierumaki, Finland. BEVEN K., 2002: Rainfall runoff modeling The Primer. John Wiley and Sons. BLAŽKOVÁ Š., 1993: Srážkoodtokové modelování založené na principu jednotkového hydrogramu. Výzkumný ústav vodohospodářský T.G. Masaryka. Praha. BUCHTELE, J. 1972: Kategorizace povodňového režimu na tocích Vltavské kaskády, Sborník prací hydrometeorologického ústavu v Praze, svazek 18, s. 64 – 139, HMÚ v Praze BURNASH, R., FERRAL, L., 1996: Conceptualization of the Sacramento Soil Moisture Accounting Model. NWSRFS Users Manual, Part II.3, National Weather Service, NOAA, DOC, Silver Spring.. DAŇHELKA J. et kol., 2003: Posouzení vhodnosti srážko-odtokových modelů s ohledem na simulaci povodňových stavů pro lokality na území ČR. ČZU a ČHMU, Praha. DAWSON C. W, ABRAHART R . J., SEE L. M., 2007: HydroTest: A web-based toolbox of evaluation metrics for the standardised assessment of hydrological forecasts. Environmental Modelling & Software 22: 1034-1052. DESCROIX L., NOUVELOT G.F., VAUCLIN M., 2002: Evaluation of an precipitation index to model runoff yield in the western Sierra Madre (North–west Mexico). Journal of Hydrology 263: 114-130. DHI, 2000: Mike SHE, An Integrated Hydrological Modelling System. DHI Water & Environment, Stern Allé 5, DK-2970 Horsholm, Denmark. Nofo print. DIERMANSE F. L. M., 2001: Physically based modelling of rainfall-runoff processes. Ph.D. thesis, Delft University Press, The Netherlands. DINGMAN L.S., 2002: Physical hydrology. Prentice Hall. EXPERIMENTÁLNÍ POVODÍ MODRAVA, http://www.kvhem.cz/vyzkum/povodimodrava/

- 81 -

FEDORA M.A., BESCHTA R.L, 1988: Storm runoff simulation using an antecedent precipitation index (API) model. Journal of Hydrology 112: 121-133. HORÁČEK S., 2006: Application of chosen event based rainfall runoff models on the small mountain catchment Modrava 2. In , 20. 9. 2006 Luxembourg. Centre de Recherche Public - Gabriel Lippmann. Belvaux: Centre de Recherche Public - Gabriel Lippmann. HRÁDEK F. et KUŘÍK, P., 2008: Hydrologie. ČZU, Praha. KUČEROVÁ A., HRSTKA O.,: 2004: Improvements of real coded genetic algorithms based on diffirential operators preventing premature convergence. Advances in Engineering Software 35: 237-246. HUDEČKOVÁ K., 2008: Vyhodnocení hydropedologického průzkumu experimentálního povodí Modrava 2. Diplomová práce. Nepublikováno, Dep.: KVHEM FŽP ČZU, Praha. JAČKA L., 2009: Stanovení vybraných hydropedologických charakteristik na povodí Modrava 2. Diplomová práce. Nepublikováno, Dep.: KVHEM FŽP ČZU, Praha. KINEROS 2, SOFTWARE DOCUMENTATION, http://www.tucson.ars.ag.gov/kineros/ KLABZUBA J., 2001: Aplikovaná meterologie a klimatologie. V.díl, Bilance tepla na aktivní povrchu, teplota půdy, vzduchu a vody. ČZU, Praha. KOHLER M.A, LINSLEY R.K, 1951: Predicting the runoff from storm rainfall. Weather Bureau, US Department of Commerce, Research Paper No.34, Washington. KOVÁŘ P., 1990: Využití hydrologických modelů pro určování maximálních průtoků na malých povodích. Vysoká škola zemědělská v Praze. KUDRNOVÁ P., 2007: Řešení infiltrace pomocí vybraných postupů. Diplomová práce. Nepublikováno, Dep.: KVHEM FŽP ČZU, Praha. KUNA P., KUŘÍK P., 1998: Hydrometeorologická měření na Šumavě. Zprávy lesnického výzkumu č.03-04. KŘOVÁK F., KUŘÍK P., 2000: Vliv lesních ekosystémů na odtokové poměry krajiny. In: MÁNEK J. [ed]: Aktuality šumavského výzkumu, sborník z konference: 75–79.

KURÁŽ M., 2009: The Residency Report of Michal Kuráž at ITCR Costa Rica - the Evaluation of the Saturated Hydraulic Conductivity of the Maritza, Catchment, Guanacaste, and the BOUSSMO Code Documentation – Rainfall Runoff Numerical Model. Nepublikováno, Dep.: IDLI FSV ČVUT, Praha. LEVÝ O., 2008: Geofyzikální průzkum povodí Modrava 2. INSET s.r.o., Nepublikováno, Dep.: KVHEM FŽP ČZU, Praha.

- 82 -

PAVLÁSEK J., 2005 : Aplikace vybraných matematických modelů pro řešení odtoku z povodí. Dizertační práce. Nepublikováno, Dep.: KVH FLE ČZU, Praha.

PAVLÁSEK J., 2008 : Vyhodnocení povodňové události na povodí Modrava 2 z 8.8. 2008. In: MÁCA P., NECHVÁTAL M., KULHAVÝ Z., SOUKUP M. [eds]: Monitoring a vyhodnocení extrémních odtokových poměrů v povodí drobných vodních toků z hlediska prevence a zmírňování povodňových škod. Sborník workshopu grantového projektu NAZV 1G46040, ČZU a VÚMOP, Praha: 10. ŘEDINOVÁ J., 2004: Analýza maximálních průtoků v povodích Modrava 2, Morávka a Spůlka. Diplomová práce. Nepublikováno, KVH FLE ČZU, Praha ŘIČICOVÁ P. ET KOL., 2003: Povodeň v srpnu 2002 v Jizerských horách. VTEI VÚV T. G. M., SIS ČHMÚ. Praha

SMITH M. B ET KOL., 2000: Evaluation of advantages of the Continuous SAC-SMA model over an API model.15th Conference on Hydrology, AMS,January 9-14, 2000, Long Beach, CA NOAA/NWS, Silver Spring, MD 20910. SOBÍŠEK, B. ET KOL., 1993: Meteorologický slovník, výkladový a terminologický. Academia. Praha. TESAŘ, M., M. ŠÍR, E. ZELENKOVÁ, Ľ. LICHNER, 2003: Vodní a teplotní režim lesa, paseky a holiny ve vegetační sezóně, Seminář s mezinárodní účastí „Hydrologie půdy v malém povodí“, Šír, M., Ľ. Lichner, M. Tesař [Eds], Praha. VALENTOVÁ J., 2007: Hydraulika podzemní vody. ČVUT, Praha. YAPO P.O, GUPTA H.V.,SOROOSHIAN S., 1996: Automatic calibration of conceptual rainfall-runoff models:sensitivity to calibration data. Journal of Hydrology 181: 23-48. YU, P. S., LIU, C. L., LEE, T. Y. 1994: Application of transfer function model to a storage-runoff process. In Hipel K.W.: Time Series Analysis in Hydrology and Environmental Engineering. Kluwer academic Publisher, Dordrecht, s.87-97 ZEMAN E., 1994: Hydroinformatika a hydrologické modely. Habilitační práce, ČVUT, Praha.

- 83 -

9. PŘÍLOHY Příloha č.1: Konfigurační soubor modelu Boussmo: # ############# BOUSSMO 1.0beta ######################### # a BOUSSinesque equation solver for MOdrava type catchments, #algorithm to saporate surface runoff is based on material properties # and reservoir approach. The surface runoff is calculated by #kinematic wave equation. # The catchment is expressed by the reservoir, the outflow from #subsurface flow is calculated by Boussinesque # equation, the outflow from surface runoff reservoir is calculated by #the kinematic wave equation. # a program for evaluating the rainfall outflow event on a small #steep mountain catchments # originaly Modrava catchment was suppose to be a model scenario for #this approach # later applied on Maritza catchment, Guanacaste, Costa Rica. # Unsaturated zone is in recent version neglected, appliable for rainy #season.

# # # # # # # # # #

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see .

#-----------------------program configuration-------------------------------

#hydraulic conductivity 3.47743e-05 m/s # depth of nonpermeable layer under the river basin 30.6667 m #depth of river basin 1.66772 m #top boundary condition 0.0 m # slope of the plane 0.100046 #initial unsaturated flow 0 m/s

- 84 -

#lenght of the plane 450.0 m #width of the plane 200 m #grid size for boussinesque (higher precission required) prostorova 0.23 #grid size for kinematic wave (possible lower precission) 1

#effective porosity 1.91258e-05 #null size specific discharge (precission) 0 #manning roughness coefficient of hillslope (SI units) 0.12 #minimal time step for kinematic wave 1e-5 #end time simulation - konec hydrogramu 2.628e4 s #cathment volume 550000000000 m3

- 85 -

Příloha č.2: Porovnání výsledků modelu Sacramento s modelem APIc (Smith 2000).

- 86 -

Příloha č.3: Thomsonův přeliv.

Thomsonův měrný profil na Modravě 2. (Experimentální povodí Modrava)

- 87 -

Příloha č.4: Lokalizace povodí Modrava 2

Poloha experimentálního povodí Modrava 2. (Experimentální povodí Modrava)

Zarůstající paseka na Modravě 2. (Experimentální povodí Modrava)

- 88 -

Příloha č.5: Průzkumné geofyzikální profily (Levý 2008)

Liniové profily pro vytyčení geofyzikálního průzkumu (Levý 2008)

- 89 -

Příloha č.6: Skripty v programu R: Spojení všech srážek do jedné události: a<-read.table('in1.txt') b<-read.table('in2.txt') c<-read.table('in3.txt') d<-read.table('in4.txt') e<-read.table('in5.txt') f<-read.table('in6.txt') g<-read.table('in7.txt') T=length(b[,1]) TT=length(c[,1]) TTT=length(d[,1]) TTTT=length(e[,1]) TTTTT=length(f[,1]) TTTTTT=length(g[,1]) a[length(a[,1])+1,1]=length(a[,1])*300+300 a[length(a[,1]),2] =0 b[1,1]<-30000 b[,1]<-seq(b[1,1],b[1,1]+(T-1)*300,by=300) b[length(b[,1])+1,1]=length(b[,1])*300+30000 b[length(b[,1]),2]=0 c[1,1]<-60000 c[,1]<-seq(c[1,1],c[1,1]+(TT-1)*300,by=300) c[length(c[,1])+1,1]=length(c[,1])*300+60000 c[length(c[,1]),2]=0 d[1,1]<-90000 d[,1]<-seq(d[1,1],d[1,1]+(TTT-1)*300,by=300) d[length(d[,1])+1,1]=length(d[,1])*300+90000 d[length(d[,1]),2]=0 e[1,1]<-120000 e[,1]<-seq(e[1,1],e[1,1]+(TTTT-1)*300,by=300) e[length(e[,1])+1,1]=length(e[,1])*300+120000 e[length(e[,1]),2]=0 f[1,1]<-150000 f[,1]<-seq(f[1,1],f[1,1]+(TTTTT-1)*300,by=300) f[length(f[,1])+1,1]=length(f[,1])*300+150000 f[length(f[,1]),2]=0 g[1,1]<-180000 g[,1]<-seq(g[1,1],g[1,1]+(TTTTTT-1)*300,by=300) g[length(g[,1])+1,1]=length(g[,1])*300+180000 g[length(g[,1]),2]=0

- 90 -

merge<-c(a[,1],b[,1],c[,1],d[,1],e[,1],f[,1],g[,1]) rain<-c(a[,2],b[,2],c[,2],d[,2],e[,2],f[,2],g[,2]) plot(merge,rain,ylim = c(max(rain),0),xlab='time step [300 sec]',ylab='srážky [m/s]',type='h',col='red',axes=FALSE,) axis(2) axis(1) mat<-c(merge,rain) write.table(matrix(mat,nrow=length(merge),ncol=2),row.names=FALSE,col.name s=FALSE,sep= " ",file="4compare.in.txt") Celkový odtokový koeficient: a<-read.table('in1.txt') # pro srazky b<-read.table('in2.txt') c<-read.table('in3.txt') d<-read.table('in4.txt') e<-read.table('in5.txt') f<-read.table('in6.txt') g<-read.table('in7.txt') h=a$V2 hh=sum(h) ch=b$V2 chh=sum(ch) i=c$V2 ii=sum(i) j=d$V2 jj=sum(j) k=e$V2 kk=sum(k) l=f$V2 ll=sum(l) m=g$V2 mm=sum(m) aa<-read.table('out1.txt') # pro odtok bb<-read.table('out2.txt') cc<-read.table('out3.txt') dd<-read.table('out4.txt') ee<-read.table('out5.txt') ff<-read.table('out6.txt') gg<-read.table('out7.txt') hh2=aa$V2 hhh=sum(hh2)/170000 chh2=bb$V2 chhh=sum(chh2)/170000 ii2=cc$V2 iii=sum(ii2)/170000

- 91 -

jj2=dd$V2 jjj=sum(jj2)/170000 kk2=ee$V2 kkk=sum(kk2)/170000 ll2=ff$V2 lll=sum(ll2)/170000 mm2=gg$V2 mmm=sum(mm2)/170000 k1=hhh/hh k2=chhh/chh k3=iii/ii k4=jjj/jj k5=kkk/kk k6=lll/ll k7=mmm/mm ww=sum(k1,k2,k3,k4,k5,k6,k7) www=round(ww,7) write.table(www,file="sum_odtok_koef.txt",col.names=FALSE,row.names=FALSE)

- 92 -

Příloha č.7: Regresní analýza pro různé typy API Pro lineární regresní analýzu: API_30 120 Řada1 Lineární (Řada1)

100

API

80

60

40

20

y = -1.0173x + 82.836 2

R = 0.5549 0 0

20

40

60

80

V (poč.ztráty)

Pro mocninnou regresní analýzu: API_30 120 Řada1 Mocninný (Řada1) 100

API

80

60

40 -0.5334

y = 277.96x 2

R = 0.6743

20

0 0

20

40

60

V (poč.ztráty)

- 93 -

80

Pro exponenciální regresní analýzu:

API_30 120 Řada1 Exponenciální (Řada1) 100

API

80

60

40

y = 91.174e

20

-0.0209x

2

R = 0.6855 0 0

20

40

60

80

V (poč.ztráty)

Pro lineární regresní analýzu:

API_5 120 Řada1 Lineární (Řada1)

100

API

80

60

40 y = -0.6501x + 43.927 20

2

R = 0.2062

0 0

20

40

60

V (poč.ztráty)

- 94 -

80

Pro lineární regresní analýzu: API_5 + T_5 120 Řada1 Lineární (Řada1)

100

A PI

80 60 y = -0.6504x + 53.869

40

2

R = 0.2258 20 0 0

10

20

30

40

50

60

70

V (poč.ztráty)

Pro mocninnou regresní analýzu: API_5 120 Řada1 Mocninný (Řada1)

100

API

80 60 40 20 -0.6484

y = 166.2x 0

2

R = 0.2527 0

20

40 V (poč.ztráty)

- 95 -

60

80

Pro mocninnou regresní analýzu:

API_5 + T5 120

Řada1 100

Mocninný (Řada1)

API

80

60

40 -0.4371

y = 130.94x 20

2

R = 0.2359

0 0

10

20

30

40

50

60

70

V (poč.ztráty)

Pro exponenciální regresní analýzu:

API_5 120 Řada1 Exponenciální (Řada1) 100

API

80

60

40

20 y = 43.626e

0

2

0

20

40

60

V (poč.ztráty)

- 96 -

80

-0.0263x

R = 0.2739

Pro exponenciální regresní analýzu:

API_5 + T5 120

Řada1 Exponenciální (Řada1)

100

A PI

80 60 40 20

y = 51.863e 2

-0.0165x

R = 0.2231 0 0

10

20

30

40

50

V (poč.ztráty)

- 97 -

60

70

Příloha č.8: Průběh konvergence chyby:

Mocninný model APi(30)+T(30) epsilon -56.384 -31.347 -30.118 -29.793 -29.497 -28.871 -28.856 -28.791 -28.573 -28.349 -27.994 -27.935 -27.777 -27.746 -27.713 -27.657 -27.611 -27.595 -27.593 -27.447

počet iterací 148 841 3798 4071 4653 6216 7112 10765 11670 12790 13989 16632 17000 19512 20291 23075 23655 24660 25671 31012

Lineární model APi(30)+T(30) počet epsilon iterací -62.834 126 33.167 886 -30.774 1318 -30.535 2761 -29.771 2973 -29.606 3678 -29.584 4284 -29.544 4806 -29.356 5215 -28.644 5478 -28.093 6498 -27.581 7647 -27.119 12005 -27.114 27240 -27.087 39694 -27.029 54793

- 98 -

Mocninný model API (30) počet epsilon iterací -28.791 10765 -28.573 11670 -28.349 12790 -27.994 13989 -27.935 16632 -27.777 17000 -27.746 19512 -27.713 20291 -27.657 23075 -27.611 23655 -27.595 24660 -27.593 25671 -27.447 31012

Příloha č.9: Validace nevykazující dobrou shodu:

- 99 -

- 100 -

- 101 -

- 102 -

- 103 -