ČESKÁ ZEMĚDĚLSKÁ UNIVERZITA V PRAZE FAKULTA ŽIVOTNÍHO PROSTŘEDÍ Katedra vodního hospodářství a environmentálního modelování
DIPLOMOVÁ PRÁCE
Matematické modelování transportních procesů v porézním prostředí
Autor diplomové práce: Lucia Lóžiová Vedoucí diplomové práce: Ing. Michal Kuráž, PhD.
© Praha 2013
Čestné prohlášení Prohlašuji, že jsem diplomovou práci na téma „Matematické modelování transportních procesů v porézním prostředí“ vypracovala zcela sama, a že jsem všechny uvedené materiály a zdroje uvedla v seznamu literatury.
….……………………… V Praze dne 22. dubna. 2013
Lucia Lóžiová
Poděkování Ráda bych zde poděkovala všem, kteří mi byli při tvorbě této práce nápomocni. Především vedoucímu mé diplomové práce Michalu Kurážovi za poskytnuté rady a čas. Dále také rodině a přátelům za podporu.
Abstrakt Tato diplomová práce se zabývá modelováním podpovrchového odtoku na svahu malého horského povodí Modrava 2. Proudění porézním prostředím je popsáno deterministickým modelem s Richardsovou rovnicí. Pro experimentální povodí byl vytvořen konceptualizovaný 2D model příčného půdního profilu v programu Hydrus. V práci jsou popsány potřebné půdní charakteristiky, zákonitosti pohybu vody v porézním prostředí, hlavní rovnice proudění pro nasycenou i nenasycenou zónu a možnosti jejího řešení s počátečními a okrajovými podmínkami. Na základě těchto znalostí je sestaven konceptuální model, který simuluje odtok z povodí skrz definované hranice půdního profilu. Z modelovaných odtokových hodnot je posuzován srážkoodtokový proces na povodí.
Klíčová slova: proudění vody, konceptualizace, Richardsova rovnice, malé horské povodí, Modrava, Hydrus 2D/3D, srážkoodtokový proces, 2D modelování
Abstract This diploma work concerns modeling of subsurface outflow from a hill slope of a small mountainous catchment Modrava 2. Flow through a porous medium is described by Richards’ equation. Conceptualized 2D model of cross soil profile for experimental catchment was created in Hydrus program. Thesis describes required soil characteristics, physical laws for flow in porous media, main flow equation for saturated and unsaturated porous medium and possible ways to solve it with initial and boundary conditions., Conceptual model is built on these basis, which simulates outflow through defined soil profile boundaries. Rainfall-runoff process on the experimental catchment is valuated with these modeled outflow values.
Key words: water flow, conceptualization, Richards’ equation, small mountainous catchment, Modrava, Hydrus 2D/3D, Rainfall-runoff process, 2D modeling
OBSAH 1. ÚVOD ................................................................................................................. 10 2. CÍLE ................................................................................................................... 11 3. LITERÁRNÍ REŠERŠE ................................................................................... 12 3.1. VODA V PŮDĚ .................................................................................................... 13 3.2. PORÉZNÍ PROSTŘEDÍ ........................................................................................ 15 3.2.1. ZRNITOSTNÍ SLOŽENÍ....................................................................................... 15 3.2.2. PÓROVITOST.................................................................................................... 15 3.2.3. OBJEMOVÁ VLHKOST ...................................................................................... 16 3.2.4. RETENČNÍ ČÁRY PŮDNÍ VLHKOSTI ................................................................... 16 3.2.5. ANALYTICKÉ FUNKCE PRO VYJÁDŘENÍ RETENČNÍCH ČAR................................ 18 3.2.6. HOMOGENITA A HETEROGENITA, IZOTROPIE A ANIZOTROPIE........................... 20 3.2.7. HYDRAULICKÁ VODIVOST ............................................................................... 20 3.2.8. RETC A ROSETTA........................................................................................ 22 3.3. PROUDĚNÍ VODY V PORÉZNÍM PROSTŘEDÍ ...................................................... 24 3.3.1. ROVNICE KONTINUITY..................................................................................... 24 3.3.2. DARCYHO ZÁKON............................................................................................ 24 3.3.3. DARCY-BUCKINGHAMŮV ZÁKON .................................................................... 26 3.3.4. RICHARDSOVA ROVNICE ................................................................................. 26 3.4. ŘEŠENÍ ROVNICE PROUDĚNÍ............................................................................. 28 3.4.1. POČÁTEČNÍ PODMÍNKY .................................................................................... 28 3.4.2. OKRAJOVÉ PODMÍNKY..................................................................................... 28 3.5. KONCEPTUÁLNÍ MODEL ................................................................................... 30 3.5.1. KONCEPTUALIZACE PŘÍČNÉHO PROFILU POVODÍ ............................................. 30 3.5.2. KALIBRACE A VALIDACE ................................................................................. 31 4. CHARAKTERISTIKA ZÁJMOVÉHO ÚZEMÍ ........................................... 32 4.1. MODRAVSKÁ EXPERIMENTÁLNÍ POVODÍ ......................................................... 32 4.2. MALÁ MOKRŮVKA ........................................................................................... 33 4.2.1. POKRYV .......................................................................................................... 34 4.2.2. HYDROPEDOLOGICKÉ PODMÍNKY .................................................................... 35 4.2.3. SKLONITOST POVODÍ ....................................................................................... 36
5. METODIKA....................................................................................................... 38 5.1. POPIS DAT ......................................................................................................... 38 5.2. SESTAVENÍ MODELU ......................................................................................... 39 5.2.1. DOMÉNA A JEJÍ VÝPOČETNÍ SÍŤ........................................................................ 39 5.2.2. MATERIÁLOVÉ ROZDĚLENÍ.............................................................................. 40 5.2.3. POČÁTEČNÍ PODMÍNKA .................................................................................... 42 5.2.4. OKRAJOVÉ PODMÍNKY..................................................................................... 42 5.2.5. KALIBRACE A VALIDACE ................................................................................. 43 6. VÝSLEDKY ....................................................................................................... 44 7. DISKUZE ........................................................................................................... 47 8. ZÁVĚR ............................................................................................................... 48 9. POUŽITÁ LITERATURA ............................................................................... 49 10. PŘÍLOHY ......................................................................................................... 52
1.
ÚVOD
Tato diplomová práce se zabývá problematikou odtoku malého Šumavského povodí. Horské lesní oblasti podstatně ovlivňují srážkoodtokový proces svou infiltrační a retenční schopností. Forma a rychlost odezvy je při odtoku z horských povodí výrazně ovlivněna podpovrchovými půdními horizonty, informace o jejich schopnosti vést vodu jsou tedy pro predikci zásadní. Modelováním podpovrchového odtoku se zabývají např. Bayer a kol. (2004) na Šumavském povodí Liz, Pavelková a kol. (2012) na povodí Uhlířská v Jizerských horách či Tesař a kol. (2001) na dalších Šumavských povodích. Infiltrované vstupní srážky na malém horském povodí obvykle v profilu odtékají jen z poloviny rovnou do toku, část odteče do zásoby podzemní vody, odkud v sušším období dotuje základní odtok (Bayer a kol., 2004). Transformace srážky je dvoustupňová: v perkolační fázi srážka půdou pouze proteče rovnou do transportního kolektoru a v akumulační fázi se voda v půdě zdrží a do podloží neodtéká (Tesař a kol., 2001). Literární rešerše přehledně uvádí užívané půdní vlastnosti porézního prostředí. Dále jsou zde uvedeny základní fyzikální zákonitosti proudění v nasyceném i nenasyceném prostředí, používané pro deterministické modelování a způsob jejich řešení. V poslední části rešerše je popsána konceptualizace modelu jako způsob řešení proudění vody v porézním prostředí. V další kapitole je popsáno zájmové území z hlediska hydropedologických a geodetických charakteristik důležitých pro konceptualizaci povodí. Experimentální povodí bylo popsáno v několika publikacích, hydrologická stanice je zde vybudována již od roku 1998. V oblasti probíhají různé výzkumy Katedry vodního hospodářství
a
environmentálního
modelování,
zejména
hydrologického
hydropedologického charakteru. Povodí je tedy detailně prozkoumáno z hlediska vegetačního pokryvu, složení půdního profilu a kontinuálně se zde měří hodinové průtoky v litrech za sekundu, teploty vzduchu ve stupních Celsia a úhrny spadlých srážek v milimetrech s příslušným datem a hodinou měření. V metodice jsou použitá vstupní data popsána spolu s postupem sestavení modelu v programu Hydrus a zpracováním jeho výstupů.
-10-
2.
CÍLE
Cílem této práce je konceptualizace experimentálního povodí Malá Mokrůvka. K vypracování modelu je nejprve nezbytné shromáždit hydropedologické a geodetické charakteristiky zájmové oblasti, které mají na srážkoodtokové reakce povodí zásadní vliv. Na jejich základě se poté vytvoří model příčného profilu. Proces proudění vody v půdních horizontech je třeba popsat Richardsovou rovnicí pomocí deterministického modelu, do kterého vstupují měřené srážkové hodnoty. Cílem je upravit parametry modelu porovnáváním výstupných odtoků s měřenými pomocí kalibračních kritérií a posoudit srážkoodtokový proces.
-11-
3.
LITERÁRNÍ REŠERŠE
V literární rešerši je nejprve vysvětleno názvosloví a rozdělení podpovrchové vody a její důležitost v hydrologickém cyklu. Ve druhé podkapitole jsou uvedeny veškeré hydropedologické charakteristiky zásadní pro proudění vody. Vlastnosti půdy slouží pro modelování jako vstupní hodnoty, bez nichž se nelze obejít. V této části jsou také popsány důležité analytické funkce k predikci půdní vlhkosti a hydraulické vodivosti a možnosti určení jejich parametrů pomocí softwaru. Další část obsahuje základní rovnice, ze kterých se skládá řídící rovnice proudění v porézním prostředí, tedy rovnice kontinuity, Darcyho zákon a DarcyBuckinghamův zákon. Richardsova rovnice je zde uvedena ve dvou hlavních formách, v kapacitním a difúzním tvaru. V předposlední části o řešení rovnice proudění je vysvětleno nutné vymezení oblasti numerického řešení počátečními a okrajovými podmínkami s popisem jejich nejpoužívanějších typů. Poslední část rešerše popisuje princip konceptuálního modelu a jeho zjednodušujících zásad, které umožňují řešit podpovrchový odtok.
-12-
3.1. Voda v půdě Termín podpovrchová voda zahrnuje veškerou vodu nacházející se v půdě, tedy pod zemským povrchem. Prostředí, ve kterém voda proudí, se dělí na zónu nasycenou a nenasycenou, oddělené hladinou podzemní vody (Valentová, 2001). Nenasycená neboli vadózní zóna či zóna aerace je zespoda ohraničena hladinou podzemní vody nebo nepropustným podložím, horní hranici tvoří terén. V této zóně se pohybuje infiltrující a perkolující srážková voda, eventuálně vzlínající podzemní voda z nasycené zóny. Voda zaplňuje póry částečně nebo i úplně, zbytek je vyplněn vzduchem. Směr proudění je ve vadózní zóně převážně vertikální. (Císlerová & Vogel, 1998). Oblast nad hladinou podzemní vody se dále člení na kořenové, přechodné a kapilární pásmo (obrázek č. 3.01). Jednotlivá pásma mají proměnlivý rozsah i meze v závislosti na vlhkosti, hloubce hladiny podzemní vody a dosahu kořenů vegetace. V případě vysoké hladiny podzemní vody může kapilární voda dosáhnout kořenové zóny a přechodné pásmo tak neexistuje. Výška kapilárního pásma závisí na vzlínání, tedy na velikosti půdních pórů, těsně nad hladinou podzemní vody je prostředí téměř zcela nasycené (Valentová, 2001).
Obr. 3.01: Rozdělení vody ve vertikálním profilu (Valentová, 2001).
V nasycené zóně neboli zóně saturace se nachází voda podzemní, množství zadržené vody je zde podmíněno pouze objemem prostoru pórů (Vogel, 2004). Zóna je shora ohraničená hladinou podzemní vody nebo nepropustnou vrstvou. Horninové prostředí, ve kterém díky jeho propustnosti proudí voda je nazýváno zvodeň, kolektor či aquifer. Prostředí ve srovnání téměř nepropustné se označuje izolátor (Valentová, 2001). Kolektory se dále rozlišují na ty s napjatou hladinou shora i zdola -13-
ohraničené izolátorem a kolektory s volnou hladinou, které mají nepropustnou vrstvu pouze zespoda, nahoře se tvoří hladina podzemní vody (Císlerová & Vogel, 1998). Podzemní voda je součástí hydrologického cyklu zobrazeného na obrázku č. 3.02. Srážky ve formě deště nebo sněhu částečně odtečou jako povrchový odtok a část se infiltruje do půdy. Odtok je dále způsoben evapotranspirací skládající se z evaporace z povrchu půdy a vodní plochy a transpirace z rostlin (Valentová, 2001). Všechny tyto prvky jsou členy pro rovnici hydrologické bilance, která hydrologický cyklus vyjadřuje (Tourková, 1996). Půdní voda a její proudění je často rozhodujícím činitelem pro celý systém koloběhu vody v přírodě. Půda zadržuje cca desetinásobek vody v nádržích, má tedy v hydrologickém cyklu dominantní postavení (Kutílek, 1978). Celkově je toto množství odhadováno na 67 000 km3. Z hlediska sladké vody se jedná o velmi účinný zásobní systém (Lal & Shukla, 2004). Z důvodu velikosti této zásoby a pomalého proudění v horninovém prostředí jsou akumulační změny v nasyceném prostředí relativně malé, pomalé a značně komplikovanější v porovnání s rychlou povrchovou odezvou (Bear & Cheng, 2010).
Obr. 3.02: Schéma hydrologického cyklu (Bear & Cheng, 2010).
-14-
3.2.
Porézní prostředí
Pórovitý materiál zahrnuje půdní částice (pevnou fázi) prostoupené složitým systémem vzájemně propojených pórů. Tyto prostory jsou zaplněny kapalnou a plynnou fází (Kutílek, 1984), jejich geometrie je velice obtížně popsatelná, je zapotřebí zjednodušujících úmluv. K popisu porézního prostředí vzhledem k proudění vody jsou definovány hydropedologické charakteristiky, jako zrnitostní složení, pórovitost, permeabilita a z ní vycházející hydraulická vodivost. Tyto vlastnosti slouží pro modelování jako nepostradatelné vstupní údaje a na jejich správném určení značně závisí výsledky simulací (Císlerová & Vogel, 1998). Reálné prostředí je prostorově i časově vysoce variabilní, což popis půdních charakteristik komplikuje. Před samotným vytvořením modelu je nutností podrobná analýza půd v zájmovém území (Císlerová, 1989).
3.2.1.
Zrnitostní složení
Zrnitost určuje velikost poloměru zrn a jejich podíl ve zkoumaném vzorku půdy. Popisuje se zrnitostní křivkou, ta je při výzkumu základním ukazatelem, jelikož je možné jí jednoduše zjistit z porušeného vzorku. Zrnitostní složení je pro tvar a velikost pórů jedním z nejpodstatnějších faktorů (Císlerová & Vogel, 1998). Podle zastoupení frakcí jemnozemě (< 2mm) je prováděno hodnocení zrnitosti. Půdní částečky jemnozemě se dle velikosti dělí na jíl (< 0.002mm), prach (0.0020.05mm) a písek (0.05-2mm). Pomocí trojúhelníkového klasifikátoru se půdy rozdělují na prach, prachovou hlínu, hlínu, jílovitou hlínu, prachový jíl, jíl, písčitý jíl, písčitou jílovitou hlínu, písčitou hlínu, hlinitý písek a písek (Němeček, 2001).
3.2.2.
Pórovitost
fází v celkovém objemu. Je formulována vzorcem 3.01 (Kutílek, 1984).
Pórovitost vyjadřuje podíl pórů, tedy prostorů zaplněných plynnou a kapalnou
kde:
je pórovitost
je objem pórů je objem neporušeného vzorku půdy -15-
3.01
Efektivní pórovitost Pro proudění vody a transport látek mají význam pouze propojené póry, které se na proudění podílejí. Pro tento účel je zavedena efektivní pórovitost, která zahrnuje pouze objem pórů účastnících se proudění (Mls, 1984). Se stoupající velikostí zrna pórovitost obvykle klesá a efektivní pórovitost stoupá (Tourková, 1996). To, že voda v pórech neproudí, bývá někdy způsobeno také zaslepeným koncem, kde je voda prakticky nehybná, nebo velmi malým průřezem póru, např. při zrnech velmi malé velikosti zamíchaných mezi velmi velké (Bear & Cheng, 2010).
3.2.3.
Objemová vlhkost
objemovou. Objemové vyjádření, uvedené v rovnici 3.02, je pro bilancování půdní Množství vody v půdě, tedy vlhkost lze definovat jako hmotnostní, nebo
vody a různé nasycenosti prostředí vhodnější (Kutílek, 1984).
kde:
3.2.4.
3.02
je objemová vlhkost je objem vody
je objem neporušeného vzorku půdy
Retenční čáry půdní vlhkosti
Retenční křivka je grafickým znázorněním vztahu objemové půdní vlhkosti a sacím tlakem neboli tlakovou výškou. Zobrazuje retenční schopnost materiálu, tedy schopnost zadržovat vodu proti působení vnější síly (Kutílek, 1984). Její tvar závisí především na geometrických vlastnostech pórového systému. Body retenční čáry lze zjistit měřením vlhkostí odpovídajících tlakovým výškám po krocích, její určení je však složité a nejednoznačné. Průběh čáry je pro každou půdu odlišný, je třeba jí vždy stanovit zvlášť (Kutílek & kol., 2004).
-16-
Hystereze retenční křivky Při plnění a prázdnění pórů vodou má retenční křivka odlišný tvar, jak je vidět na obrázku 3.03 projevuje se zde hystereze. Tento jev je způsoben především vzduchem uzavíraným ve „slepých“ pórech a rozdílnou hodnotou smáčecího úhlu na suchém a smáčeném povrchu (Kutílek, 1984). U hrubozrnných materiálů může být hystereze významná, její matematický popis je značně složitý. Ve většině případů je hystereze zanedbána (Císlerová & Vogel, 1998).
Obr. 3.03: Hystereze retenční křivky, rozdíl mezi drenážní a zvlhčovací větví (Lal & Shukla, 2004).
Vodní kapacita vodní kapacitu, která popisuje sklon retenční křivky, její vyjádření 3.03 uvádí
Z retenční křivky lze odvodit další charakteristickou hydraulickou vlastnost,
např. Kutílek (1984), Císlerová & Vogel (1998).
kde:
je vodní retenční kapacita
3.03
je objemová vlhkost tlaková výška
Specifická storativita Retenční vodní kapacita má v nasyceném prostředí nulovou hodnotu, v tomto případě se na jejím místě uvádí specifická storativita. Storativita je vlastnost prostředí přijímat a uvolňovat určité množství vody. Specifickou storativitou se označuje objem vody, který zvodnělá vrstva uvolní při jednotkovém poklesu piezometrické výšky (Pech, 2010). -17-
3.2.5.
Analytické funkce pro vyjádření retenčních čar
Pro popis retenčních čar existují analytické funkce, kterými je možné přibližně používá efektivní vlhkost vyjádřena rovnicí 3.04 vyjadřující stupeň nasycení půdy
vystihnout hledanou závislost vlhkosti a tlakové výšky. Pro výpočet se obvykle
(např. Císlerová, 1989, Kodešová, 2012).
kde:
3.04
je efektivní vlhkost
je nasycená (saturovaná) vlhkost je reziduální vlhkost
Rovnice Brookse & Coreyho (1964), retenční křivka této funkce je zobrazena na obrázku 3.04. Rovnice 3.05, Jednodušší ze dvou nejčastěji uváděných výrazů navrhli Brooks & Corey
pojmenovaná po svých autorech obsahuje dva parametry: vyjadřuje rozdělení Druhý parametr je , původně Brooksem & Coreym (1964) definovaný jako
velikosti pórů, podle Císlerové (1989) úzce souvisí se sklonem zrnitostní křivky.
probublávací tlak (bubbling pressure), což je minimální kapilární tlak, při kterém začíná probublávat vzduch. Ve skutečnosti se však jedná o vstupní hodnotu vzduchu, která je obvykle nižší (Císlerová, 1989).
kde:
|| 1
!
je efektivní vlhkost je tlaková výška
, jsou parametry rovnice
-18-
"
#
3.05$
Rovnice van Genuchtena Tento výraz je vyjádřen v rovnici 3.06 pomocí tří parametrů ', ( a ). Parametr (
Druhou často zmiňovanou aproximační funkci uvedl van Genuchten (1980).
je obvykle vypočten z parametru ): ( 1 , ) , 1. Retenční křivka funkce van *
Genuchtena je znázorněna na obrázku 3.04.
kde:
+
1 -1 . '||+ / 1 je efektivní vlhkost
', ( a )
"0
#0
$
3.06
je tlaková výška jsou parametry rovnice
Obrázek 3.04: Retenční křivka dle van Genuchtena (vlevo) a dle Brookse & Coreyho (vpravo) (RETC)
Rozdíl mezi uvedenými funkcemi se nejvíce projevuje v blízkosti nasycené vlhkosti. Rovnice Brookse & Coreyho je méně flexibilní a nevhodná pro materiály se širokým zrnitostním spektrem, jejichž retenční čáry se esovitě prohýbají. Van Genuchtenova rovnice je vzhledem k většímu počtu parametrů flexibilnější, esovitý tvar tedy vystihne daleko lépe. Nevýhodou je však nulová vstupní hodnota vzduchu, která je pro většinu půd nereálná (Císlerová, 1989).
-19-
3.2.6.
Homogenita a heterogenita, izotropie a anizotropie
V homogenním prostředí se vlastnost se změnou pozice nemění. Pokud je na poloze uvnitř závislá, jedná se o prostředí heterogenní. Vlastnosti se mohou měnit skokově, nebo postupně (Pech, 2010). Prostředí může být izotropní nebo anizotropní vždy vzhledem k nějaké vlastnosti, např. permeabilitě a z ní vycházející hydraulické vodivosti. Pokud se vlastnost v různých směrech nemění, je prostředí izotropní, jestliže je na směru závislá tak je anizotropní. V případě vrstevnatých materiálů je podél vrstev obvykle větší vodivost. Daná charakteristika se v anizotropním prostředí vyjadřuje tenzorem stupně dle dimenze řešeného modelu (Mls, 1984).
3.2.7.
Hydraulická vodivost
Propustnost
Propustnost neboli permeabilita 0 charakterizuje schopnost porézního prostředí
hydraulická vodivost 1 3.07, kterou je zapotřebí stanovit pro každou kapalinu propouštět vodu bez ohledu na vlastnosti kapaliny. Z propustnosti vychází
zvlášť (Pech, 2010).
kde:
1
034 5 6
1
0
je hydraulická vodivost
6
je gravitační zrychlení
34 5
3.07
je permeabilita je hustota kapaliny
je dynamická viskozita
Nasycená hydraulická vodivost
Koeficient hydraulické vodivosti 17 vyskytující se v Darcyho rovnici závisí jak
na vlastnostech prostředí, tak na vlastnostech proudící kapaliny (Kutílek & kol., 2004). Nasycenou hydraulickou vodivostí je míněna schopnost prostředí vést vodu ve stavu, kdy jsou všechny póry vyplněné vodou (Lal & Shukla, 2004).
-20-
V izotropním prostředí je 1 skalár, pro anizotropní tenzor. Pro trojrozměrné
prostředí lze hydraulickou vodivost zapsat maticí 3.08 (Bear & Cheng, 2010). 1:: 1 91;: 1<:
1:; 1;; 1<;
1:< 1;< = 1<<
3.08
Nenasycená hydraulická vodivost V nenasyceném prostředí proudí voda jen v některých, menších pórech. Křivku hydraulické vodivosti je tedy možné předpovědět na základě teorie kapilárního modelu z retenční křivky půdní vlhkosti. Tato teorie aproximuje pórový systém na svazek rovnoběžných kapilár kruhového průřezu, ve kterých probíhá pouze laminární proudění. Hlavním předpokladem je, že distribuční funkci rozdělení poloměru pórů je možné zjistit z retenční křivky (Kodešová, 2012). nejčastěji používají vztahy odvozené Burdinem 3.09 a Mualemem 3.10, ze Pro určení hodnot hydraulické vodivosti různě nasyceného prostředí se
kterých lze v kombinaci s tradičními funkcemi pro retenční čáry odvodit analytické vyjádření relativní hydraulické vodivosti.
? 1 ? * @C ? @C
AB
? ? 1 C.D E F * @C ?
@C
AB
kde:
1 je relativní hydraulická vodivost
3.09 3.10
je efektivní vlhkost je tlaková výška
3.05 a 3.06 s výrazy 3.09 a 3.10. Pokud vycházíme z rovnice Brookse & Předpověď hydraulické vodivost se provádí v různých kombinacích rovnic
Coreyho, vznikne výraz 3.11, pro van Genuchtena je to rovnice 3.12. Do funkcí
se dosazují hodnoty parametrů G, H a I podle vybraného Mualemova nebo Burdinova modelu (Císlerová, 1989).
-21-
kde:
1 1 J!K
1 je relativní hydraulická vodivost 1
3.11
je nasycená hydraulická vodivost
G 2, H 3 pro model Burdina
je parametr rovnice Brookse & Coreyho
G 2, H 2.5 pro model Mualema
kde:
(, )
J ⁄/ /
1 1 L1 M1 *
O P
jsou parametry rovnice van Genuchtena, ( 1 I ⁄)
G 1, H 2, I 2 pro model Burdina
3.12
G 2, H 0.5, I 1 pro model Mualema
Křivka nenasycené hydraulické vodivosti je stejně jako retenční křivka zatížena hysterezí, její průběh je v reálném prostředí odlišný pro odvodňování a zvlhčování. Laboratorní stanovení hydraulické vodivosti nenasycené půdy je obtížné, používá se zejména její nepřímé určení (Kodešová, 2012).
3.2.8.
RETC a ROSETTA
Pro komplexní vyšetření parametrů hydraulických funkcí je možné použít program RETC, zdarma dostupný na http://www.pc-progress.com. Rosetta je implementován v programu RETC i HYDRUS. RETC používá modely Brookse & Coreyho a van Genuchtena pro stanovení retenčních křivek a kapilární modely Burdina a Mualema k předpovědi průběhu nenasycené hydraulické vodivosti, tyto modely lze v RETC různě kombinovat. Program je možné využít pro nalezení nenasycených hydraulických funkcí z již zjištěných parametrů retenční křivky, k predikci nenasycené hydraulické vodivosti a difuzivity ze změřených bodů retenční křivky či naopak k predikci retenční čáry z hydraulické vodivosti a difuzivity (van Genuchten & kol., 1991). Parametry hydraulických modelů je také možné předpovědět pedotransferovými funkcemi implementovanými v programu Rosetta. Hydraulické vlastnosti jsou pro modelování půdního prostředí nepostradatelné, jejich měření je však někdy z časových či finančních důvodů nemožné. -22-
Na výběr je 5 modelů umožňujících odhad parametrů retenční čáry a nasycené i nenasycené hydraulické vodivosti. Predikci lze realizovat z různého množství informací. Nejjednodušší variantou je odhad ze zrnitostní třídy, Rosetta pak pouze vygeneruje průměrné hodnoty parametrů pro danou třídu. Složitější modely jsou založeny na neuronových sítích, vstupním údajem je procentuální zastoupení zrnitostních frakcí, které lze rozšířit o objemovou hmotnost, polní kapacitu objemová vlhkost při tlakové výšce 3,3 m a bod vadnutí - objemová vlhkost při tlakové výšce 150 m (Schaap & kol., 2001).
-23-
3.3.
Proudění vody v porézním prostředí
Proudění vody v půdě, tedy v porézním prostředí lze popsat spojením dvou fundamentálních rovnic: rovnicí kontinuity, Darcyho zákonem v případě nasyceného prostředí a Darcy-Buckinghamovým zákonem pro nenasycené prostředí (Kutílek, 1978, Císlerová, 1989).
3.3.1.
Rovnice kontinuity
Formulací zákona zachování hmoty je rovnice kontinuity vyjadřující bilanci s konstantní hustotou je rovnice kontinuity vyjádřena ve tvaru 3.13. (Císlerová & objemové vlhkosti porézního prostředí v čase. Pro nestlačitelnou kapalinu
Vogel, 1998).
Q. RS
3.3.2.
T TU
3.13
Darcyho zákon
výšce # 0 proudí voda ve všech pórech a nasycená objemová vlhkost se rovná Pro zjednodušení řešení nasyceného proudění je předpokládáno, že při tlakové
pórovitosti . Ve skutečnosti je však v některých pórech uzavřený vzduch (Kutílek,
v nasyceném prostředí. Rovnice 3.14 byla odvozena na základě experimentálního
2004). V roce 1956 publikoval Henry Darcy vztah pro jednorozměrné proudění
sledování průtoku ve válci naplněném pískem a později rozšířena i na ostatní
propustné materiály (Hálek & Švec, 1973). Tento zákon představuje lineární závislost rychlosti proudění na hydraulickém gradientu, který charakterizuje změnu energie pohybující se vody (Kazda, 1983).
kde:
R
1 V
R 1V
je objemový průtok válcem je hydraulická vodivost
* h? W
* ? W
je plocha průtočného průřezu je rozdíl piezometrických výsek
je vzdálenost profilů -24-
3.14
Obecně pak pro osu Y kladnou směrem vzhůru s počátkem na spodní hladině
vody platí rovnice pro trojrozměrné proudění 3.15 (Kutílek & kol., 2004). kde:
R
1
RS 1Z
je objemový průtok
3.15
je hydraulická vodivost je hydraulická výška
Hydraulická výška je definována jako součet tlakové výšky a vektoru
0 0 geodetické výšky YS [0\, pro osu Y kladnou směrem vzhůru pak platí ZY [0\ Y 1 rozepsat do jednotlivých složek 3.16, (Říha, 1997).
(Kuráž, 2011). S použitím tenzorového vyjádření hydraulické vodivosti lze průtok T T T 1:; 1:< . 1 T] T^ TY T T T 1;; 1;< . 1 R; 1;: TY T] T^ T T T R< 1<: 1<; 1<< . 1 T] T^ TY R: 1::
3.16
Meze platnosti Darcyho zákona lineární. Taková závislost však platí jen v určitém rozsahu Reynoldsova čísla _`
Darcyho zákon vyjadřuje závislost rychlosti na hydraulickém gradientu jako
vyjádřeného výrazem 3.17 (Mls, 1984). S rostoucí rychlostí dochází k pozvolnému
postlineární. Obvykle se jeho platnost uvažuje do _` 1 až 20 (Kazda, 1983). Tato
odklonu od Darcyho zákona vlivem setrvačných sil, oblast tohoto proudění se nazývá
hodnota je často mylně považována za rozhraní laminárního a turbulentního proudění vody v půdě. Laminární proudění se však může zachovat při rychlostech daleko
vyšších, než linearita Darcyho zákona. Existuje i spodní hranice platnosti, kdy k odklonu od linearity dochází vlivem zvýšených molekulárních sil, oblast se nazývá se prelineární. Tento jev je obvykle úplně zanedbáván, odchylka je totiž velmi malá (Hálek & Švec, 1973).
-25-
kde:
a
b
3.3.3.
_` je rychlost proudění
a b
3.17
je charakteristický rozměr zrna je kinematická viskozita
Darcy-Buckinghamův zákon
Edgar Buckingham odvodil rovnici 3.18 pro proudění v nenasyceném
prostředí později, nezávisle na Darcyho zákoně. Je podobná rovnici pro nasycené prostředí, do výpočtu je však zahrnuta funkce nenasycené hydraulické vodivosti v závislosti na objemové vlhkosti (Císlerová & Vogel, 1998).
kde:
R
1
3.3.4.
RS 1Z
je objemový průtok
3.18
je hydraulická vodivost je objemová vlhkost je hydraulická výška
Richardsova rovnice
Richardsova rovnice je použitelná za zjednodušujících předpokladů: geometrie pórového prostředí může být ztvárněna kapilárním modelem nehybných a nedeformovatelných pórů, tlak vzduchu v pórech je atmosférický a neměnný, proudění vzduchu nemá žádný vliv na proudění vody a proudící voda je nestlačitelná (Císlerová & Vogel, 1998).
Řídící rovnice 3.19 pro řešení nestacionárního proudění formuluje základní
vztah mezi objemovou vlhkostí a průtokem v čase a prostoru spojením Darcy-
Buckinghamova nebo Darcyho zákona a rovnice kontinuity (Kutílek & kol., 2004). Z. M1Zh . ZzO
-26-
T TU
3.19
Kapacitní tvar
Jednou z forem Richardsovy rovnice je kapacitní tvar 3.20 vyjadřující
hodnoty v závislosti na tlakové výšce . Tento tvar také zahrnuje retenční vodní
kapacitu C definovanou v rovnici 3.03. (Císlerová & Vogel, 1998). Q. 1Q .
T1<< T TY TU
3.20
T1<< T V7 TY TU
3.21
V nasyceném prostředí se vodní kapacita nahradí specifickou storativitou V7 ,
vznikne tak rovnice pro proudění podzemní vody 3.21 (Kuráž, 2011). Q. 1Q .
Difúzní tvar
formu Richardsovy rovnice 3.22. V tomto tvaru se objevuje difuzivita půdní vody,
Pokud je neznámou funkcí v rovnici objemová vlhkost, jedná se o difúzní
kterou lze vyjádřit výrazem 3.23 (Kutílek & kol., 2004). Q. dQ .
d 1
T1<< T TY TU 1
-27-
3.22
3.23
3.4.
Řešení rovnice proudění
Matematické modelování vyžaduje mnoho zjednodušení, vytvořený model se materiálových charakteristik, viz kapitola 3.5. (Vogel, 2004). Kromě charakteristik pak nazývá konceptuální. Aproximace se týkají např. režimu proudění čí
popsaných v předchozích kapitolách je řešení problému nezbytné vymezit
počátečními a okrajovými podmínky (Kazda, 1983).
3.4.1.
Počáteční podmínky
3.24 předepisuje hodnotu proměnné pro všechny body modelované oblasti Ω
Počáteční podmínkou je určena výchozí hodnota neznámé funkce. Podmínka
v počátečním čase UC , nebo je možné určit počáteční vlhkostní pole podmínkou 3.25 (Bear & Cheng, 2010).
e, UC C e fg he i Ω e, UC C e fg he i Ω
3.4.2.
3.24 3.25
Okrajové podmínky
Na všech hranicích oblasti je nutné určit okrajové podmínky, je třeba definovat průběh hodnot pro celou dobu trvání simulace. U rovnice proudění vody se používají dva typy okrajových podmínek, a to Dirichletova a Neumannova (Císlerová, 1989). Dirichletova okrajová podmínka Tato podmínka je často nazývána také stabilní, z fyzikálního hlediska se jedná určité části hranice řešené oblasti proudění Ω. V případě známé tlakové výšky Г na o tlakovou okrajovou podmínku. U tohoto typu je známá hodnota neznámé funkce na
části hranice Г se podmínka definuje předpisem 3.26, podmínku lze definovat také
pomocí objemové vlhkosti jako výraz 3.27 (Císlerová & Vogel, 1998). e, U Г fg he, U i Г j k$0, l$ e, U Г fg he, U i Г j k$0, l$
3.26 3.27
Tlaková okrajová podmínka se používá pro charakterizaci výtopové infiltrace
s konstantní hladinou na horním okrají půdního profilu, k určení hladiny podzemní vody uvnitř půdního profilu (Kodešová, 2012), nebo může znamenat hranici za kterou je řeka či nádrž (Bear, 1972). -28-
Neumannova okrajová podmínka Nestabilní neboli toková okrajová podmínka je použita, známe-li rychlost
proudění v kolmém směru přes hranici oblasti Гm . Voda může do systému přitékat,
nebo z něj odtékat. Neumannova okrajová podmínka se definuje předpisem 3.28
vztaženým k vnější normále hranice )nS, tok musí být vyjádřen pomocí rovnice proudění s proměnnou h nebo (Císlerová & Vogel, 1998). 1 o
Te, U . )p q q Г fg he, U i Г j k$0, l$ T)nS
3.28
Neumannova podmínka může vyjadřovat srážky nebo výpar na okraji oblasti vyjadřujícím terén, nebo měřený odtok ve dně (Kodešová, 2012). podobě jako nulový tok „no flow“ 3.28, kdy je výraz roven nule. Tato podmínka
V inženýrské praxi se toková okrajová podmínka používá také v homogenní
může reprezentovat například nepropustné geologické podloží (Mucha & Šestakov, 1987).
Te, U 1 o . )p q 0 fg he, U i Г j k$0, l$ T)nS
3.28
Další používanou variantou je také volná drenáž „free drainage“ 3.29, kde je
derivace tlakové výšky vůči hranici nulová, voda tak odtéká pouze vlivem geodetického gradientu. Tok touto hranicí je tedy roven hodnotě nenasycené hydraulické vodivosti (Kuráž, 2011).
T e, U 0 fg he, U i Г j k$0, l$ T)nS
3.29
Dále se používá podmínka výronová plocha „seepage face“, ta je vhodná například pro výtok podzemní vody do nádrže nebo toku. Tato podmínka je nulovým tokem až do chvíle kdy dojde k nasycení a voda pak touto hranicí vysakuje volně na povrch (Bear, 1972). Fyzikálně se tedy v nenasyceném stavu chová jako Neumannova okrajová podmínka „no flow“ a při nasycení se změní na Dirichletovu matematické vyjádření viz výrazy 3.30 (Kuráž, 2011).
ve formě nulové tlakové výšky, v průběhu simulace se mezi nimi „přepíná“, pro T e, U . )p e 0 fg he, U " 0 fg he, U i Г j T)nS T e, U . )p e " 0 fg he, U i Г j e, U 0 fg h T)nS -29-
k$0, l$
k$0, l$
3.30
3.5.
Konceptuální model
Jedním z nejefektivnějších způsobů ve výzkumu komplexního celku je jeho konceptualizace. Skutečný komplikovaný systém je fakticky neřešitelný, podle tohoto přístupu ho lze zastoupit zjednodušeným fyzikálním jevem, který je v rámci dosažitelných nástrojů řešitelný (Bear, 1972). Tato metodologie snižuje složitost simulačního modelu, zatímco se zachová validita výsledků. Je ale nutné respektovat způsob použití modelu a jeho náležitosti (Pachepsky a kol., 2006). Konceptualizace spočívá ve schematizaci řešené oblast zjednodušení a podmínek proudění. Je třeba vymezit režim proudění na ustálené či neustálené, laminární, formulovat konstitutivní vztahy definující porézní materiál, určit dimenzi modelu, zahrnout předpoklad homogenity či anizotropie, vybrat časové intervaly pro neustálené proudění a vhodný matematický popis jevu (Mucha & Šestakov, 1987).
3.5.1.
Konceptualizace příčného profilu povodí
Vyřešit exaktní fyzikální model horského lesního povodí je prakticky nemožné, vzhledem k mnohým nepravidelnostem a heterogenitě půdního prostředí. Realizuje se tedy model značně zjednodušený s přibližným řešením, který lze použít k objemové predikci. Možnou aproximací je konceptuální model s reálným podpovrchovým tokem. Vstupující srážka je v půdním profilu rozdělena na gravitační vertikální odtok a odtok po svahu. Povrchový odtok je obvykle realizován ve vrchní velice propustné vrstvě profilu (Vogel, 2004). Řídící rovnice proudění je analyticky řešitelná jen u velmi jednoduchých úloh. Pro složitější problémy je nezbytné zvolit numerické řešení formou simulačního programu. Vymezující součástí modelu jsou předem zjištěné a dle konceptuálního přístupu
zjednodušené
charakteristiky.
Výstupem
jsou
pak
informace
o
časoprostorovém vývoji zkoumaných hodnot v modelované oblasti a zadaném období s informacemi a celkové bilanci a objemovém toku přes hranice oblasti (Císlerová & Vogel, 1998).
-30-
3.5.2.
Kalibrace a validace
Přesnost modelu a jeho vypovídající schopnost je třeba hodnotit také podle míry schematizace problému s vlivem přijatých předpokladů a zjednodušení. Důležitou roli hraje kalibrace modelu, která sestává z upravování a změn určitých parametrů až d dosažení přijatelné shody mezi měřenými a simulovanými hodnotami. Metody optimalizace mohou být dle počtu upravovaných parametrů statistické, ale i založené na intuici řešitele s metodou pokus omyl. Validací modelu je myšleno ověření parametrů získaných kalibrací (Říha, 1997).
-31-
4. 4.1.
CHARAKTERISTIKA ZÁJMOVÉHO ÚZEMÍ Modravská experimentální povodí
Fakulta lesnická a environmentální ČZU vybudovala po kůrovcové kalamitě v letech 1994–1995 z prostředků grantu MŽP ČR 14-053/97 v roce 1998 tři hydrologické stanice v oblasti Modravy a Kvildy. Povodí byla vybrána tak, aby měly co nejpodobnější základní fyziogeografické vlastnosti, liší se pouze porostem. Jednalo se o odumřelý les, zdravý les a pokácený les. Pozorování na modravských povodích začalo z důvodu hledání optimálního způsobu hospodaření se smrkovými porosty zasaženými kůrovcovou kalamitou, tedy jestli je vhodnější postiženou oblast vytěžit, nebo ponechat bez zásahu. Cílem výzkumu bylo zejména srovnání retenční schopnosti při povodních a hydrologické bilance na povodích s odlišnou kulturou, i vyjádření významnosti zásahů člověka. Výzkum se původně zaměřoval pouze na měření teploty, srážek a průtoků v letních měsících z důvodu technické nedostupnosti elektrické rozvodné sítě v těchto odlehlých lokalitách (Kuna & Kuřík, 1998). Poloha všech povodí je vyznačena na obrázku 4.01 V současné době jsou povodí spravována Katedrou vodního hospodářství a environmentálního modelování Fakulty životního prostředí (KVHEM, 2011).
Obr. 4.01: Mapa s vyznačením polohy modravských povodí (KVHEM, 2013).
-32-
4.2.
Malá Mokrůvka
Povodí popisuje například Pavlásek & kol. (2006, 2009), Jačka (2009), (Jačka & Pavlásek, 2010), Jačka & kol. (2012) a mnoho studentů Fakulty životního prostředí s jím zabývá ve svých bakalářských a diplomových pracích (např. Bašta, 2008, Jindrová, 2011, Petrželková, 2011). Území je tedy podrobně prozkoumané a k jeho další analýze a modelování je dostupných mnoho použitelných dat. Experimentální povodí Modrava 2 neboli „Malá Mokrůvka“, kterým se práce zabývá, je součástí pramenné oblasti Ptačího potoka, je hydrologické pořadí je 1- 0801-002. Nachází se 5km jižně od Filipovy Huti při hranici s Bavorskem, místní název lokality je „Medvědí doupě“. Plocha povodí je 0,16km2 (Pavlásek & kol., 2006). Rozkládá na severovýchodním svahu Malé Mokrůvky a severozápadním svahu Mrtvého vrchu. Oblast je vyznačena na obrázku 4.02 Zeměpisné souřadnice závěrného profilu jsou 48°58’27,3’’ severní šířky a 13°30’42,5’’ východní délky (Bašta, 2008). Nejvyšším místem povodí je vrchol Malé Mokrůvky, 1 330 m n. m.. Nejníže položeným bodem je uzávěrový profil, 1 188 m n. m.. (Jačka & kol., 2012).
Obr. 4.02: Mapa s vyznačením povodí Malé Mokrůvky (http://www.mapy.cz/)
Na povodí se měří průtok, srážky a teplota. V uzávěrových profilech povodí jsou vybudovány měrné přelivy. Průtok je počítáván z přepadové výšky na trojúhelníkovém průřezu ostrohranné přelivné stěny typu Thomson, která je automaticky snímána tlakovým čidlem. Měrný přeliv je vidět na obrázku 4.03. Srážky se měří v blízkosti přelivu překlopným srážkoměrem se záchytnou plochou -33-
200 cm2 ve výšce 1 m nad zemským povrchem. Teplota je měřená ve výšce 2 m nad zemským povrchem a konduktivita vody na měrném přelivu (Pavlásek & kol., 2009). Od roku na povodí Modrava 2 probíhá také monitoring základních charakteristik sněhové pokrývky (výška sněhu a objemová hmotnost sněhu), který byl v roce 2007 rozšířen i o monitoring kvalitativních charakteristik sněhové pokrývky, povodí je vybaveno devatenácti sněhoměrnými latěmi (KVHEM, 2013).
Obr. 4.03: Fotografie měrného přelivu Modravy 2 (KVHEM, 2011)
4.2.1.
Pokryv
V lokalitě byl vytěžen původní smrkový porost zasažený kůrovcem, starý až 160 let. Po těžbě zde zůstaly tlející větve, pařezy, půda udusaná těžebním zařízením (Pavlásek & kol., 2006). Zůstatky po těžbě mohou ovlivňovat srážkoodtokový proces a hydrologický režim povodí vytvářením preferenčních tras proudění vody po povrchu (Jačka & kol., 2012). V současné době je oblast pokryta mladými 10-15tiletými dřevinami. Jedná se zejména o smrk ztepilý (Picea abies). Dále se zde vyskytuje javor klen (Acer pseudoplatanus), jeřáb ptačí (Sorbus aucuparia) (Jačka & kol., 2012).
-34-
Bylinné patro vegetačního pokryvu tvoří zejména třtina chloupkatá (Calamagrostis villosa), a brusnice borůvka (Vaccinium myrtillus), dále papratka horská (Athyrium distentifolium), bika lesní (Luzula sylvatica), metlička křivolaká (Avenella flexuosa) a ploník obecný (Polytrichum commune) (Jindrová, 2011). Na obrázku 4.04 je fotografie pokryvu na povodí Modrava 2.
Obr. 4.04: Fotografie pokryvu Modravy 2 (KVHEM, 2011)
4.2.2.
Hydropedologické podmínky
Převládá zde pudní typ podzol, subtypem je podzol modální. Organický horizont nadložního humus O má hloubku 5 cm, překrývá se místy se slabě vyvinutým humusovým horizontem Ah o hloubce 1 až 5 cm. Následuje vybělený eluviální vrstva horizontu E v rozsahu 8 až 15 cm, iluviální černorezivý humusoseskvioxidický Bhs o tloušťce 15cm a rezivý seskvioxidický horizont Bs měřící 25cm (Jačka & kol., 2012). Pórovitost je až 86% v případě profilu O, E má 51%, Bhs 38%. Půdy jsou hlinitopísčité a písčitohlinité Půdní profil má celkově hloubku až 80 cm (Jačka & Pavlásek, 2010).
-35-
Rozbor půd v zájmové oblasti, zpracovaný Petrželkovou (2011) ukazuje na celkově malý obsah jílu ve všech horizontech (max. 6%). Veškeré zastoupené půdy lze označit jako hlinité písky až písky s menším obsahem skeletu. Zrnitostní složení jednotlivých horizontů v procentech je uvedeno v tabulce 4.05. Tab. 4.05: Zrnitostní složení jednotlivých půdních horizontů
E Bhs Bs C
skelet 6,11 17,73 18,85 2,91
jíl 1,9 2,1 1 6
jemnozem prach písek 22,1 76 10,9 87 15 84 17 77
Nasycená hydraulická vodivost Ks určená Jačkou & Pavláskem (2010) Guelphským permeametrem byla průměrně 3,1 x 10-5 ms-1. Pro jednotlivé horizonty je to hodnota 3,3 x 10-5 ms-1 společně pro O, Ah a E, 1,6 x 10-5 ms-1 pro horizont Bhs a 2,0 x 10-6 ms-1 pro Bs. Spodní horizont s pískovými sedimenty C má tuto hodnotu 1,2 x 10-5 ms-1, tedy číslo vyšší než hodnota předchozí. Dle Jačky & kol. (2012) při měření bez ohledu na vegetační pokryv jednoválcovou infiltrační metodou má Ks průměrnou hodnotu 4,211 x 10-5 ms-1, a Guelphským permeametrem 3,275 x 10-5 ms1
, tedy hodnota lišící se jen nevýznamně od výsledků publikovaných v roce 2010.
4.2.3.
Sklonitost povodí
Hodnoty průměrné sklonitosti povodí napravo i nalevo od údolnice jsou zjištěné Baštou (2008). Pravá mírnější část má sklon 12,5%, je také méně členitá než levobřežní strana. Oblast vlevo od údolnice je průměrně sklonitá 22,6%, což je navýšeno severovýchodním svahem Malé Mokrůvky. Průměrná hodnota sklonu pro celé území je 19%. Povodí bylo tachymetricky změřeno prostřednictvím totální stanice, polygonový pořad je zobrazen na obrázku v příloze 10.01, plastický reliéf celého povodí s vyznačenými body a uzávěrovým profilem je na obrázku 4.06. Nadmořské výšky a souřadnice změřených bodů získané z diplomové práce Bašty (2008) jsou uvedeny v příloze 10.02.
-36-
Obr. 4.06: Plastický reliéf zájmové oblasti s body polygonového pořadu (Bašta, 2008)
-37-
5. 5.1.
METODIKA Popis dat
Hodnoty jsou na povodí zaznamenávány každé 2 minuty. Pro modelování byly použity pouze průtoky a srážky. Časové řady z experimentálního povodí jsou pro další použití zpracovány do hodinových časových kroků. Vstupní data obdržená od Katedry vodního hospodářství a environmentálního modelování Fakulty životního prostředí byly ve formátu .xls, jejich ukázka je na obrázku 5.01. Obsahují datum, hodinu, referenční čas, teplotu ve stupních Celsia, průtok v litrech za sekundu a srážky v milimetrech.
Obr. 5.01: Data v hodinových krocích (Modrava 2, 1998).
Na povodí se měří od roku 1998 bez zimních pozorování, od roku 2006 už kontinuálně. Měření však provází mnohé výpadky způsobené poruchami přístrojů, vyčerpáním kapacity data-loggeru, vybitím baterek nebo nehodami, jako např. poškození místní zvěří (Pavlásek, 2011, in verb). Jako vstupní data byly vybrány dva různé datové soubory bez přestávek o délce 100 dní, tedy 2400 hodin. Pro kalibraci se jedná o období 8. 6. 2005 až 15. 9. 2005, pro validaci je to 22. 4. 2007 až 30. 7. 2007. Další hodnoty potřebné pro vytvoření modelu, jako půdní a geodetické charakteristiky jsou již popsány v předchozí kapitole.
-38-
5.2.
Sestavení modelu
Pro vytvoření 2D modelu příčného profilu povodí byl použit program HYDRUS 2D/3D, verze 1.12.0070, 3D-Standart, 2006. Hydrus je simulační software pro řešení proudění vody, tepla a rozpuštěných látek ve dvou nebo třídimenzionálním relativně nasyceném porézním prostředí. Proudění vody je řešeno Richardsovou rovnicí. Program zahrnuje MESHGEN2D pro tvorbu základní geometrie oblasti a její diskretizace na síť konečných prvků (Šimůnek & kol., 2006). K výpočtu hydraulických vlastností byla vybrána rovnice van Genuchtena v kombinaci s Mualemovým modelem.
5.2.1.
Doména a její výpočetní síť
Základní tvar profilu povodí s jeho sklonitostí byl určen dle geodetických bodů a jejich souřadnic uvedených v diplomové práci Bašty (2006). K vytvoření profilu byly použity body R2, 25, 26, 27, 28 a 9.0. Pro jejich umístění na povodí a souřadnice viz přílohy 10.01 a 10.02. Schéma s vyznačenými body a souřadnicemi je na obrázku 5.02. Celková délka domény je 429.015 metrů, s průměrným sklonem 21,4%. Hloubka profilu byla zvolena s ohledem na změřenou hloubku půdního profilu a tloušťky jednotlivých horizontů. Tento rozměr byl v průběhu sestavování modelu několikrát upravován, původních 0,8 m bylo při vytváření sítě pro stabilní výpočet a zkoušení různého rozložení spodní úrovně profilu rozšířeno až na 2 m, toto však bylo z důvodu nereálnosti a nadbytečnosti zavrženo a ve finální verzi má profil tloušťku 1 m.
Obr. 5.02: Schéma domény modelu
-39-
Výpočetní síť byla zpočátku také mnohokrát pozměňována. Prvotní ideou bylo, aby měl model co nejpodrobnější síť vystihující jednotlivé půdní horizonty, zvláště u povrchu byla bodová síť zvolena příliš podrobně. Takové řešení je ale zcela zbytečné, výpočet je zdlouhavý a často nestabilní. Při doméně o rozměrech téměř půl kilometru není možné mít uzly např. se vzdáleností 1 cm. Základní pravidla pro vytváření sítě uvádí Šimůnek & Šejna (2009). Síť tedy byla vytvořena jednodušeji, každý půdní horizont má ve své šířce 2 uzly sítě, na povrchu je tak vzdálenost uzlů 26 cm a na dně 62.5 cm, což pro tuto doménu úplně stačí. Celkem má vytvořená síť 7668 uzlů, schéma části sítě je na obrázku 5.03.
Obr. 5.03: Část výpočetní sítě
5.2.2.
Materiálové rozdělení
Parametry materiálových charakteristik jednotlivých půdních horizontů byli pomocí Rosetty implementované v Hydrusu stanoveny ze zrnitostního složení, které uvádí Petrželková (2011). Do těchto vypočítaných parametrů byly dosazeny změřené hodnoty nasycené hydraulické vodivosti a pórovitosti již zmiňované v kapitole 4.2.2. o hydropedologických charakteristikách. Profil byl nejprve rozdělen do čtyř horizontů Ah/E, Bs, Bhs, C se změřenými hodnotami nasycené hydraulické vodivosti. V této variantě byl simulovaný odtok do koryta zlomkem měřených hodnot, jelikož veškerá voda odtekla dolů volnou drenáží a nikoliv po svahu. Byla tedy přidána spodní málo propustná vrstva, její hydraulická vodivost byla upravována tak, aby propustila jen určité množství vody a simulovaný -40-
odtok po svahu odpovídal měřenému. Takovýto model ovšem ve svém složení zdaleka neodpovídal realitě na povodí. Podle Vogela (2004) je vhodné v modelu porézního prostředí vytvořit vrchní značně propustnou vrstvu, která odpovídá skutečné organické vrstvě o vysoké pórovitosti a hydraulické vodivosti na povrchu povodí, kterou prakticky nelze změřit (Jačka, 2013, in verb) a zároveň nahrazuje povrchový odtok. Konečný model má 4 horizonty: spodní C v hloubce cca 60 až 100 cm., nad ním je cca 30cm Bs+Bhs, dále horizont E o tloušťce přibližně 22cm s řádově nižší propustností, vrchních 8 cm představuje horizont O, ve kterém byla uvažována zmíněná vysoká hydraulická vodivost. Hodnoty vypočtených a zadaných parametru jsou uvedeny v tabulce 5.04. Tab. 5.04: Parametry hydraulických funkcí jednotlivých horizontů
horizont
θr
θs
α
n
KS
1 2 3 4
0,0277 0,0294 0,0354 0,0386
0,86 0,51 0,45 0,3866
4,82 4,64 2,61 4,12
1,5383 1,5308 1,4117 1,5743
K1 K2 0,00648 0,0432
Rozložení materiálu u vedených v tabulce 5.04 v uzlech výpočetní sítě je na obrázku 5.05.
Obr. 5.05: Materiálová distribuce ve výpočetní síti
-41-
5.2.3.
Počáteční podmínka
Zadaná počáteční podmínka, matematicky vyjádřena výrazem 5.01 definuje
vlhkostí prostředí 35%. Tato hodnota byla jen odhadem, na základě několika
simulací s různou počáteční vlhkostí bylo zjištěno, že ve výsledcích nehraje téměř žádnou roli. Její vliv se projevuje pouze cca prvních 20 dní, které do hodnocení výsledků a porovnání s měřeným odtokem nejsou zahrnuty. e, UC 0.35 fg he i Ω
5.01
Čas simulace byl nastaven na 2400 hodin, s maximálním krokem 1 hodina, minimálním 0,0001 hodin a počátečním 0,001 hodin. Pro každou hodinu jsou zadány měřené srážky.
5.2.4.
Okrajové podmínky
Dále bylo zapotřebí definovat okrajové podmínky. Vyznačení jejich umístění je na schématu v obrázku 5.06.
Obr. 5.06: Vyznačení okrajových podmínek
Okrajová podmínka Γ* na horní hranici je matematicky Neumannovou
podmínkou v čase proměnnou, v Hydrusu se nazývá atmosférická. Její hodnota se mění dle zadaných měřených vstupních srážek každou hodinu.
Na pravé spodní straně Γ? je výronová plocha, která reprezentuje odtok
z povodí do koryta Malé Mokrůvky.
-42-
Na dně byla zvolena homogenní Neumannova podmínka Γp volná drenáž,
kterou voda odtéká spodem na hladinu podzemní vody, nikoliv po svahu. Vlevo je oblast uzavřena nulovým tokem Γu .
Charakteristika a matematické vyjádření použitých okrajových podmínek je popsáno v kapitole 3.4.2. Okrajové podmínky.
5.2.5.
Kalibrace a validace
Čas simulace byl nastaven na 2400 hodin, s maximálním krokem 1 hodina, minimálním 0,0001 hodin a počátečním 0,001 hodin. Pro každou hodinu jsou zadány měřené srážky. vodivosti ve dvou vrchních horizontech 1* a 1? . Jsou sestaveny tři varianty modelu, Kalibrace byla provedena manuální parametrizací nasycené hydraulické
v nichž byla nastavena různá vodivost druhé vrstvy 1? a podle ní potom kalibrován
parametr 1* tak, aby rezidua měřeného a simulovaného souboru dat byly minimální. nejvyšších hodnotách. Parametr 1? byl určen ze změřené hodnoty a postupně řádově
Z nich byl poté vybrán nejvhodnější model vystihující hodnoty odtok i v jeho
zmenšován.
Výsledné simulované hodnoty odtoku byly získány z průtoku hranicí výronové plochy. Ten samozřejmě jako výstup z 2D modelu vychází v jednotkách m2/hod. Povodí je tímto způsobem konceptualizace zjednodušeno na obdélník, jehož strana má délku modelovaného 2D profilu 429.015 m a hodnoty je nutné vynásobit stranou druhou 372.95 m, která je vypočtena ze známé plochy zájmového území 0.16 km2. Výsledky byly pomocí R porovnány střední chybou ME a odmocninou střední kvadratické chyby RMSE (odhad rozptylu reziduí) a zpracovány do grafů výsledných hodnot a reziduí.
-43-
6.
VÝSLEDKY
Celkem byly sestaveny tři varianty modelu s různou nasycenou hydraulickou hodinu 1* a 1? jsou uvedeny v tabulce 6.01.
vodivostí ve dvou vrchních horizontech. Různé hodnoty parametrů v metrech za
Tab. 6.01: Hodnoty kalibrovaných parametru jednotlivých variant modelu
varianta
K1[m/hod]
K2[m/hod]
A B C
0.001188 0.01188 0.1188
17.065 88.910 298.710
Jako nejvhodnější byl vybrán model A, tedy s nejméně propustnou druhou vrstvou. Tento jako jediný vystihl maximální měřené průtoky a kritéria RMSE i ME vyšly pro validaci tohoto modelu nejlépe. Grafické výstupy modelu A zobrazené přímo v Hydrusu jsou na obrázcích 6.02 až 6.05. První obrázek, 6.02 zobrazuje 2D tok přes hranici výronové plochy v čase, který byl později pro porovnání vynásoben rozměrem povodí.
Obr. 6.02: Tok přes hranici výronové plochy
-44-
Obrázek 6.03 zobrazuje odtok spodem domény do podloží, tedy tok skrz hranici volné drenáže.
Obr. 6.03: Tok přes hranici volné drenáže
Na obr. 6.04 a 6.05 jsou pak kumulativní hodnotu odtoků výronovou plochou a volnou drenáží.
Obr. 6.04: Kumulativní tok přes hranici výronové plochy
-45-
Obr. 6.05: Kumulativní tok přes hranici volné drenáže
Vstupem do systému byli měřené srážky v metrech, graf jejich hodnot je zobrazen v příloze 10.03 pro kalibrační soubor z doby 28. 6. 2005 až 15. 9. 2005 a v příloze 10.06 pro soubor validační v období 12. 5. 2007 až 30. 7. 2007. Prvních 20 dní simulace je z porovnání vyjmuto. Výstupní hodnotou je odtok z povodí v metrech krychlových za hodinu, tedy hodnota již vynásobena rozměrem povodí. V grafické příloze 10.04 a 10.07 jsou vyneseny simulované i měřené odtoky, v příloze 10.05, 10.08 pak jejich rezidua s vypočtenými kritérii ME a RMSE. Další varianty modelu B a C nevystihují vrcholy hodnot průtoku tak jako A. Jejich grafický výstup pro kalibrační i validační data je zobrazen v přílohách 10.09 až 10.16. Hodnotící kritéria jednotlivých modelů jsou v tabulce 6.04 uvedena pro kalibraci i validaci. Tab. 6.04: Hodnoty kritérií kalibrace a validace
kalibrace validace kalibrace model B validace kalibrace model C validace model A
-46-
ME 0,000351 0,106756 -0,000236 0,275878 -0,000007 0,228679
RMSE 0,063178 0,077332 0,062570 0,101327 0,061306 0,107915
7.
DISKUZE
Hlavním cílem této práce bylo vytvořit konceptuální model povodí ve formě 2D příčného profilu. Profil byl sestaven na základě informací o geodetickém a hydropedologickém charakteru zájmové oblasti. Po získání výsledků simulace a jejich úpravě bylo možné posoudit srážkoodtokový proces na povodí. Objem infiltrovaných vstupních srážek se v půdním profilu rozdělí na odtok po svahu a odtok dolů volnou drenáží. Graf na obrázku 6.03 zobrazuje odtok spodem domény do podloží, obrázek 6.02 množství vody, které odteče po svahu a pak výronovou plochou do toku. Odtok do podloží má daleko plynulejší průběh. Dle kumulativních hodnot na obrázku 6.04 a 6.05 je odtok do podloží několikanásobně vyšší než ten po svahu. Podle Bayera a kol. (2004) se tyto objemy rozdělují cca na poloviny. Poznamenal však, že poměr odtoků se podle probíhající fáze vodního režimu půd může měnit. Výsledky modelu byly hodnoceny na základě kritérií ME (střední chyba) a RMSE (odmocnina střední kvadratické chyby). Pomocí ME byl model kalibrován, shoda validačních měřených a simulovaných dat byla podle RMSE i ME nejlepší pro model A, v modelech B a C byla střední chyba ME více než dvakrát větší. Chování modelu lze posoudit také vizuálně z grafického zpracování v příloze 10.05-10.16. Model je schopen vystihnout způsob reakce povodí na probíhající srážku, varianta A dosahuje maximální měřené hodnoty průtoku uzávěrovým profilem. Na povodí jsou však preferenční cesty a do celkových hodnot zasahuje i povrchový odtok, který se může projevovat jinak a daleko prudčeji než průtok porézním prostředím.
-47-
8.
ZÁVĚR
Experimentální povodí Malá Mokrůvka neboli Modrava 2 bylo deterministicky modelováno Richardsovou rovnicí v programu Hydrus 2D. Ze zjištěných hydropedologických a geodetických charakteristik byl sestaven konceptuální model příčného profilu povodí. Metoda konceptualizace je v případě proudění v porézním prostředí velice efektivní, protože skutečný pórový systém je reálně neřešitelný a není ani možné provést natolik detailní výzkum, aby vystihl jeho komplikovanost. Schematizací oblasti i proudění dojde ke zjednodušení, díky kterému je možné proudění řešit. Nezbytnou
součástí
konceptuálního
modelování
je
správná
aplikace
počátečních a okrajových podmínek, určení příslušných parametrů konstitutivních vztahů definujících porézní materiál a vhodná prostorová i časová diskretizace. Výsledky jsou tedy přibližné, lze je však použít pro posouzení objemové bilance půdního profilu díky informaci o objemovém toku definovanými hranicemi a k vyhodnocení srážkoodtokových procesů na povodí. Shoda pozorovaných a modelovaných odtoků u povodí není perfektní, nicméně podle zobrazených výsledků je 2D model půdního profilu schopen simulovat základní charakter reakce svahu experimentálního povodí pomocí vrchní více propustné vrstvy oddělené velmi málo propustným horizontem.
-48-
9.
POUŽITÁ LITERATURA
BAŠTA P., 2008: Digitální model terénu povodí Modrava 2. Diplomová práce ČZU, Praha, vedoucí práce P. Máca. BAYER T., TESAŘ M., ŠÍR M., 2004: Runoff formation in a small mountaineous watershed. Aktuality šumavského výzkumu II: s. 56-62. BEAR J., CHENG A. H-D., 2010: Modeling groundwater flow and contaminant transport - Theory and applications of transport in porous media. Springer, London, 850s., ISBN 14-020-6682-1. BEAR, J., 1972: Dynamics of Fluids in Porous Media. American Elsevier, New York, 764 s. BROOKS R. H., COREY A. T., 1964: Hydraulic properties of porous media. Hydrology papers 3. Colorado State University, Fort Collins in Colorado, 27s. CÍSLEROVÁ M., 1989: Inženýrská hydropedologie. ČVUT, Praha, 156s. CÍSLEROVÁ M., VOGEL T., 1998: Transportní procesy. ČVUT, Praha, 182s, ISBN 80-01-01866-0. HÁLEK V., ŠVEC J., 1973: Hydraulika podzemní vody. Academia, Praha, 376 s. JAČKA L., PAVLÁSEK J., 2010: Vybrané hydropedologické charakteristiky podzolů v centrální oblasti NP Šumava. Vodohospodářské technicko-ekonomické informace (VTEI) 52/5: s. 17-19. JAČKA L., PAVLÁSEK J., JINDROVÁ M., BAŠTA P., ČERNÝ M., BALVÍN A., PECH P., 2012: Steady infiltration rates estimated for a mountain forest catchment based on the distribution of plant species. Journal of forest science 58/12: s. 536-544. JINDROVÁ M., 2011: Vymezení charakteristických hydropedologických jednotek na povodí Modrava 2 na základě fytocenologického průzkumu. Diplomová práce ČZU, Praha, vedoucí práce J. Pavlásek. KAZDA I., 1983: Proudění podzemní vody. Řešení metodou konečných prvků. SNTL, Praha, 232 s. KODEŠOVÁ R., 2012: Modelování v pedologii, 2. vydání. ČZU, Praha. 152s, ISBN 80-213-1347-1. KUNA P., KUŘÍK P., 1998: Hydrometeorologická měření na Šumavě. Zprávy lesnického výzkumu 43/3-4: s. 37-39. -49-
KURÁŽ M., 2011: Numerical solution of the flow and transport equations in porous media with the dual permeability conceptual approach. VeRBuM, Zlín, 178s, ISBN 978-80-87500-12-5. KUTÍLEK M., 1978: Vodohospodářská pedologie, 2. přepracované vydání. SNTL, Praha, 295 s, ISBN 80-213-1347-1. KUTÍLEK M., 1984: Vlhkost pórovitých materiálů. SNTL, Praha, 210s. KUTÍLEK M., KURÁŽ V., CÍSLEROVÁ M., 2004: Hydropedologie 10, 2. přepracované vydání. ČVUT, Praha, 176s, ISBN 978-80-01-02237-5. LAL R., SHUKLA K. M., 2004: Principles of Soil Physics. Marcel Dekker Inc., New York, 682 s, ISBN 0-203-02123-1. MLS J., 1984: Hydraulika podzemní vody. ČVUT, Praha, 157s. MUCHA I., ŠESTAKOV V., 1987: Hydraulika podzemných vôd. Alfa, Bratislava spolu se SNTL, Praha, 344 s. NĚMEČEK J., 2001: Taxonomický klasifikační systém půd české republiky. ČZU spolu s VÚMOP, Praha. 79s. ISBN 80-238-8061-6. PACHEPSKY Y., GUBER A. K., VAN GENUCHTEN M. TH., NICHOLSON T. J., CADY R. E., ŠIMŮNEK J., SCHAAP M. G., 2006: Model abstraction techniques for soil-water flow and transport. Division of Fuel, Engineering & Radiological Research, Office of Nuclear Regulatory Research, US Nuclear Regulatory Commission PAVELKOVÁ H., DOHNAL M., VOGEL T., 2012: Hillslope runoff generation – comparing different modeling approaches. Journal of Hydrology and Hydromechanics 60/2: s. 73-86. PAVLÁSEK J., 2010: Retenční schopnosti malého horského povodí při extrémních srážkoodtokových . Vodohospodářské technicko-ekonomické informace (VTEI) 52/5: s. 12-14. PAVLÁSEK J., MÁCA P., ŘEDINOVÁ J., 2006: Analýza hydrologických dat z Modravských povodí. Journal of Hydrology and Hydromechanics 54/2: s. 207-216. PAVLÁSEK, J., ŘEDINOVÁ, J., SKALSKÁ P., 2009: Evaluation of Monitoring on Modrava Catchments. Soil & Water Resources 4: s. 66-74. PECH P., 2010: Speciální případy hydrauliky podzemních vod. ČZU ve VÚV T. G. Masaryka, Praha, 103s, ISBN 978-80-87402-04-7. -50-
PETRŽELKOVÁ V., 2011: Možnosti stanovení zrnitostního složení půd. Bakalářská práce ČZU, Praha, vedoucí práce J. Pavlásek. ŘÍHA J., 1997: Matematické modelování hydrodynamických a disperzních jevů. VUT, Brno, 185 s, ISBN 80-214-0827-8. SCHAAP M. G., LEIJ F. J., VAN GENUCHTEN M. TH., 2001: Rosetta: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. Journal of Hydrology 251/3-4: s 163-176. ŠIMŮNEK J., ŠEJNA M., 2009: Notes on Spatial and Temporal Discretization (when working with Hydrus). Dostupné z http://www.pc-progress.com. ŠIMŮNEK J., VAN GENUCHTEN M. TH., ŠEJNA M., 2006: The HYDRUS Sofware Package for Simulating the Two- and Three-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media - Technical Manual, Version 1.0. Dostupné z http://www.pc-progress.com. TESAŘ M., ŠÍR M., SYROVÁTKA O., PRAŽÁK J., LICHNER L., KUBÍK F., 2001: Soil Water Regime In Head Water Regions–Observation, Assessment And Modelling. Journal of Hydrology and Hydromechanics 49/6: s. 355-375. TOURKOVÁ, J. 1996: Hydrogeologie. ČVUT, Praha, 165 s., ISBN 80-01-03101-2. VALENTOVÁ J., 2001: Hydraulika podzemní vody, 2. přepracované vydání. ČVUT, Praha, 174 s, ISBN 80-010-2404-0. VAN GENUCHTEN M. TH., 1980: A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44: s. 892-898. VAN GENUCHTEN M. TH., LEIJ F. J., YATES S. R., 1991: The RETC Code for Quantifying the Hydraulic Functions o Unsaturated Soils. Version 1.0. EPA Report 600/2-91/065, U.S. Salinity Laboratory, USDA, ARS, Riverside, California, 93s. VOGEL T., 2004: Modeling Flow and Transport in Natural Porous Media. ČVUT, Praha, 30 s, ISBN 80-01-02970-0.
-51-
10. PŘÍLOHY Seznam příloh 10.01 Schéma polygonového pořadu 10.02 Souřadnice a nadmořské výšky bodů polygonového pořadu 10.03 Graf měřených srážek kalibračního souboru 10.04 Měřené a simulované odtoky kalibračního souboru 10.05 ModelA:Rezidua měřených a simulovaných hodnot kalibračního souboru 10.06 ModelA:Graf měřených srážek validačního souboru 10.07 ModelA:Měřené a simulované odtoky validačního souboru 10.08 ModelA:Rezidua měřených a simulovaných hodnot validačního souboru 10.09 ModelB:Měřené a simulované odtoky kalibračního souboru 10.10 ModelB:Rezidua měřených a simulovaných hodnot kalibračního souboru 10.11 ModelB:Měřené a simulované odtoky validačního souboru 10.12 ModelB:Rezidua měřených a simulovaných hodnot validačního souboru 10.13 ModelC:Měřené a simulované odtoky kalibračního souboru 10.14 ModelC:Rezidua měřených a simulovaných hodnot kalibračního souboru 10.15 ModelC:Měřené a simulované odtoky validačního souboru 10.16 ModelC:Rezidua měřených a simulovaných hodnot validačního souboru
-52-
Příloha 10.01 Schéma polygonového pořadu (Bašta, 2008)
Příloha 10.02 Souřadnice a nadmořské výšky bodů polygonového pořadu (Bašta, 2008)
-53-
10.03 Graf měřených srážek kalibračního souboru
10.04 ModelA:Měřené a simulované odtoky kalibračního souboru
-54-
10.05 ModelA:Rezidua měřených a simulovaných hodnot kalibračního souboru
10.06 Graf měřených srážek validačního souboru
-55-
10.07 ModelA:Měřené a simulované odtoky validačního souboru
10.08 ModelA:Rezidua měřených a simulovaných hodnot validačního souboru
-56-
10.09 ModelB:Měřené a simulované odtoky kalibračního souboru
10.10 ModelB:Rezidua měřených a simulovaných hodnot kalibračního souboru
-57-
10.11 ModelB:Měřené a simulované odtoky validačního souboru
10.12 ModelB:Rezidua měřených a simulovaných hodnot validačního souboru
-58-
10.13 ModelC:Měřené a simulované odtoky kalibračního souboru
10.14 ModelC:Rezidua měřených a simulovaných hodnot kalibračního souboru
-59-
10.15 ModelC:Měřené a simulované odtoky validačního souboru
10.16 ModelC:Rezidua měřených a simulovaných hodnot validačního souboru
-60-
