Česká zemědělská univerzita v Praze Fakulta lesnická a environmentální Katedra ekologie a životního prostředí
VYUŽITÍ MATEMATICKÉHO MODELOVÁNÍ PŘI HODNOCENÍ EKOLOGICKÝCH ŠKOD NA PODZEMNÍCH VODÁCH (na příkladu vybrané rafinérské společnosti)
DIPLOMOVÁ PRÁCE
Vedoucí diplomové práce: Doc. Ing. RNDr. Ivana Landa, DrSc. Diplomant: Jan Baier
2007
1
Prohlášení
Prohlašuji, že jsem diplomovou práci na téma „Využití matematického modelování při hodnocení ekologických škod na podzemních vodách (na příkladu vybrané rafinérské společnosti)“ vypracoval samostatně pod vedením Doc. Ing. RNDr. Ivana Landy, DrSc. a s použitím uvedené literatury.
V Praze dne ………………..
Podpis diplomanta ………………
2
Poděkování
Děkuji rodině za podporu, Ivanu Landovi za poskytnutí podkladových materiálů, odborné vedení a jeho čas strávený při dlouhých konzultacích.
3
OBSAH
Seznam obrázků ...................................................................................................5 Seznam grafů........................................................................................................8 Seznam tabulek ....................................................................................................9 1. Úvod................................................................................................................10 2. Cíle diplomové práce ......................................................................................11 3. Rozdělení a popis modelů...............................................................................11 3.1 Fyzikální (měřítkové)................................................................................11 3.2 Analogové................................................................................................11 3.3 Matematické principy ...............................................................................12 3.3.1 Popis matematických modelů ..........................................................12 3.3.2 Rozdělení matematických modelů...................................................12 3.3.3 Popis a druhy jednotlivých modelů ..................................................13 3. 4 Popis programu Processing Modflow Pro ...............................................14 3.4.1 Historie a popis Processing Modflow Pro ........................................14 3.4.2 Součásti Processing Modflow Pro ...................................................15 4. Numerické metody ..........................................................................................16 4.1 Metoda konečných diferencí ....................................................................17 4.2 Metoda konečných prvků .........................................................................20 5. Základní rovnice a děje popisující proudění podzemní vody a transportu znečištění.......................................................................................................21 5.2 Darcyho zákon.........................................................................................21 5.3 Rovnice kontinuity a třírozměrného proudění ..........................................22 5.4 Šíření znečištění ......................................................................................23 5.4.1 Rovnice popisující šíření znečištění ................................................24 5. 5 Okrajové a počáteční podmínky .............................................................25 5.5.1 Hranice s předepsanou hodnotou hydraulické výšky.......................25 5.5.2 Hranice s předepsaným tokem........................................................25 5.5.3 Polopropustná hranice.....................................................................25 5.5.3 Počáteční podmínky ........................................................................26 6. Testování použitelnosti programu. Lokalita Kralupy........................................26 6. 1 Ekologie petrochemických podniků.........................................................26 6.2 Podmínky provozu podniku......................................................................28
4
6. 3 Ekologické a přírodní poměry areálu Kaučuk Group a.s., Kralupy nad Vltavou a v širším okolí ...........................................................................28 6.3.1 Vymezení areálu Kaučuk Group a.s. ...............................................28 6.3.2 Topologické poměry ........................................................................30 6.3.3 Klimatické a hydrologické poměry ..................................................30 6.3.4 Geomorfologické a geologické poměry ...........................................30 6.3.5 Proudění podzemních vod...............................................................31 6.3.7 Využití vodních zdrojů .....................................................................32 6.3.8 Zdroje a druhy kontaminace ............................................................32 7. Modelové řešení..............................................................................................33 7.1 Informační zabezpečení modelu a vyhodnocení dat................................33 7.2 Problémy při sestavování modelu............................................................34 7.3 Vymezení modelovaného území..............................................................35 7.4 Sestavení rámcového modelu .................................................................35 7.5 Stanovení okrajových a počátečních podmínek.......................................37 7.6 Diskretizace parametrů prostředí.............................................................38 7.7 Modelování stávajícího stavu...................................................................38 7.8 Kalibrace a optimalizace modelu ............................................................39 7.9 Modelování šíření znečištění ...................................................................43 7.9.1 Popis Modelu...................................................................................43 7.9.2 Definované vlastnosti modelu..........................................................43 7.9.3 Výsledky modelu MT3DMS .............................................................44 8. Modelované situace ........................................................................................45 8.1 Zobrazení a popis modelovaných situací.................................................46 9. Závěr...............................................................................................................70 10. Seznam literatury ..........................................................................................73 11. Seznam Příloh...............................................................................................75
5
Seznam obrázků Obr.1 Darcyho pokus Obr.2 Areál Kaučuk Group a.s.,Kralupy nad Vltavou Obr.3 Modelované území Obr 4. Dialogové okno definování výpočetní sítě modelu Obr 5. Modelované území (obrázek je ukázkou jednoho z výstupu modelování) Obr. 6 Izolinie hladiny podzemní vody Obr. 7 Hladina podzemní vody ovlivněná HOPV Obr. 8 Znečištění po ročním kontinuálním úniku Obr.9 hydroizohypsy bez čerpání z HOPV Obr.10 Izolinie šíření kontaminantů bez čerpání HOPV po 1 měsíci Obr.11 Izolinie šíření kontaminantů bez čerpání HOPV po 6 měsících Obr.12 Izolinie šíření kontaminantů bez čerpání HOPV po 12 měsících Obr.13 Izolinie šíření kontaminantů bez čerpání HOPV po 36 měsících Obr.14 hydroizohypsy při 100 % čerpání z HOPV Obr.15 Izolinie šíření kontaminantů při 100 % čerpání HOPV po 1 měsíci Obr.16 Izolinie šíření kontaminantů při 100 % čerpání HOPV po 6 měsících Obr.17 Izolinie šíření kontaminantů při 100 % čerpání HOPV po 12 Obr.18 Izolinie šíření kontaminantů při 100 % čerpání HOPV po 24 měsících Obr.19 hydroizohypsy při 80 % čerpání z HOPV Obr.20 Izolinie šíření kontaminantů při 80 % čerpání HOPV po 1 měsíci Obr.21 Izolinie šíření kontaminantů při 80 % čerpání HOPV po 6 měsících Obr.22 Izolinie šíření kontaminantů při 80 % čerpání HOPV Obr.23 Izolinie šíření kontaminantů při 80 % čerpání HOPV Obr.24 Hydroizohypsy při 50 % čerpání z HOPV Obr.25 Izolinie šíření kontaminantů při 50 % čerpání HOPV po 1 měsíci Obr.26 Izolinie šíření kontaminantů při 50 % čerpání HOPV po 6 měsících Obr.27 Izolinie šíření kontaminantů při 50 % čerpání HOPV po 6 měsících Obr.28 Izolinie šíření kontaminantů při 50 % čerpání HOPV Obr.29 Hydroizohypsy při 20 % čerpání z HOPV Obr.30 Izolinie šíření kontaminantů při 20 % čerpání HOPV po 1 měsíci Obr.31 Izolinie šíření kontaminantů při 20 % čerpání HOPV po 6 měsících Obr.32 Izolinie šíření kontaminantů při 20 % čerpání HOPV po 12 měsících Obr.33 Izolinie šíření kontaminantů při 20 % čerpání HOPV po 24 měsících Obr.34 hydroizohypsy při 150 % čerpání z HOPV Obr.35 Izolinie šíření kontaminantů při 150 % čerpání HOPV po1 měsíci Obr.36 Izolinie šíření kontaminantů při 150 % čerpání HOPV po 6 měsících Obr.37 Izolinie šíření kontaminantů při 20 % čerpání HOPV po 6 měsících Obr.38 Izolinie šíření kontaminantů při 150 % čerpání HOPV po 6 měsících Obr.39 Hydroizohypsy při čerpání pouze z kralupské větve HOPV Obr.40 Izolinie šíření kontaminantů čerpání kralupské větve HOPV po1 měsíci
6
Obr.41 Izolinie šíření kontaminantů čerpání kralupské é větve HOPV po 6 měsících Obr.42 Izolinie šíření kontaminantů čerpání kralupské větve HOPV po12 měsících Obr.43 Izolinie šíření kontaminantů čerpání kralupské větve HOPV po 24 měsících Obr.44 Hydroizohypsy při čerpání pouze z veltruské větve HOPV Obr.45 Izolinie šíření kontaminantů čerpání veltruské větve HOPV po1 měsíci Obr.46 Izolinie šíření kontaminantů čerpání veltruské větve HOPV po 6 měsících Obr.47 Izolinie šíření kontaminantů čerpání veltruské větve HOPV po12 měsících Obr.48 Izolinie šíření kontaminantů čerpání veltruské větve HOPV po 24 měsících Obr.49 Hydroizohypsy při čerpání 150 % z veltruské větve HOPV Obr.50 Izolinie šíření kontaminantů čerpání 150 % z veltruské větve HOPV po 1 měsíci Obr.51 Izolinie šíření kontaminantů čerpání 150 % z veltruské větve HOPV po 6 měsících Obr.52 Izolinie šíření kontaminantů čerpání 150 % z veltruské větve HOPV po 12 měsících Obr.53 Izolinie šíření kontaminantů čerpání 150 % z veltruské větve HOPV po 24 měsících Obr.54 Hydroizohypsy při čerpání 100 % z kladrubské a 80 % z veltruské větve HOPV Obr.55 Izolinie šíření kontaminantů čerpání 100 % z kladrubské a 80 % veltruské
větveHOPV
po1 měsíci Obr.56 Izolinie šíření kontaminantů čerpání 100 % z kladrubské a 80 % veltruské větve HOPV po 6 měsících Obr.57 Izolinie šíření kontaminantů čerpání 100 % z kladrubské a 80 % veltruské větve HOPV po 12 měsících Obr.58 Izolinie šíření kontaminantů čerpání 100 % z kladrubské a 80 % veltruské větve HOPV po 24 měsících Obr.59 Hydroizohypsy při čerpání 100 % z veltruské a 80 % z kladrubské větve HOPV Obr.60 Izolinie šíření kontaminantů čerpání 80 % z kladrubské a 100 % veltruské větve HOPV po1 měsíci Obr.61 Izolinie šíření kontaminantů čerpání 80 % z kladrubské a 100 % veltruské větve HOPV po 6 měsících Obr.62 Izolinie šíření kontaminantů čerpání 80 % z kladrubské a 100 % veltruské větve HOPV po 12 měsících Obr.63 Izolinie šíření kontaminantů čerpání 80 % z kladrubské a 100 % veltruské větve HOPV po 24 měsících Obr.64 Hydroizohypsy při simulaci zvýšené hladiny podzemní vody Obr.65 Izolinie šíření kontaminantů po 1 měsíci při simulaci zvýšeni HPV Obr.66 Izolinie šíření kontaminantů po 6 měsících při simulaci zvýšeni HPV Obr.67 Izolinie šíření kontaminantů po 12 měsících při simulaci zvýšeni HPV Obr.68 Izolinie šíření kontaminantů po 24 měsících při simulaci zvýšeni HPV
7
Seznam grafů Graf 1. Porovnání pozorovaných a vypočítaných hodnot v při zadání odhadovaných hodnot parametrů. Graf 2. Porovnání pozorovaných a vypočítaných hodnot při Cbot = 0.000001 m.s
-1
Graf 3. Porovnání vypočítaných a vypočítaných hodnot při zadání nakalibrovaných parametrů Graf 4. Porovnání pozorovaných a vypočítaných hodnot k 17.2. 2006 Graf 5. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech bez čerpání z HOPV Graf 6. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech při čerpání 100 % z objektů HOPV Graf 7. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech pří čerpání 80 % z objektů HOPV Graf 8. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech pří čerpání 50 % z objektů HOPV Graf 9. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech pří čerpání 20 % z objektů HOPV Graf 10. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech pří čerpání 150 % z objektů HOPV Graf 11. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech při čerpání pouze z kralupské větve HOPV Graf 12. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech při čerpání pouze z veltruské větve HOPV Graf 13. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech při čerpání 150 % z veltruské větve HOPV Graf 14. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech při čerpání 100 % z kladrubské a 80 % z veltruské větve HOPV Graf 15. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech při čerpání 100 % z veltruské a 80 %z kladrubské větve HOPV Graf 16. Změny hladin v pozorovacích vrtech při zvýšení hladiny podzemní vody
8
Seznam tabulek Tab. 1 Průměrné měsíční úhrny srážek (1931 – 1960) Tab. 2 Definované parametry řeky Tab. 3 Definované parametry GBH Tab. 4 Parametry prostředí Tab. 5 Pozorované a spočítané hodnoty při Cbot = 0.000001 m.s
-1
Tab. 6 Odhadované parametry pomocí PEST Tab. 7 Vypočítané a naměřené hodnoty Tab. 8 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech bez čerpání z HOPV Tab. 9 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech při čerpání 100 % z objektů HOPV Tab. 10 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech pří čerpání 80 % z objektů HOPV Tab. 11 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech pří čerpání 50 % z objektů HOPV Tab. 12 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech pří čerpání 20 % z objektů HOPV Tab. 13 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech pří čerpání 150 % z objektů HOPV Tab. 14 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech při čerpání pouze z kralupské větve HOPV Tab. 15 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech při čerpání pouze z veltruské větve HOPV Tab. 16 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech při čerpání 150 % z veltruské větve HOPV Tab. 17 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech 100 % z kladrubské a 80 % z veltruské větve HOPV Tab. 18 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech 100 % z veltruské a 80 % z kladrubské větve HOPV
9
1. Úvod
Vlivem růstu počtu obyvatel dochází k rozsáhlému antropogennímu ovlivnění a nepříznivým změnám v životním prostředí. S tím souvisí i nadměrné využívání přírodních zdrojů, které vede k poklesu zásob až k jejich vyčerpání. Mezi ekologicky nevýznamnější a tudíž i nejvyužívanější přírodní zdroj patří voda. Její nadměrné využívání vede na mnoha místech planety nejen k jejímu nedostatku, ale také k rozsáhlému znečištění vodních zdrojů. Z ekologického a ekonomického hlediska patří k nejvyužívanějším zdrojům pitné vody, voda podzemní, neboť se vyznačuje relativně stabilní kvalitou. Proto je důležité regulovat její využívání a zvýšit spolehlivost její ochrany. Mezi největší znečišťovatele podzemních vod v průmyslových zemích patří průmysl, a to jak na úrovni lokální (chemický), subregionální (těžba nerostných surovin), tak i regionální (zemědělství). Jedním z největších znečišťovatelů povrchových i podzemních vod je chemický průmysl, kdy dochází jak k produkci značného
množství
odpadních
vod,
tak
k nekontrolovatelným
únikům
znečišťujících látek do podpovrchových vod. Chemický průmysl působí nejen jako jeden z největších znečišťovatelů, ale také odběratelů povrchových i podzemních vod. Existuje mnoho prognózních metod pro posuzování vlivu úniků znečištění na jakost vod. Mezi nejdůležitější patří metody matematického modelování. Existuje mnoho matematických programů, lišících se svými algoritmy a použitými
výpočetními
metodami
(postupy),
uživatelským
prostředím,
možnostmi a podmínkami jejich použití. Často se stává rozhodující jejich cena. Výhodné je, že kromě relativně drahých programů, které jsou z finančního hlediska pro řadu pracovišť téměř nedostupné, existují i takzvané free programy1, mezi které patří také Processing modflow Pro, který je od verze 5.3 na stránkách www.pmwin.net volně ke stažení. Program vychází z verze MODFLOW vyvíjené od sedmdesátých let v programovacím jazyce Fortran. K řešení úlohy ovšem používám školní verzi Processing modflow Pro 7.0.0 z roku 2002, a 3D nadstavbu 3D Master jež dohromady stojí 1595.00 $ 2.
1 2
Programy volně šířené, volně dostupné, např. na internetu Což odpovídá ceně cca 30 000 Kč
10
2. Cíle diplomové práce
V rámci diplomové práce se soustředím na možnost využití metodiky matematického modelování hydraulických a ekologických úloh na příkladu provozu reálného systému ochrany podzemních vod pomocí hydraulické bariéry. Jde o areál petrochemického areálu Kaučuk Group a.s., Kralupy nad Vltavou. K tomuto úkolu jsem využít program Processing modflow Pro a na základě výsledků modelování navrhnout a doporučit režim provozu hydraulické bariéry. Ze zadání též vyplývá, že cílem je i zhodnocení vlivu klimatických podmínek a dále optimalizace čerpaná množství v závislosti na úniku znečištění z rafinérie. Práce je zaměřena na zhodnocení čerpaných množství vody z hydraulické bariéry v závislosti na již zjištěném poklesu využitelného množství podzemních vod v objektech využívaných pro individuální zásobování pitnou vodou. V rámci řešení této úlohy aplikovaného výzkumu jsem měl za cíl dokonale ovládnout program Processing modflow Pro a popsat jeho výhody a nedostatky při tvorbě modelu.
3. Rozdělení a popis modelů
3.1 Fyzikální (měřítkové)
Fyzikální
modely
simulují
části
přirozeného
světa
v laboratorních
podmínkách. Používají se například pro simulaci srážko – odtokových vztahů, nebo v říční hydraulice při simulaci zatížení vodních staveb. (Dingman, 1998)
3.2 Analogové Analogové modely jsou fyzikální simulace studovaného systému. Pomocí analogového modelu se například může simulovat tok kapaliny nahrazený elektřinou nebo teplem. Možnost simulace vychází z podobnosti Darcyho a Ohmova
a
Furierova
zákona
pro
šíření
fyzikálních
polí
a
to
pole
piezometrického, elektrického a tepelného. Tato analogie umožňuje využít i jednotlivých polí pro simulaci polí analogických, např. elektrického a tepelného pole pro řešeni hydrodynamických úloh. 11
3.3 Matematické principy
3.3.1 Popis matematických modelů
V posledních desetiletích došlo ve většině inženýrských disciplin k rychlému rozvoji moderních vědních oborů založených na počítačové technice. V souhlase se současnými trendy se využívají jako moderní prostředky k popisu vodní bilance v hydrologii simulační modely. Na rozdíl od klasických metod umožňují komplexní přístup k řešené problematice, podrobnou analýzu působení všech důležitých faktorů a vyzkoušení mnoha variant řešení včetně očekávaných dlouhodobých výsledků. Simulační modely transportních procesů na tyto modely navazují a je jimi možno simulovat migraci znečištění, jeho rozpad a posuzovat mnoho variant ekologických scénářů. (Vogel, Císlerová, 1998) Matematické modely jsou explicitní sekvence souboru na sebe navazujících rovnic, numerických a logických kroků, které transformují numerické vstupy do numerických výstupů (rovnice, veličiny, parametry, algoritmus). (Dingman, 1998)
3.3.2 Rozdělení matematických modelů
Matematické modely jak popisují (Hrádek, Kuřík 2002, Kovář 1990, Beven 2001), se dělí na modely stochastické, deterministické a smíšené. •
Stochastické
modely
představují
skupinu
modelů,
které
se
dají
charakterizovat absencí vazebnosti mezi příčinou a následkem popisovaného jevu v rámci charakterizovaného systému. •
Deterministické modely mají za úkol popisovat pomocí matematických
vztahů fyzikální systém. Přesnost popisu fyzikálního systému modelem se může zvyšovat s ohledem na kvalitu vstupních dat, protože se stoupající přesností popisu stoupají i nároky na vstupní data. Podle kvantity a kvality pozorovaných proměnných
a
odvozených
parametrů
se
ustálilo
základní
rozdělení
deterministických modelů do dvou skupin: a) hydrologické modely (také označovány jako parametrické nebo srážkoodtokové) b) hydrodynamické modely (Deterministic, hydrodynamic Laws - DL modely, 12
White Box) fyzikálně popisují realitu nejvěrněji. Respektují principy zachování hmoty, hybnosti a energie. Jsou to modely s geometricky rozdělenými
parametry,
které
popisují
řešené
procesy
pomocí
diferenciálních rovnic. Struktura systému je u hydrodynamických modelů vložena přímo do základních rovnic. Modely mohou popisovat vybrané dílčí
hydrologické
procesy
(komponentní
modely)
hydrologické procesy v povodí (komplexní modely).
nebo
všechny
Do této skupiny
modelů patří i mnou požívaný program Processing modflow Pro. •
Mezi stochastickými a deterministickými modely dochází k jistému překrývání a dostáváme tím modely smíšené. Tyto modely pak obsahují submodely stochastické i deterministické povahy a sestavují se pro zdokonalení výstupů deterministického modelu.
3.3.3 Popis a druhy jednotlivých modelů
AQUA3D Řeší trojrozměrné proudění podzemní vody a šíření kontaminantů pomocí metody konečných prvků – GFEM. Řeší systém rovnic popisující proudění podzemní vody a šíření pro homogenní izotropní zvodnělé vrstvy i pro nehomogenní a anizotropní. Model šíření kontaminantů je plně integrován s modelem proudění - používá stejné datové vstupní a výstupní formáty. GMS – Groundwater Modeling Systém GMS představuje sadu nástrojů pro každou fázi simulace proudění podzemní vody. Model umožňuje charakterizovat lokalitu, kalibraci parametrů, následné zpracování a vizualizaci výsledků. GMS podporuje, jak metodu konečných diferencí, tak metodu konečných prvků v 2D i 3D zobrazení, zahrnující MODFLOW 2000, MODPATH, MT3DMS/RT3D, SEAM 3D, ART3D, UTCHEM, FEMWATER, PEST, MODAEMA a SEEP2D. FEFLOW Představuje jeden z nejpropracovanějších souborů programů pro modelování proudění podzemní vody a transportních procesů v porézním prostředí za nasycených a nenasycených podmínek. FEFLOW je simulační systém založený na metodě konečných prvků, který obsahuje interaktivní grafické rozhraní, nástroje pro regionalizaci a vizualizaci dat a účinné numerické techniky.
13
Poskytuje nástroje pro tvorbu sítě konečných prvků, přiřazení modelových parametrů a okrajových podmínek, průběh simulace a vizualizaci výsledků.
SWMS_3D SWMS_3D je počítačový program pro simulaci pohybu vody a rozpuštěných látek v
trojrozměrném proměnném saturovaném prostředí. Program řeší
proudění podzemní vody v nasyceném, nenasyceném prostředí pomocí Richardsonovy rovnice a šíření rozpuštěných látek pomocí konvekčně-disperzní rovnice. Řídící rovnice pro transport a proudění jsou numericky řešeny pomocí Galerkinovy metody lineárních konečných prvků. Modflow surfact Modflow Surfact je souhrnný nástroj pro 3D modelování na základě metody konečných diferencí, založených na USGS modflow kódu, který představuje nejvíce používaný kód pro proudění podzemní vody. Do tohoto nástroje byly dále zpracovány další výpočetní moduly, což výrazně zvyšuje možnosti simulací. Modflow Surfact v sobě integruje moduly proudění i šíření u kontaminantů. Je podporován ve dvou grafických rozhraních: Visual Modflow Pro a Groundwater Vistas. Visual Modflow Pro Visual
Modflow představuje
ověřený
standart
pro
profesionální
3D
modelování proudění podzemní vody a transportu kontaminantů za použití komponent Modflow-2000, Modpath, MT3MS, RT3D. Groundwater Vistas Groundwater Vistas představuje propracované Windows grafické rozhraní pro 3D modelování podzemní vody a transportu kontaminantů, které v sobě spojuje účinný nástroj pro design s obsáhlými grafickými analytickými nástroji. Rozšířená (stochastická) verze představuje nástroj pro vyhodnocení rizik.
3. 4 Popis programu Processing Modflow Pro
3.4.1 Historie a popis Processing Modflow Pro
Program Processing Modflow Pro (dále PMWIN Pro) byl původně vyvinutý pro sanační projekt skládky odpadů v severním Německu v roce 1989, pouze jeden rok po prvním oficiálním vydání modelu řešícího proudění podzemních 14
vod MODFLOW. První verze s vlastním prostředím běžící v systému MS-DOS, napsaly Chiang a Kinzelbach v roce 1991. Ačkoli tato verze podporovala pouze programy MODFLOW – 88 a MODPATH a grafický výstup byl omezený na hydrogramy a izolinie piezometrických výšek, reprezentoval technický průlom v rámci grafických výstupů z modelů popisujících proudění podzemní vody. S rozvojem výpočetní techniky se vyvíjel i PMWIN Pro. Poslední verze již podporuje výpočetní modul MODFLOW – 2000 a mnoho výpočetních i grafických podprogramů (MODFLOW, PEST, UCODE, MODPATH,. 3D MASTER)
3.4.2 Součásti Processing Modflow Pro
Model MODFLOW jak uvádí McDonald (1998) a Vlnas (2003) je založen na konceptu zvodní s napjatou hladinou a zvodní s volnou hladinou, to znamená, že řeší vodorovné proudění podzemní vody v jednotlivých zvodních zvlášť. Vzájemná interakce vrstev je vyjádřena vertikálním přetokem z jedné vrstvy do druhé, který je buď přímo zadán, nebo je vyčíslen z vertikálních hydraulických vodivostí sousedních vrstev. Dalšími fyzikálními a hydraulickými parametry prostředí vstupujícími do modelu jsou horizontální hydraulické vodivosti, storativita nebo pórovitost. Pro řešení základní diferenciální rovnice popisující proudění podzemní vody je užita metoda konečných diferencí (kapitola 4.1) s uzly umístěnými do středů bloků pravoúhlé sítě. Oblast proudění je jednoznačně vymezena polohou spodního a horního okraje každé zvodně a hranicí zájmové oblasti. V každém uzlu sítě je možné zadat okrajovou podmínku 1.-3. typu. Pro uživatele programu jsou připraveny tyto moduly: název studny, drenáž, evapotranspirace, infiltrace, tok bez kontroly průtoku, tok s kontrolou průtoku, obecná tlaková okrajová podmínka. Je možné si vybrat ze tří možných iteračních způsobů řešení iteračních rovnic: metoda implicitní, metoda konjugovaných gradientů a relaxační metoda. Model simuluje jak stacionární tak i tranzientní proudění podzemní vody. Základními výstupy jsou mapy izolinií hydraulických výšek a mapy izolinií snížení hydraulických výšek pro jednotlivé zvodně, při tranzientním proudění pro jednotlivé tlakové a časové úrovně. Pro daný tlakový stav, časový krok a vymezenou lokalitu lze vyhodnotit vodní bilanci. To znamená, že lze například v 15
případě, kdy podzemní voda přímo komunikuje s vodami povrchovými, zjistit kolik vody je ve vymezeném úseku recipientem drénováno, resp. kolik vody je z něj infiltrováno do podzemních vod. Tyto údaje mají velký význam při kalibraci modelu. Další součástí PMWIN Pro je program PEST (Parameter Estimation), který byl poprvé zveřejněn v roce 1994. Během procesu odhadu parametrů PEST hledá optimální hodnoty parametrů pro které je suma odchylek čtverců mezi pozorovanými a vypočítanými hodnotami minimální. Odhady parametrů jsou řízeny Gauss-Marquardt-Levenberg algoritmem. Modulární program MT3D řešící transport kontaminantů byl vytvořen Zheng v roce 1990, byl dále vyvíjen až do verze MT3DMS vytvořené Zheng and Wang v roce 1998. Simuluje změny koncentrace kontaminantu v podzemních vodách, ke kterým dochází vlivem advekce, disperze, difuse a chemických reakcí. V MT3DMS lze použít mnoho druhů okrajových podmínek, externích zdrojů a propadů. MT3DMS umí řešit transport kontaminatů jak ve zvodních s volnou hladinou, napjatou, tak i variabilní volno/napjatou hladinou podzemní vody. Dále řeší transport kontaminantů z externích zdrojů, vrtů, drénů, vodních toků a plošného znečištění. Výsledky programu MODFLOW je možné využít jako vstupní údaje pro program MODPATH, pomocí kterého je možné vypočítat trajektorie částic v dané oblasti. Program WeTbech 360 3D Master vytvořen v roce 2002 týmem autorů Wen Hsing – Chiang, Jeff
Chen, Jeff Lin. 3D Master je vizualizační nástroj, který je
schopný zobrazovat a animovat vektorová data, letecké snímky a numerické výsledky simulačních modelů. 3D zobrazování dat výrazně ulehčuje simulaci a interpretaci výsledků numerických modelů.
4. Numerické metody
Parciální diferenciální rovnice popisující trojrozměrné nestacionární proudění podzemní vody jsou v jejich obecné formě analyticky jen obtížně řešitelné a proto se v praxi používají různá zjednodušení (např. hydraulický přístup – zanedbání vertikální složky rychlosti proudění a převedení prostorového proudění na rovinné) umožňující alespoň přibližná řešení konkrétních problémů. Analytické řešení je většinou možné jen v případě, že zájmová oblast má 16
jednoduchý tvar, prostředí je homogenní izotropní, počáteční podmínka je definována konstantní hodnotou v celé oblasti a na hranicích platí jednoduché okrajové
podmínky.
Pokud
charakter
úlohy
je
v souladu
s uvedenými
zjednodušujícími předpoklady, pak lze nalézt analytické řešení i poměrně složitého problému. S rozvojem numerické matematiky a výpočetní techniky se do popředí v řešení úloh popisovaných parciálními diferenciálními rovnicemi dostaly numerické metody. Pomocí numerických metod, z nichž jsou pro řešení proudění podzemní vody nejčastěji používané metoda konečných diferencí (dále MKD) a metoda konečných prvků dále (MKP), se vytvářejí numerické modely, které umožňují na počítači simulovat požadovaný děj. Výhodou simulačních modelů je, že nevyžadují pravidelný tvar hranice řešené oblasti, prostředí nemusí být homogenní ani izotropní, na různých částech hranice mohou platit různé okrajové podmínky, uvnitř modelované oblasti se mohou vyskytovat zdroje a propady s časově proměnou hodnotou dotace či odběru apod. (Valentová, 1998)
4.1 Metoda konečných diferencí
Metoda konečných diferencí je pravděpodobně nejstarší numerickou metodou, kterou našla své využití ještě před zavedením výkonných počítačů. Princip MKD spočívá v nahrazení parciálních derivací vyskytujících se v základních
řídících
rovnicích
algebraickými
výrazy
vyjadřujícími podíl
konečných diferencí závislé a nezávislé proměnné. Jako příklad je uvedena aplikace MKD při časové a prostorové diskretizaci rovnice popisující dvojrozměrné nestacionární tlakové proudění podzemní vody v horizontální rovině v homogenní anizotropní zvodni. Uvedené proudění je popsáno rovnicí:
kx ⋅
∂2H ∂ 2 H S p ∂H + k ⋅ = y ∂ x2 ∂y 2 h ∂t
(1.1)
s okrajovými podmínkami: _
H(t) = H (t), t > 0, kx ⋅
na Γ1
∂H ∂H ⋅ nx + k y ⋅ ⋅ ny = 0 ∂x ∂y
(1.2) na Γ2
(1.3)
počáteční podmínkou 17
v Ω v čase t = 0,
H(x,y,0) = H0 (x,y) h
...
mocnost zvodnělé vrstvy [m],
H(x,y,t)
...
piezometrická výška [m],
_
H (t) H(x,y,0)
… předepsaná piezometrická výška na části hranice Γ1 [m], ... piezom. výška v čase t = 0 na náhradní oblasti Ω [m],
Kx,ky
...
koeficienty hydraulické vodivosti ve směru os x a y[m.s ],
nx,ny
...
směrové kosiny vnější normály k části hranice Γ2,
Sp
...
koeficient pružné zásobnosti (storativity),
t
...
čas [s],
x, y
...
prostorové proměnné [m],
Γ1
... část hranice oblasti se zadanou hodnotou piezometrické výšky,
Γ2
...
nepropustná hranice,
Ω
...
obast řešení.
-1
Pro přibližné řešení rovnice (1.1) se použije metoda konečných diferencí pro pravoúhlou síť se stejným krokem ∆x a ∆y. Její podstata záleží v tom, že se: •
operátory parciálních derivací nahradí diferenčními operátory;
•
úloha převede na řešení soustavy lineárních algebraických rovnic pro diskrétní body na časové ose.
Prostorové parciální derivace ve směru os x a y z rovnice (1.1) je pak možné vyjádřit následovně: ∂H i , j ∂x
≈
H i +1, j − H i −1, j
∂H i , j
,
2∆x
∂y
≈
H i +1, j − H i −1, j 2∆y
(1.5)
a ∂2H ∂x 2
≈
H i +1, j − 2 ⋅ H i , j + H i −1, j ∆x 2
∂ 2 H H i +1, j − 2 ⋅ H i , j + H i −1, j ≈ , ∂y 2 ∆y 2
(1.6)
Levou stranu rovnice (1.1) je pak možné zapsat následovně: kx ⋅
H i +1, j − 2 ⋅ H i , j + H i −1, j ∆x 2
+ ky
H i +1, j − 2 ⋅ H i , j + H i −1, j ∆y 2
.
(1.7)
Derivaci v čase na pravé straně rovnice lze zapsat s použitím konečných diferencí následovně: S ∂H i , j S H i , j (t + ∆t ) − H i , j (t ) ≈ ⋅ . h ∂t h ∆t
(1.8)
Vztah (1.8) obsahuje v uzlu (i,j) hodnoty piezometrické výšky v časech t a t+∆t . Zaveďme značení t1 = t0 + ∆t, t2 = to+2*∆t ,…, tn = to+n*n∆t, kde t0 je
18
čas odpovídající t = 0. Existuje několik možností, která z těchto časových úrovní bude pro aproximaci v prostoru. Nejjednodušší způsob je explicitní metoda. Předpokládá, že všechny hodnoty piezometrické výšky při prostorové aproximaci ve výrazu (1.7) uvažovány od začátku časového intervalu. Po dosazení pravé strany (1.7) a (1.8) do rovnice (1.1) obdržíme: t1 t0 H it−01, j − 2 ⋅ H it,0j + H it+01, j H1t,0j −1 − 2 ⋅ H it,0j + H it,0j +1 S H i, j − H i, j = kx ⋅ + k ⋅ (1.9) y h ∆t ∆x 2 ∆y 2
V rovnici (1.9) je pouze jedna neznámá hodnota piezometrické výšky
H it,1j , kterou lze snadno vyjádřit a vypočítat tak pro všechny uzly oblasti nové hodnoty piezometrické výšky v čase t1. Postup se opakuje pro časovou úroveň t2 s použitím výsledků řešení t1. Tento postup může být aplikován jak při jednorozměrné schematizaci, tak při obecně trojrozměrných úlohách. Explicitní metoda je vhodná pouze pro velmi krátké časové kroky, což pro většinu úloh hydrodynamiky není vhodné. Omezení týkající se velikosti časového kroku vede ke značnému zvýšení jejich počtu. Implicitní metoda řeší problém nestability vzniklé při použití explicitního řešení. Využívá aproximace hodnot na konci časového intervalu, nebo obecně ve vhodném mezilehlém bodě daného časového intervalu:
H i , j = ε ⋅ H it,0j + (1 − ε ) ⋅ H it,1j , ε ∈ 〈 0;1〉
(1.10)
kde ε je interpolační parametr. Pro ε = 1 vede vztah (1.1) na explicitní schéma, pro ε = 0 vede na plně implicitní schéma. Velmi často se používá tzv. Crank-Nicholsonovo schéma s hodnotou ε = 0,5. Pro plně implicitní schéma má výsledná rovnice tvar: t1 t0 H it−11, j − 2 ⋅ H it,1j + H it+11, j H1t,1j −1 − 2 ⋅ H it,1j + H it,1j +1 S H i, j − H i, j = kx ⋅ + k ⋅ (1.11) y h ∆t ∆x 2 ∆y 2
Ve vztahu (1.11) se objevuje pět neznámých hodnot piezometrické výšky v čase t1 = t0 + ∆t v pěti uzlech (i-1, j), (i, j), (i+1, j),(i, j-1), (i, j+1). Tyto hodnoty závisejí pouze na jedné známé hodnotě piezometrické výšky v čase t0. Proto je třeba v rovnici (1.11) sestavit pro všechny uzly oblasti a tím získat uzavřenou soustavu rovnic, kterou je možno řešit přímou (finitní) metodou, nebo nepřímou (iterační) metodou. (Říha, 1997) Mezi hlavní výhody této metody patří její snadná algoritmizace a nevznikají problémy při konvergenci úloh. Hlavní nevýhodou MKD je, že její aproximativní řešení nedává spojitou funkci. Další nepřesnosti vznikají při derivacích 19
diskrétních hodnot (např. určení tlakového gradientu) a nemožnosti přesně vystihnout tvar oblasti.
4.2 Metoda konečných prvků
Ve srovnání s metodou konečných diferencí představuje metoda konečných prvků (dále MKP) novější numerickou metodu, která pro svou výhodnost doznala v modelování proudění podzemní vody velkého rozšíření. Zatímco v MTD hledáme řešení pouze v izolovaných bodech (uzlech sítě). V MKP je hledaným řešením spojitý nebo po částech spojitý průběh neznámé veličiny v celé řešené oblasti, která je předem rozdělena na konečné prvky. Na rozdíl od MKD, která vyžaduje ortogonální síť, není při tvorbě sítě konečných prvků nutno dodržovat žádnou pevnou strukturu, síť je možno přizpůsobit složitým tvarům řešené oblasti a je možné ji lokálně zahušťovat. Konečné prvky mohou mít tvar obecného trojúhelníku nebo čtyřúhelníku s různým počtem uzlů (ve vrcholech i na stranách), je dokonce možné použít i prvky s křivočarými stranami. (Valentová, 1998) Výhodou MKP oproti MKD je, že každý konečný prvek může mít obecně různé fyzikální vlastnosti, které je během výpočtu možné měnit na základě získaných mezivýsledků. Výsledná matice soustavy algebraických lineárních rovnic je symetrická a pásová s dominantní diagonálou. Dále je v MKP výrazně snazší realizace okrajových podmínek (Říha, 1997) Při aplikaci metody konečných prvků se nejčastěji vychází ze dvou principů: a)Variační princip b)Princip vážených reziduí Při použití variačního principu se řešená úloha nejprve převádí na variační problém, tj. na problém nalezení funkce, která udílí extrémní hodnotu určitému funkcionálu. Metoda vážených reziduí spadá mezi tzv. přímé variační metody, které vycházejí přímo z diferenciálních rovnic popisující řešený problém. Metodu lze použít i pro problémy, pro které není klasická variační formulace známa. Při řešení těchto rovnic se používá např. Galerkinova metoda. (Valentová, 1998) 20
Hlavní nevýhody MKP
jsou složité algoritmizace popisovaných úloh a
problémy s divergencí v případě smíšených problémů.
5. Základní rovnice a děje popisující proudění podzemní vody a transportu znečištění
5.2 Darcyho zákon
Mezi základní zákony popisující proudění podzemní vody patří Darcyho zákon, který je založen na principu proudění vody ve válci naplněném pískem. Písek je zcela nasycen vodou a množství vody vtékající do válce za jednotku času se rovná množství vody z válce vytékajícímu.
Obr. 1 Darcyho pokus
Závěry mohou být vyjádřeny pomocí rovnice označené jako Darcyho zákon:
Q=
K ⋅ S ⋅ (H 1 − H 2 ) L
kde Q
(2.1) 3
-1
… proleklé množství vody [m .s ], 2
S
… průřezová plocha [m ]
L
… délka válce [m]
H1–H2
…
ztráta hydraulické výšky [m]
21
H1
…
hodnota hydraulické výšky na začátku vzorku [m]
H2
…
hodnota hydraulické výšky na začátku vzorku [m]
K
… koeficient hydraulické vodivosti [m.s ].
-1
Zobecněný Darcyho zákon zapsaný v diferenciální formě Darcyho rychlost:
v = −K
H1 − H 2 ∆H = −K ⋅ , L ∆l
(2.2)
nebo v diferenciální formě: v = −K
dH dl
(2.3)
kde dH/dl = gradient hydraulické výšky. Pro trojrozměrné proudění v homogenním prostředí, kdy je hydraulická vodivost konstantní, platí rovnice: v = − KJ = − KgradH
(2.4)
kde v je vektor rychlosti se složkami vx, vy, vz ve směru souřadnicových os x, y, z ( v x = − K
∂H ∂H ∂H , vy = − K , vz = − K ,), J je hydraulický gradient se ∂x ∂y ∂z
složkami Jx, Jy, Jz ( J x = −
∂H ∂H ∂H , Jy = − , Jz = − ). ∂x ∂y ∂z
5.3 Rovnice kontinuity a třírozměrného proudění
Rovnici třírozměrného proudění podzemní vody (3.2) vychází ze spojení zobecněného tvaru Darcyho zákona – pohybové rovnice a rovnice kontinuity ( 3.1 ). Rovnice kontinuity vyjadřuje poměr množství vody akumulované v objemu za čas ∆t. Množství vody do objemu přiteklé je rovno množství vody vyteklé za stejný čas. V nestlačitelném prostředí je S = 0, což nám vykrátí pravou stranu rovnice a dostáváme Laplaceovu rovnici (3. 3); konstanta hydraulické vodivosti je v homogenním izotropním prostředí konstantní. (Valentová, 1998)
∂ ρv x ∂ ρv y ∂ ρv z + + − ∂y ∂z ∂x
∂H = ρ ⋅ S , ∂t
( 3.1 ) 22
∂2 H ∂2 H ∂2 H ∂H =S K 2 + + 2 2 ∂y ∂z ∂t ∂x
( 3.2)
∂2 H ∂2 H ∂2 H + + =0 ∂x2 ∂ y2 ∂z2
( 3.3)
S
… specifická storativita
5.4 Šíření znečištění
Některé látky vstupující do pórovitého prostředí jsou s vodou nemísitelné, nebo jen částečně mísitelné (např. nafta neb oleje). Souběžný pohyb vody a těchto látek se označuje jako vícefázové proudění. Pokud se naopak jedná o rozpustné látky, které se mísí s vodou, nazývá se jejich pohyb v půdních a horninových vrstvách šíření rozpuštěných látek. Složitý komplex procesů ovlivňujících pohyb rozpuštěných látek pod povrchem je možno rozdělit na ty, které jsou určující pro pohyb látek a ty, které způsobují jejich transformace. V této souvislosti mluvíme u látek a látkové přeměně. Pokud dochází pouze k šíření nereagujících nebo málo reagujících látek, jedná se o konzervativní transport. O látkovou přeměnu jde v případech, kdy jsou chemické reakce rozpuštěných látek nezanedbatelné, nebo dokonce dominantní, a je nutné je do popisu zahrnout. V tomto případě se jedná o nekonzervativní transport. Matematický popis transportních procesů vychází stejně jako v případě proudění fází z makroskopického popisu založeného na mechanice kontinua. Soubor tekutých fází je nahrazen sérií fiktivních fázových makrokontinuí, jednotlivé rozpuštěné látky (složky směsí) jsou nahrazeny složkovými makrokontinui. (Vogel, Císlerová, 1998)
23
5.4.1 Rovnice popisující šíření znečištění
a) Konzervativní transport Jednorozměrná migrace nesorbující se látky v 1-D filtračním poli je definována obecnou diferenciální rovnicí ( 4.1).
Dx
∂ 2C ∂C ∂C − vxp = 2 ∂x ∂t ∂x
I
II
(4.1)
III
Jednotlivé členy této rovnice vyjadřují : I. Disperzivní transport, II. Advektivní transport, III.
Nestacionaritu procesu v čase
b) Nekonzervativní transport 3-D migrace sorbující a reagující látky v 3 - D filtračním poli je definována obecnou diferenciální rovnicí (4.2). ∂ (vx C) ∂ (v y C) ∂ (v C) ∂C ∂2 C ∂2 C ∂2 C z = n D + + n o + n D + n D ± o hx o hy o hz 2 2 2 ∂ t ∂ x ∂ y ∂ z ∂x ∂y ∂z
I
II
∂ Cs ± no ∂t
IV
±
III
∂C ± Qi ⋅ C + − Qi C ± no ∂ t
(4.2)
VI
V no
... efektivní pórovitost
C
... koncentrace rozpuštěné látky [mg.l ]
Dh
... koeficient hydrodynamické disperze [m .s ]
v
... rychlost [m.s-1 ]
Cs
…
koncentrace látky sorbovaná na povrchu pevné fáze[mg.l ]
Qi
...
průtok [m .s ]
-1
2
-1
-1
3
-1
I) Změna koncentrace rozpuštěné látky v bodě [x,y,z]. II) Rychlost šíření látky z nebo do bodu [x,y,z] vlivem advekce. III) Rychlost šíření látky z nebo do bodu [x,y,z] vlivem disperze. IV) Rychlost sorbce rozpuštěnné látky v bodě [x,y,z]. V) Rychlost změn koncentrace rozpuštěné látky vyvolané reakcemi. VI) Přítok / odběr rozpuštěné látky z vnějšího zdroje / propadu. (www.fle.czu.cz)
24
5. 5 Okrajové a počáteční podmínky
Při řešení hydrogeologických úloh šíření znečištění, jak uvádí Valentová (1998) rozlišujeme následující okrajové a počáteční podmínky:
5.5.1 Hranice s předepsanou hodnotou hydraulické výšky
Jestliže ve všech bodech hranice řešené oblasti nebo na její části známe hodnotu hydraulické výšky po celou dobu zkoumaného procesu, jedná se o hranici s předepsanou hodnotou hydraulické výšky – okrajová podmínka prvního typu, nazývaná také Dirichletova. Tuto podmínku lze vyjádřit pomocí zápisu: H = f1 (x,y,z) nebo H = f2 (x,y,z,t) na S ,
(5.1)
kde f1 a f2 jsou známe funkce. První případ vyjadřuje stacionární podmínku, zatímco ve druhém případě je okrajová podmínka závislá na čase. Okrajové podmínky tohoto typu se vyskytují vždy tam, kde je oblast proudění ve styku s otevřenou volnou hladinou – řekou, jezerem apod.
5.5.2 Hranice s předepsaným tokem
Jestliže ve všech bodech hranice je známa hodnota toku ve směru kolmém na hranici, jedná se o hranici s předepsaným tokem vn = f (x,y,z,t) na S ,
(5.2)
kde vn je složka rychlosti kolmá k hranici oblasti a f (x, y, z, t) je známá funkce. Tato okrajová podmínka se nazývá také jako okrajová podmínka druhého typu, nebo Neumannova okrajová podmínka. Speciálním případem této okrajové podmínky je nepropustná hranice. V případě nepropustné hranice je složka hustoty toku kolmá k hranici rovna nule vn = 0. 5.5.3 Polopropustná hranice
Tento typ okrajové podmínky se vyskytuje tam, kde oblast proudění je v kontaktu s otevřeným vodním zdrojem (nebo jiným porézním prostředím), ale je oddělena polopropustnou vrstvou. Pro tento typ okrajové podmínky se používá 25
označení třetího typu nebo smíšená okrajová podmínka, Cauchyho okrajová podmínka, nebo okrajová podmínka třetího typu. Předpokládejme, že H je hodnota hydraulické výšky uvnitř řešené oblasti a H0 hodnota hydraulické výšky vně oblasti. Velikost toku kolmo k hranici řešené oblasti lze vyjádřit pomocí vztahu : vn = (H0 – H ) / c, kde c = B / K,
(5.3)
kde B je šířka polopropustné vrstvy, K je její hydraulická vodivost, c je odpor vrstvy.
5.5.3 Počáteční podmínky
Počátečními podmínkami rozumíme stav charakterizující proudění v čase t0 = 0. Stanovení těchto podmínek nám umožňuje řešit nestacionární proudění. Schematicky lze podmínku vyjádřit jako funkci f () souřadnic x, y, z v prostoru v čase t0, kdy známe hydraulickou výšku H. Průběh hydraulické výšky se s časem mění. (Mucha, Šestakov, 1987) H = f ( x, y, z , t 0 )
(5.4)
6. Testování použitelnosti programu. Lokalita Kralupy
V rámci ověření použitelnosti programu PMWIN pro jsem se zaměřil na využití metody matematického modelování na podmínky firmy Kaučuk Group, kde jsou dlouhodobě vytvářeny hydrodynamické podmínky zamezující šíření znečištění. Jako hlavní ochranná metoda je použito čerpání podzemních vod na tzv. hydraulické bariéře (dále HOPV). Firma Kaučuk patří k nejvýznamnějším petrochemickým podnikům v České republice.
6. 1 Ekologie petrochemických podniků
Petrochemie patří k jednomu z hlavních zdrojů znečištění životního prostředí. K úniku nežádoucích organických sloučenin dochází při jejich přepravě, manipulaci a zpracování a následně při nakládání s odpady a odpadními vodami. Do skupiny tzv. ropných produktů jsou zařazovány a) nafta, b) benzin c) mazací oleje, d) topné oleje, e) speciální chemické výrobky a meziprodukty, 26
které se používají v chemickém průmyslu. Ke znečišťujícím látkám patří: a) uhlovodíky (nasycené i nenasycené, alifatické i aromatické), b)organické sloučeniny síry (thioalkoholy, thiofenoly, sulfidy), c)organické sloučeniny kyslíku (alkoholy, fenoly, kresoly, furfurol, aldehydy, ketony, naftenové kyseliny, estery), d) organické sloučeniny dusíku (aminy,
deriváty
pyridinu),
e)
chlorované
uhlovodíky
(dichlorethan,
etylendichlorid), f) další látky (oleje, tenzory, oxidy kovu, NAOH, H2S, NH+4. Při řešení znečištění v areálu Česká Rafinérská A. S., Kralupy nad Vltavou, se budu zabývat hlavně znečištěním nesorbujících se látek. Protože se v rafineriích spotřebovává obrovské množství vody, leží tyto závody zpravidla v blízkosti velkých řek (Obr. 2). Díky velkému množství spotřebované vody, vzniká velké množství vod odpadních, což může mít za následek řadu ekologických problémů. Proto většina rafinérských podniků má ve svém areálu vlastní čistírnu odpadních vod včetně laboratoří na jejich rozbor.
Obr.2 Areál Kaučuk Group a.s.,Kralupy nad Vltavou (www.mapy.cz)
27
6.2 Podmínky provozu podniku
Výroba syntetického styren-butadienového kaučuku (SBR) byla zahájena v roce 1963. Na základě monomerů butadienu a styrenu, které se oba připravují také v podniku, se vyrábí i druhý základní výrobek – polystyrenové plasty (PS) v různých variantách jako houževnatý, krystalový, zpěňovatelný PS a kopolymer ABS. Součástí podniku je vlastní energetika, ve které se vyrábí potřebná pára a elektrický proud zejména pro vlastní spotřebu, ale i pro prodej. V roce 1975 byl podnik rozšířen o rafinerii, ve které se na palivářské výrobky (topné plyny, benziny, nafta, petrolej a topný olej) zpracovává ropa. Od 1. 1. 1996 byla rafinerie v rámci restruktualizace petrochemického a rafinérského průmyslu oddělena a začlenila se do nové společnosti ČESKÁ RAFINÉRSKÁ, a.s., Od 1. 7. 1997 je KAUČUK, a.s., součástí skupiny UNIPETROL, a.s.
(Zpráva o vlivu, 2000) 6. 3 Ekologické a přírodní poměry areálu Kaučuk Group a.s., Kralupy nad Vltavou a v širším okolí
6. 3. 1 Vymezení areálu Kaučuk Group a.s.
Zájmový areál Kaučuku se nachází 20 km severně od Prahy, na pravém břehu řeky Vltavy, která je od rafinérie vzdálena 0,8 km. Areál tvoří severovýchodní okraj města Kralupy nad Vltavou. Při svém severním okraji vlastní areál téměř navazuje na další město – Veltrusy. Pozemek Kaučuku je protáhlého tvaru severojižní orientace s rozměry cca 2,0 x 1,0 km a rozlohou 1,62 km2. Ze správního hlediska se základní areál Kaučuk Group a.s. nachází v jihozápadním výběžku okresu Mělník na několika katastrálních územích: •
k.ú. Lobeček
•
k.ú. Chvatěruby
•
k.ú. Veltrusy
Přehled měst a obcí v blízkém okolí areálu Kaučuk Group a.s. : •
Kralupy nad Vltavou, navazují na areál směrem na jihozápad
•
Nelahozeves, 1,5 km směrem na západ (dále Z)
•
Hleďsebe, 1,3 km směrem na severozápad (dále SZ) 28
•
Veltrusy, 0,5 km směrem na sever (dále S)
•
Úžice, 3,0 km směrem na východ (dále V)
•
Kozonoh, 2,7 km směrem na jihovýchod (dále JV)
•
Chvatěruby, 1,5 km směrem na jihojihovýchod (dále JJV)
Přehled Území se zvláštní ochranou: •
NPR Kopeč – stepní porosty vzácné Lipnice Bádenské, rozloha 2 ha nachází se cca 6,5 km směrem na V
•
PP Otrokovická skála – stepní společenstva na algonkických břidlicích, nachází cca 5 km na JV
•
PP Minická skála – krajinná dominanta, skála s teplomilnou květenou, rozloha 0,37 ha, nachází se cca 3,5 km na jihozápad (dále JZ)
•
PP Sprašovská rokle u Zeměch – teplomilná květena, rozloha 1,49 ha, nachází se cca 4 km směrem na JZ
•
PP Netřeská slaniska – zbytek slanomilných společenstev, rozloha 1,01ha, nachází se cca 7 km směrem na V (Závěrečná zpráva, 1997)
Obr.3 Modelované území (www.mapy.cz)
29
6. 3. 2 Topologické poměry
Vlastní lokalita areálu se nachází na okraji rovné údolní terasy a zaujímá i výrazný terénní stupeň vysoký cca 12 m na východní hranici areálu. Směrem k západu od tohoto stupně je terén mírně zvlněný. Nadmořská výška areálu je 175 až 176 m n.m., oblast nad terénním stupněm pak 187 m n.m. Samotný terénní stupeň má úklon 15-19°k západu. (Záv ěrečná zpráva, 1997) Od hranic areálu až k řece se terén mírně svažuje až na 170 m n.m. Významné je, že celá oblast je zjednodušeně ze tří stran omezena řekou Vltavou a z východní pak sedimenty křídového a permokarbonského stáří.
6.3.3 Klimatické a hydrologické poměry
Z hlediska klimatických charakteristik se jedná o teplou oblast, suchou s mírnou zimou. Průměrná roční teplota je 8 – 9 °C. Pr ůměrné roční srážky na základě pozorování ve srážkoměrných stanicích Kralupy nad Vltavou a Všetudy z období 1931-1960 jsou 477-493 mm. Zájmová oblast spadá do povodí řeky Vltavy s číslem hydrologického pořadí 1-12-02. Průměrný průtok z údajů limnigrafické stanice ve Vesňanech (cca 7 km po proudu od Kralup je 149,91 m3.s-1 a Q355 je 26,39 m3.s-1. Do Vltavy se v zájmovém území vlévá Mlýnský potok.
Srážkoměrá Stanice Kralupy nad Vltavou
I
II
III
IV
V
VI
23
22
23
31
55
67
měsíc VII 75
VIII
IX
X
XI
XII
59
35
37
25
25
Celkem Za rok [mm] 477
Tab.1 Průměrné měsíční úhrny srážek (1931 – 1960)
6.3.4 Geomorfologické a geologické poměry
Zájmové území je zařazeno do Středolabské tabule a jejího podcelku Mělnické kotliny. Jedná se o erozně denudační sníženinu v širší oblasti soutoku Vltavy s Labem a při dolním toku Vltavy. Ploché dno je charakterizováno akumulačním reliéfem údolních niv, mladopleistocenních a středopleistocenních
30
říčních teras. Vzácněji se vyskytuje denudační reliéf zarovnaného slínovcového povrchu – kryopediment. Na řešeném území byla vrtnými pracemi prokázána přítomnost písků a štěrkových písků pod terénním skokem. Ty náleží terase spodního wurmu s mocností,
která
se
plynule
ubývá
ve
směru
od
severozápadního
k jihozápadnímu (směrem k terénnímu skoku) s maximem 15 m a vykliňuje na bázi terénního skoku. Tyto písčité sedimenty jsou místně překryty tenkou vrstvou půdy nebo navážek nicméně často vycházejí až napovrch. Písčité nebo štěrko-písčité sedimenty, které jsou relikty riiské terasy, se nacházejí v území nad terénním skokem. Jejich mocnost je podstatně nižší a rychle ubývají směrem od terénního skoku. (Závěrečná zpráva, 1997)
6.3.5 Proudění podzemních vod
Proudění podzemních vod je ovlivněno tokem Vltavy, kdy v přirozeném stavu je směr proudění souhlasný s řekou, tj. cca k severu. (Závěrečná zpráva, 1997) V období před provozem hydraulické ochrany podzemních vod (dále HOPV) byl celý proud podzemní vody ve směru k Veltrusům, tvořený infiltrací z nadjezí Miřejovického jezu do údolní nivy a přítokem z vyšších teras, odvodňováno do úseku Vltavy v podjezí. Po zapojení HOPV, a to již od samého počátku, je tento proud ochuzen o čerpání vody z jejích objektů (minimálně o cca 40 l.s-1, maximálně až o více než 80 l.s-1) což má za následek ovlivnění úrovní hladin vody i ve studnách ve Veltrusech. Občany je ovlivnění registrováno vždy až v souběhu se sucho periodou (Kliner, 2000). Vrtný průzkum provedený v cenomanských sedimentech prokázal, že hladina podzemní vody je volná a že nedochází k proudění směrem do centra křídové pánve, ale naopak směrem k Vltavě a čerpané hydraulické cloně. V dané oblasti se nachází jedna hlavní aluviální zvodeň, další zvodně se nacházejí ve větších hloubkách a nemají vliv na aluviální zvodeň a ani nejsou zvodní ovlivňovány. Výjimku tvoří kontakt aluviální struktury na východním okraji oblasti, kde dochází s velkou pravděpodobností k příronu. Mocnost
zvodnění terasy je
do
značné
míry ovlivněna
morfologií
karbonského podloží terasy. Vyskytují se mocnosti od 0 m (pod terénním stupněm) až do 10 m. Při východním okraji rafinerie se nachází výrazná elevace karbonských sedimentů, která způsobuje výrazné snížení mocnosti kvartérní 31
zvodně. V několika monitorovacích vrtech hladina podzemní vody nedosáhla báze kvartérních sedimentů a zůstala zakleslá do permokarbonských hornin. (Závěrečná zpráva, 1997)
6.3.7 Využití vodních zdrojů
Povrchové vody jsou využívány ve velkém rozsahu jako chladící vody. Ty jsou jako surovina odebírány z toku Vltavy proti směru proudění. Odběrné místo se nachází při vtoku Vltavy do Kralup nad Vltavou ve vodárně u potrubního mostu. Použité vody (odpadní) jsou zpětně vypouštěny otevřeným kanálem zpět do Vltavy poblíž Městské čistírny odpadních vod. Podzemní vody se v areálu Kaučuk Group.a.s. nevyužívají. Výjimkou je čerpání podzemních vod na HOPV, které však nelze považovat za využívání podzemních vod, ale za jejich ochranu. K jímání a odběrům podzemní vody ovšem dochází v domovních studnách v Kralupech nad Vltavou a ve Veltrusech. (Závěrečná zpráva, 1997)
6.3.8 Zdroje a druhy kontaminace
V celém areálu Kaučuk lze očekávat kontaminaci prostředí následujícími polutanty: ropnými látkami (NEL), polyaromatickými uhlovodíky (PAU), skupinou BTEX (benzen, toluen, etylbenzen, xylen), styrenem, methyltercbutyletherem (MTBE), chloridy, těžkými kovy a polychlorovanými bifenyly (PCB). Nejvíce znečištěné plochy byly identifikovány v oblasti uvnitř kolejiště a západně od něj, kde jsou umístěny stáčecí a plnící plochy. V těchto místech byly zjištěny vysoké koncentrace NEL v zeminách a to jak v přípovrchových partiích nesaturované zóny, tak i v hloubkovém intervalu 0,5 m nad hladinou podzemní vody. Navíc byla v této oblasti zjištěna volná fáze ropných uhlovodíků (dále NAPL) o maximální mocnosti 1,17 m. Oblast silničního distribučního střediska je zdrojem kontaminace podzemních vod BTEX. Mezi prokázané zdroje znečištění patří : a)nádržové dvory a záchytné jímky, které byly v minulosti zpevněny pouze asfaltem, b)čerpací stanice topných olejů, c) etylazice benzínů, d) netěsnící kanalizace na Silničním distribučním středisku, e) celý areál rafinerie – zaolejovaná kanalizace. (Závěrečná zpráva, 1997) 32
7. Modelové řešení
Při modelovém řešení byl použit již zmíněný program PMWIN Pro verze 7.0.0 z roku 2002. Dále byla využita data za období 20005 – 2006. Řešení bylo rozděleno do následujících kroků •
Informační zabezpečení modelu a vyhodnocení dat
•
Vymezení modelové oblasti
•
Sestavení rámcového modelu
•
Stanovení okrajových a počátečních podmínek
•
Diskretizace parametrů prostředí
•
Modelování stávajícího stavu
•
Kalibrace modelu s ohledem na stupeň znalostí hydrogeologických podmínek (vliv kolmatace řeky, vliv nerovnosti podloží, vliv režimu napjatosti hladiny, vliv přetékání ze spodních a bočních rozhraní, vliv změny režimu odběru podzemních vod, vliv vzdutí hladiny na jezu, vliv infiltrace, vliv hydraulických odporů na vrtu) a optimalizace podle konkrétních výsledků monitorování hladin (PEST)
•
Sestavení modelu šíření znečištění
•
Provedení výpočtů a zhodnocení výsledků při zadání různých čerpaných množství podzemní vody a jejich vlivu na šíření znečištění a výšku hladiny podzemní vody v domovních studnách ve Veltrusech
7.1 Informační zabezpečení modelu a vyhodnocení dat
Pro sestavení detailního modelu byly k dispozici data z řady vrtů jejíž plošné rozmístění je nerovnoměrné. Používaná data: •
výsledky režimního čerpání z objektů HOPV za rok 2005
•
měsíční stavy hladin podzemní vody v pozorovacích vrtech
•
souřadnice umístění objektů HOPV a pozorovacích vrtů
•
geologické profily v areálu rafinerie
•
mapy izolinií znečištění podzemních vod
Údaje o filtračních parametrech jsem převzal z dostupných zpráv (Závěrečná zpráva, 1997, Kliner, 2005.), přičemž pro jejich prostorové vyjádření jsem byl nucen
metodou
analogie
orientačně
posoudit
filtrační
parametry
dle 33
petrografického popisu. V případě vlastní kalibraci jsem pak řešil sérii úloh, kdy jsem měnil tyto parametry a kalibroval model z hlediska shodnosti skutečných a vypočítaných hladin podzemní vody (dále HPV).
7.2 Problémy při sestavování modelu
Přírodní procesy jsou ve své podstatě neopakovatelné. To je způsobeno vzájemným spolupůsobením příčinných (deterministických) a nahodilých (stochastických) faktorů. Z tohoto důvodu každý pokus o jejich modelování předpokládá jisté zjednodušení. (Hrádek, 1988) Na přesnost hydrodynamického modelu mají proto vliv nejen hodnoty zadávaných parametrů, ale zejména podmínky prostředí, pro které se daný model sestavuje. Přesnost mnou sestaveného modelu je významně ovlivněna heterogenitou prostředí a to jak v horizontálním, tak vertikálním směru. Dále má na modelované území vliv blízkost velké řeky, která může působit dotačně (zejména při povodňových vlnách), nebo naopak do řeky může podzemní voda přitékat z horních geologických teras. Při povodňových vlnách významně stoupá hladina
podzemní
vody
a
dochází
k vyplavování
znečišťujících
látek
z nesaturované zóny a to má za následek několikanásobné zvyšování koncentrace v podzemních vodách a následné možnosti proniknutí znečištění přes HOPV do okolního prostředí. Přestože se na daném území trvale čerpají podzemní vody a to již od sedmdesátých let minulého století, neměl jsem k dispozici prakticky žádná data, ze kterých bych mohl sestavit a nakalibrovat model ustáleného proudění podzemních vod a ten použít jako základ pro sestavování modelu neustáleného proudění.
Ještě složitější jsou podmínky týkající se znečištění. Na
modelovaném území dochází k únikům znečištění prakticky od začátku provozu rafinerie až do dnes. Zaměřil jsem se hlavně na sestavení modelu proudění podzemní vody a znečištění jsem řešil pouze z kvalitativního hlediska.
34
7.3 Vymezení modelovaného území
Nejprve jsem si na základě znalosti území a letecké mapy stanovil rozsah modelované oblasti. Jak je patrné z obr. 1 území je z leva ohraničeno Řekou Vltava a z pravé strany terénní elevací, která je považována za nepropustnou. Výsledné území má velikost 2650 x 6000 m.
Detailní vymezení pak bylo
upřesněno podle průběhu hydraulických a hydrogeologických hranic.
7.4 Sestavení rámcového modelu
První krok při vytváření matematického modelu, používající při počítání metodu konečných diferencí (kapitola 4.1) a metod z této metody odvozených je diskretizace výpočetní sítě. Ta se v modelu PMWIN Pro definuje v dialogu grid Mesch size (obr. 4) . Výpočetní síť modelové oblasti je tvořena buňkami ve třech směrech: ve směru osy x tvoří řádky, ve směru osy y sloupce a ve směru osy z vrstvy. Protože modelované území je poměrně rozsáhlé a počet modelových buněk je omezen pouze na 9000, musel jsem volit kompromis mezi hustotou buněk a počtem vrstev. Proto jsem vytvořil síť o velikosti buňky 50 x 250 m a tu zahušťoval až na výslednou velikost jedné buňky 50 x 50 m v oblasti rafinerie a 50 x 100 m na krajích počítaného území. Výsledný model je diskretizován dvěmi modelovými vrstvami pomocí 4452 výpočetních buněk v každé vrstvě.
Obr 4. Příklad dialogového okna pro definování výpočetní sítě modelu
V dalším kroku jsem definoval šířku obou vrstev. K tomuto kroku jsem využil modulu Digitizer, který je součástí programu PMWIN Pro, ale lze jej používat i 35
samostatně. V digitizéru jsem otevřel nadefinovanou výpočetní síť , kde jsem již měl známé x a y souřadnice a na bodu se známou nadmořskou výškou definoval hodnoty z souřadnice od zvolené srovnávací roviny, která byla zvolena na 150 metrů nad mořem (dále m.n.m.) . Hodnoty souřadnic jsem uložil do souboru digi_TOP1.xyz. Tento soubor jsem otevřel v dalším podporovaném modulu Field interpolator , který pomocí metody krieging interpoluje zadané hodnoty do všech buněk výpočetní sítě. To samé jsem udělal i pro druhou vrstvu, kde povrch druhé vrstvy tvoří rozhraní geologických vrstev. Jak je vidět z obr. 5 v oblasti rafinerie je výrazná terénní elevace, která má vliv na proudění podzemních vod. Model je rozdělený na dvě vrstvy přičemž první vrstva je definována jako vrstva s volnou hladinou podzemní vody a druhá jako vrstva s hladinou volno/napjatou.
Obr 5. Modelované území (obrázek je ukázkou jednoho z výstupu modelování)
36
7.5 Stanovení okrajových a počátečních podmínek
PMWIN Pro nabízí definování mnoha druhů okrajových podmínek. V modelu jsou použity okrajové podmínky I. A III. Typu tj. konstantní hladina na řece resp. vydatnost závislá na změně hladiny v řece (vliv kolmatace), konstantní hladina a přítok z vyšších geologických vrstev. 3 Okrajová podmínka River je definována pomocí nástroje zadávání liniových hranic pomocí funkce „Polyline“, u kterého se přiřazují hodnoty dvěma konečným vrcholům a do zbývajících jsou interpolovány. Zadané hodnoty jak udává Kliner (2006) jsou zobrazeny v tabulce 2. parametr Hydraulická vodivost dna Výška hladiny v řece Elevace dna řeky Šířka řeky Šířka dna řeky
hodnota 0.0001 18, 5 od jižní hranice po jez 16, 5 od jezu po severní hranici 3 80 -100 2
jednotky -1 m.s m m m m
Tab. 2 Definované parametry řeky
Protože z průzkumných prací provedených na dané lokalitě vyplývá, že hladina podzemní vody je významně ovlivněna přítokem z vyšších teras je v modelu použita okrajová podmínka třetího typu (Kap 5.5.) GBH. Hodnoty zadávané do modelu jsou uvedeny v tabulce 3. parametr Ekvivalentní hydraulická vodivost Piezometrická výška v externím zdroji
hodnota 0.0001 23
jednotky -1 m.s m
Tab. 3 Definované parametry GBH
V modelu PMWIN Pro se buňkám zadávají hodnoty : -1
pro buňky, kde
zůstává stálá hodnota zadávaného parametru (např. piezometrická výška, koncentrace), 0 pro neaktivní buňky (např. nepropustné podloží), 1 pro aktivní buňky, v nichž model počítá požadované hodnoty.
V modelu jsou zadané
v celém řešeném území 1 na severní hranici v místě terénního zlomu je zadaná hodnota 0, protože jeho působení simuluje okrajová podmínka GBH.
3
V programu použit termín River a general head boundery (GBH).
37
Aby mohl PMWIN Pro začít řešit simulaci proudění, požaduje pro model počáteční „odhady“ hodnot hladiny (počáteční podmínky). Dobrý počáteční odhad může významným způsobem zkrátit čas výpočtu. Počáteční hodnoty hladin by měly být vyšší než je nadmořská výška dna počítané vrstvy. Proto jsem jako počáteční hodnotu hladiny nastavil 22 m od srovnávací roviny 150 m.n.m.
7.6 Diskretizace parametrů prostředí
Protože jsem neměl k dispozici žádné výsledky hydrodynamických zkoušek, určil jsem hodnoty parametrů (Tab. 4) dle převládajících druhů zemin a tyto hodnoty dále kalibroval. parametr Horizontální hydraulická vodivost Vertikální hydraulická vodivost Koeficient Storativity Volná storativita Efektivní pórovitost
hodnota 1. vrstva
2. vrstva
0.0001
0.00000001
m.s
-1
0.00001 0.001 0.025 0.1
0.000000001 0.000001 0.0000025 0.001
m.s -
-1
jednotky
Tab. 4 Parametry prostředí
7.7 Modelování stávajícího stavu
Po definici okrajových a počátečních podmínek, parametrů prostředí a časového kroku výpočtu jsem spustil výpočet úrovní hladin podzemní vody neovlivněné HOPV. Jak je vidět z Obr. 6 izolinie hladiny podzemní vody jsou ovlivněny zejména tvarem a velikostí řeky, Miřejovickým jezem a přítokem vody na severní hranici území. Tyto hodnoty piezometrických výšek jsem uložil do souboru formátu ASCII. TXT a použil jako počáteční hodnoty piezometrických výšek při modelování hydrodynamiky ovlivněné HOPV. HOPV při modelování různých variant velikostí čerpaných množství podzemní vody z jejích objektů rozděluji na kralupskou (jižní) a veltruskou (severní) větev HOPV. Pro modelování situace ovlivněné provozem HOPV jsem využil data poskytnuté firmou KAREL KLINER VODNÍ ZDROJE a.s. a do objektů HOPV zadal čerpaná množství podzemní vody (Příloha 4) vztahující se k 19. 1. 2005. a do pozorovacích vrtů zadal hodnoty naměřených hladin podzemní vody pro 38
stejné datum (Tab. 5). K Porovnání výsledků modelu a pozorovaných hodnot jsem použil nástroj Head Scatter Diagram .
Graf 1. Porovnání pozorovaných a vypočítaných hodnot v při zadání odhadovaných hodnot parametrů.
Obr. 6 Izolinie hladiny podzemní vody
Jak dokládá Graf 1., naměřené a vypočítané hodnoty se významně liší, což je způsobeno tím, že hodnoty parametrů jsou pouze odhadované.
7.8 Kalibrace a optimalizace modelu
Prvním kalibrovaným parametrem byla hydraulická vodivost dna řeky (termín též filtrační odpor řeky). Její hodnotu jsem postupně měnil a následně porovnával vypočítané a pozorované hodnoty dokud se významně nelišily. Vliv změny parametru na průběh hydroizohyps je patrný v příloze 2. Nejlepší odhadovaná hodnota hydraulické vodivosit dna řeky vyšla Cbot = 0.000001 m.s-1.
39
-1
Graf 2. Porovnání pozorovaných a vypočítaných hodnot při Cbot = 0.000001 m.s
Jak je patrné z Grafu 2. a Tab. 5 vypočítané a pozorované hodnoty se již významně neliší.
-1
Tab 5. Pozorované a vypočítané hodnoty při Cbot = 0.000001 m.s
Pro upřesnění všech parametrů jsem použil program PEST (Kapitola 3.4.2). Jako odhadované parametry zvolil Horizontální hydraulickou vodivost - HK_1, Koeficient Storativity - S, Volnou storativitu – Sy a Ekvivalentní hydraulická vodivost u okrajové podmínky typu GBH – CE. Hodnoty z první kalibrace jsem 40
dosadil
do
modelu,
spustil
výpočet
a
porovnal
pozorované
hodnoty
s naměřenými. Protože byl ještě patrný rozdíl mezi hodnotami piezometrických výšek, parametry jsem mírně upravil a spustil výpočet znovu. Tento postup jsem opakoval, dokud se piezometrické výšky významně nelišily. Výsledné hodnoty optimalizovaných parametrů po kalibraci jsou uvedeny v Tab. 6 parametr hk_1 ghb_1 sy_1 s_1
hodnota 8.94E-04 1.78E-05 2.50E-02 2.68E-04
jednotky 1 m.s1 m.s-
Tab. 6 Odhadované parametry pomocí PEST
Jak je dokládá Graf 3. a Tab 7. k odchylce dochází pouze u pozorovacího vrtu HJ-45 a HJ-10, to lze vysvětlit blízkostí okrajové podmínky GBH, ale hlavně heterogenitou prostředí, protože oba vrty jsou v místě terénní elevace. Korelační koeficient pozorovaných a vypočítaných hodnot vychází 0.9967 (příloha 8).
Pozorovací vypočítaná Pozorovaná vrt hodnota hodnota SP-8 17.06122 16.96
Graf 3 . Porovnání vypočítaných a pozorovaných hodnot při zadání nakalibrovaných parametrů
SP-9
17.12238
16.96
SP-5
16.50111
16.62
SP-6
16.59772
16.71
SP-10
16.30181
16.38
SP-11
16.47163
16.46
SP-12
16.15496
16.21
SP-7
16.29899
16.33
SP-13
16.86201
15.76
SP-26
16.58594
15.78
SP-2A
16.0033
15.61
SP-19
18.79106
19.01
SP-29
16.81107
16.5
25
16.90886
15.84
HJ-19
18.79106
18.7
HJ-8
17.57983
16.53
HJ-45
20.65031
21.89
H5
17.1445
17.29
HJ-10
19.77612
17.7
Tab. 7 Vypočítané a naměřené hodnoty
41
Pro ověření modelu (validaci) jsem do čerpaných objektů HOPV zadal hodnoty vztahující se k 17. 2. 2006 (Tab. 5) a spustil výpočet. Jak je patrné z grafu 4. spočítané hodnoty se statisticky neliší a model lze považovat za nakalibrovaný.
Graf 4. Porovnání pozorovaných a vypočítaných hodnot k 17. 2. 2006
Obr. 7 Hladina podzemní vody ovlivněná HOPV
42
7.9 Modelování šíření znečištění
7.9.1 Popis Modelu
K modelování migrace kontaminantů v areálu jsem použil program MT3DMS (kapitola 3.4.2). Z hydrochemického hlediska jde o rovnovážnou kinetiku mezi znečištěním a horninovým prostředím, neboť zde dochází k úniku nežádoucích látek minimálně od roku 1970. Skutečné podmínky šíření znečištění a vznik kontaminačních
mraků
jednotlivých
látek
jsou
stanovovány
z výsledků
monitorování. Přes ekonomický a ekologický význam dané lokality považuji stupeň a hlavně rovnoměrnost získávání dat za nedostačující. Proto jsem pro relativně malou prozkoumanost byl nucen extrapolovat známé informace do okrajových oblastí, kde je těchto informací nedostatek. Z tohoto důvodu jsem hodnotil pouze kvantitativní šíření znečištění v areálu a vliv HOPV na jeho šíření.
7.9.2 Definované vlastnosti modelu
Mezi
základní
charakteristiky
prostředí
a
kontaminantů
v programu
MTR3DMS patří advekce, disperze, chemické reakce, počáteční koncentrace příslušné látky, kinetické reakce a definice okrajových podmínek. Já jsem v modelu zadával tyto charakteristiky: a) Advekce: tento modul nabízí mnoho druhů schémat řešící advektivní složku šíření kontaminantů. Každé schéma se používá pro různé druhy kontaminantů a prostředí. V modelu používám k vypočtení advekce kontaminantů metodu konečných diferencí se zadaným courantovým číslem 0.75, které udává jakou část buňky překročí znečištění v jednom časovém kroku. b) Disperze: Do modulu disperze jsem zadal dle převládajících druhů zemin hodnotu podélné disperze 10 m v první vrstvě a 1m v druhé vrstvě. Dále jsem přiřadil hodnoty horizontální disperzivity 0.1 [-], vertikální disperzivity 0.01 [-] a difúzního koeficientu 0 [m2/den] v první vrstvě a v druhé hodnoty o řád menší.
43
c) Sorbce: Protože jsem modeloval pouze účinnost HOPV , která není závislá na rychlosti migrace kontaminantů, sorbci jsem do modelu nezadával. e) Okrajové podmínky: Protože areál rafinerie působí jako kontinuální zdroj znečištění, do míst s největší pravděpodobnosti úniku znečištění jsem zadal okrajovou podmínku prvního druhu vyjadřující stálou koncentraci kontaminantů. f) Počáteční koncentrace: Počáteční koncentraci 100 mol.m-3 jsem zadal do buněk, kam byla zadána okrajová podmínka.
7.9.3 Výsledky modelu MT3DMS Výstupem z modelu jsou matice hodnot koncentrace kontaminantů v daném místě pro určitý časový krok. Tato data se mohou použít v modulu 2D visualization, který vykreslí izolinie buněk se stejnou koncentrací v daném čase.
44
0br. 8 Znečištění po ročním kontinuálním úniku
8. Modelované situace
V této kapitole analyzuji otázku šíření znečištění a změny hladin podzemních vod v domovních studnách ve Veltrusech a to v závislosti na intenzitě odběru vod z hydraulické clony (HPOV) v areálu Kaučuk Group a.s., Kralupy nad Vltavou. Výsledky řešení několika variant modelů jsou zobrazeny: •
Vliv velikosti čerpaných množství podzemní vody na hladinu podzemní vody je zobrazen na obrázcích 9, 14, 19, 24, 29, 34, 39, 44, 49, 54, 59, 64.
•
Vliv velikosti čerpaných množství na hladinu podzemní vody v domovních studnách ve Veltrusech zobrazují grafy 5 – 15 a tabulky 7 – 18.
•
Vliv velikosti čerpaných množství na šíření znečištění je zobrazeno na obrázcích 10 -13, 15 – 18, 20 – 23, 25 – 28, 30 – 33, 35 – 38, 39 – 43, 45 – 48, 50 – 53, 55 – 58, 60 – 63, 65 – 68.
Pro přehlednost následujícího textu uvádím, že příslušný text popisující každou sérii obrázků je umístěn vždy pod celou sérií, přičemž každý modelovaný příklad je v následujícím textu rozdělen do dvou částí. V první části popisuji vliv čerpaných množství podzemní vody na hydrodynamiku celé lokality, což je doloženo mapou hydroizohyps a vliv na HPV v domovních studnách ve Veltrusech, což je doloženo grafem a tabulkou. V druhé části popsuji vliv výše zmíněné hydrodynamiky na šíření znečištění z areálu rafinerie a účinnost HOPV, což dokládá série obrázků, na kterých jsou zobrazeny izolinie znečištění v různých časových krocích. Obrázky jsem vytvořil pomocí nástroje 2D Visualization v programu PMWIN Pro, nebo v programu Surfer, ve kterém jsem rovněž vytvářel obrysové mapy zobrazující situaci na dané lokalitě, které jsem potom ve formě *.dxf souboru použil při vytváření obrázků v programu PMWIN Pro.
45
8.1 Zobrazení a popis modelovaných situací
Graf 5 Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech bez čerpání z HOPV Název studny
Obr.9 hydroizohypsy bez čerpání z HOPV
Vypočítaná Pozorovaná hodnota hodnota
25
17.99675
17
17.99343
model bez 15.84 čerpání bez 15.67 čerpání
Tab. 8 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech bez čerpání z HOPV
Pro Porovnání šíření kontaminantů ovlivněný provozem HOPV, jsem považoval za nutné nejprve modelovat situaci, kdy nedochází k ovlivnění hydrodynamických podmínek provozem HOPV. Proudění podzemních vod je ovlivněno hlavně řekou Vltavou, Miřejovickým jezem, Mlýnským potokem a přítokem z horních kvartérních teras. Tuto skutečnost dokládají i výsledky modelu patrné z Obr. 6 a Obr.9. Jak dokládají modelové výsledky, provoz HOPV má zásadní vliv na HPV v domovních studnách ve Veltrusech. Bez provozu HOPV jsou hladiny cca o 2 m vyšší než při jejím provozu, jak dokládá 46
Graf 5. a Tab. 8. . Do tohoto modelu nejsou zahrnována ani čerpaná množství z domovních studní ve Veltrusech, proto je v úrovní hladin tak markantní rozdíl.
Obr.10 Izolinie šíření kontaminantů bez čerpání HOPV po 1 měsíci
Obr.12 Izolinie šíření kontaminantů bez čerpání HOPV po 12 měsících
Obr.11 Izolinie šíření kontaminantů bez čerpání HOPV po 6 měsících
Obr.13 Izolinie šíření kontaminantů bez čerpání HOPV po 36 měsících
Jak dokládají výsledky modelu šíření kontaminantů z kontinuálního zdroje viz Obr. 9-13 , bez provozu HOPV se znečištění šíří převážně vlivem advekce pouze ve směru proudění podzemní vody. Teprve přibližně po dvou letech se znečištění dostane do řeky Vltavy. Z výsledků modelu vyplývá, že provoz HOPV 47
je proto nutný. Základním problémem ovšem zůstává kde a kolik čerpat, aby nedocházelo k úniku znečištění z areálu, ale zároveň nedocházelo přitahování kontaminantů proti směru proudění podzemní vody a tím ke zvětšování znečištěné plochy.
Graf 6. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech při čerpání 100 % z objektů HOPV Název studny
Obr.14 hydroizohypsy při 100 % čerpání z HOPV
Vypočítaná hodnota
25
15.8558
17
15.59847
Pozorovaná hodnota
model 100 % 15.84 čerpání 100 % 15.67 čerpání
Tab. 9 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech při čerpání 100 % z objektů HOPV
Jak již bylo zmíněno v kapitole 6. proudění podzemních je nejvíce ovlivněno Vltavou, Miřejovickým jezem, přítokem z horních teras a provozem HOPV, což potvrzují i výsledky úlohy, zobrazené na Obr.14, při zadání konkrétních hodnot čerpaných množství podzemní vody získaných z reálných měření (Kliner, 2006). Z modelového příkladu vyplývá, že na proudění podzemních vod má i vliv čerpání z domovních studen ve Veltrusech. V dalších modelových příkladech 48
bude řešen vliv HOPV
na hladinu podzemních vod ve Veltrusech a
optimalizace čerpání, aby docházelo k co nejmenšímu ovlivnění hladiny ve studnách a zároveň k co nejmenšímu pohybu kontaminantů.
Obr.15 Izolinie šíření kontaminantů při 100 % čerpání HOPV po 1 měsíci
Obr.17 Izolinie šíření kontaminantů při 100 % čerpání HOPV po 12 měsících
Obr.16 Izolinie šíření kontaminantů při 100 % čerpání HOPV po 6 měsícíc
Obr.18 Izolinie šíření kontaminantů při 100 % čerpání HOPV po 24 měsících
Řešení šíření kontaminantů z kontinuálního zdroje při čerpání reálných objemů z HOPV zobrazené na Obr. 15 – 18 dokládá, že proti šíření bez provozu HOPV dochází k přitahování kontaminantů ke kralupské větvi HOPV a tím ke
49
zvětšování kontaminované plochy. Z výsledků modelu vyplývá, že z areálu uniká pouze jedna tisícina zadané koncentrace znečišťující látky.
Graf 7. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech pří čerpání 80 % z objektů HOPV
Vypočíta Název ná studny hodnota
Obr.19 hydroizohypsy při 80 % čerpání z HOPV
25
15.96388
17
15.77547
Pozorovaná hodnota
Model 80 % 15.84 čerpání 80 % 15.67 čerpání
Tab. 10 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech pří čerpání 80 % z objektů HOPV
V dalších modelovaných příkladech jsem řešil vliv čerpaného množství na HPV v domovních studnách a na rychlost a směr šíření znečištění. Hydroizohypsy HPV při čerpání 80 % původního objemu ze všech objektů HOPV
jsou zobrazeny na Obr. 19. Při čerpání 80 % z HOPV již dochází
k nárůstu HPV ve Veltrusech cca o 10 cm, jak dokládá graf 7 a Tab. 10.
50
Obr.20 Izolinie šíření kontaminantů při 80 % čerpání HOPV po 1 měsíci
Obr.22 Izolinie šíření kontaminantů při čerpání HOPV po 12 měsících
80 %
Obr.21 Izolinie šíření kontaminantů při 80 % čerpání HOPV po 6 měsících
Obr.23 Izolinie šíření kontaminantů při 80 % čerpání HOPV po 24 měsících
Změna čerpaného množství o 20 % na šíření kontaminantů prakticky nemá vliv, jak dokládají Obr. 20 – 23. Dochází k rychlejšímu úniku znečištění přes HOPV, ale pouze v jedné tisícině zadávaného množství. Z modelu je zřejmé, že stále dochází k přitahování kontaminantů kralupskou větví HOPV.
51
Graf 8. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech pří čerpání 50 % z objektů HOPV
Název studny
Obr.24 Hydroizohypsy při 50 % čerpání z HOPV
Vypočítaná Pozorovaná hodnota hodnota
25
16.12396
17
15.89093
model 50 % 15.84 čerpání 50 % 15.67 čerpání
Tab. 11 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech pří čerpání 50 % z objektů HOPV
Při snížení čerpání o 50 % již HPV ve Veltrusech stoupá cca o 30 cm.
52
Obr.25 Izolinie šíření kontaminantů při 50 % čerpání HOPV po 1 měsíci
Obr.27 Izolinie šíření kontaminantů při 50 % čerpání HOPV po 12 měsících
Obr.26 Izolinie šíření kontaminantů při 50 % čerpání HOPV po 6 měsících
Obr.28 Izolinie šíření kontaminantů při 50 % čerpání HOPV po 24 měsících
V této skupině úloh se ukázalo, že při snížení čerpání z HOPV na 50 % původního množství již dochází k většímu úniku kontaminantů což je zobrazeno na Obr. 25 - 28. Přibližně po šesti měsících se znečištění dostává až do Veltrus, kde může dojít k zasažení domovních studní.
Z tohoto
modelu lze usoudit, že provoz HOPV je nutný a při dalším postupu 53
optimalizace čerpaného množství se bude třeba zaměřit spíše na množství čerpané vody z jednotlivých vrtů.
Graf 9. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech pří čerpání 20 % z objektů HOPV
Název Vypočítaná studny hodnota
Obr.29 Hydroizohypsy při 20 % čerpání z HOPV
25
16.28225
17
15.90651
Pozorovaná hodnota
model 20 % 15.84 čerpání 20 % 15.67 čerpání
Tab. 12 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech pří čerpání 20 % z objektů HOPV
Hodnoty výšky HPV, při čerpání pouze 20 % množství čerpaného z objektů HOPV, ve Veltrusech jsou cca o 40 cm vyšší viz Graf 9. a Tab. 12 a již se blíží hodnotám vypočítaných modelem, ve kterém se simuluje stav bez ovlivnění HOPV.
54
Obr.30 Izolinie šíření kontaminantů při 20 % čerpání HOPV po 1 měsíci
Obr.32 Izolinie šíření kontaminantů při % čerpání HOPV po 6 měsících
20
Obr.31 Izolinie šíření kontaminantů při 20 % čerpání HOPV po 6 měsících
Obr.33 Izolinie šíření kontaminantů při 20 % čerpání HOPV po 24 měsících
Model šíření kontaminantů při 20 % provozu HOPV potvrzuje účinnost HOPV a jak je z výsledků modelů zobrazených na Obr. 30 - 33 zřejmé, při takto omezeném provozu HOPV již dochází k rozsáhlému úniku kontaminantů z areálu rafinerie.
55
Graf 10. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech pří čerpání 150 % z objektů HOPV
Název Vypočítaná studny hodnota
Obr.34 hydroizohypsy při 150 % čerpání z HOPV
25
15.64495
17
15.43505
Pozorovaná hodnota
model 150 % 15.84 čerpání 150 % 15.67 čerpání
Tab. 13 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech pří čerpání 150 % z objektů HOPV
Při zadání o 50 % většího objemu čerpaného množství dochází k dalšímu snížení HOPV o dalších cca 20 cm jak dokládá Graf 10. a Tab. 13. Z Obr. 34 je zřejmé, že v místě terénní elevace dochází k přerušení hydroizohyps, což lze vysvětlit tím, že terénní elevace je vyšší než HPV a dochází k zvětšení nesaturované
zóny
až
na
dno
první
modelové
vrstvy.
Z výsledků
hydrodynamického modelu jsou již hladiny v domovních studnách ve Veltrusech významně nižší, jak dokládá graf 10. a Tab. 13. což by mohlo mít za následek nedostatečnou zásobu podzemní vody ve Veltrusech.
56
Obr.35 Izolinie šíření kontaminantů 150 % čerpání HOPV po 1 měsíci
při
Obr.36 Izolinie šíření kontaminantů při 150 % čerpání HOPV po 6 měsících
Obr.37 Izolinie šíření kontaminantů při 150 % čerpání HOPV po 12 měsících
Obr.38 Izolinie šíření kontaminantů při 150 % čerpání HOPV po 24 měsících
Jak dokládají modelové výsledků zobrazené na Obr. 35 - 38 při čerpání o 50 % většího množství již nedochází k prakticky žádnému úniku znečištění přes HOPV, což potvrzuje důležitost HOPV při ochraně podzemních vod.
57
Graf 11. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech při čerpání pouze z kralupské větve HOPV
Název Vypočítaná studny hodnota
Obr.39 Hydroizohypsy při čerpání pouze z kralupské větve HOPV
25
16.2396
17
15.890304
Pozorovaná hodnota
model Kral. 15.84 čerpání Kral. 15.67 čerpání
Tab. 14 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech při čerpání pouze z kralupské větve HOPV
V dalších modelovaných příkladech jsem se zaměřil na omezení čerpání z některých vrtů a Porovnání jejich vlivu na hladinu ve Veltrusech a na funkci HOPV co se týče šíření kontaminantů. Při omezení čerpání pouze z kralupské větve HOPV dochází ke zvýšení HPV ve Veltrusech cca o 40 cm, z čehož vyplývá, že na snížení HPV v domovních studnách ve Veltrusech má hlavní vliv čerpání z veltruské větve HOPV. Ta má ovšem zásadní vliv na únik kontaminantů z areálu rafinerie, takže při další optimalizaci bude nutné regulovat spíše čerpaná množství z kralupské větve HOPV.
58
Obr.40 Izolinie šíření kontaminantů čerpání kralupské větve HOPV po1 měsíci
Obr.42 Izolinie šíření kontaminantů čerpání kralupské větve HOPV po12 měsících
Obr.41 Izolinie šíření kontaminantů čerpání kralupské větve HOPV po 6 měsících
Obr.43 Izolinie šíření kontaminantů čerpání kralupské větve HOPV po 24 měsících
Z modelu šíření kontaminantů při provozu pouze kralupské větve dochází k úniku kontaminantů v směru proudění podzemní vody, což je zobrazeno na Obr. 40 – 43. K jejich zadržení dochází hlavně ve směru ke kralupské větvi HOPV.
59
Graf 12. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech při čerpání pouze z veltruské větve HOPV
Obr.44 Hydroizohypsy při čerpání pouze z veltruské větve HOPV
Název Vypočítaná studny hodnota
Pozorovaná hodnota
25
16.00858
15.84
17
15.67918
15.67
model Vel. čerpání Vel. čerpání
Tab. 15 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech při čerpání pouze z veltruské větve HOPV
Z Výsledků modelu provozu pouze veltruské větve HOPV zobrazených na Obr. 44 je patrné, že oproti stavu se 100 % čerpáním z objektů HOPV je hladina v domovních studnách ve Veltrusech vyšší jak dokládá Graf 12. a Tab. 15, takže vliv na její úroveň má i čerpání z kralupské větve HOPV.
60
Obr.45 Izolinie šíření kontaminantů čerpání veltruské větve HOPV po 1 měsíci
Obr.47 Izolinie šíření kontaminantů čerpání veltruské větve HOPV po12 měsících
Obr.46 Izolinie šíření kontaminantů čerpání veltruské větve HOPV po 6 měsících
Obr.48 Izolinie šíření kontaminantů čerpání veltruské větve HOPV po 24 měsících
Jak dokládají výsledky modelu šíření kontaminantů při provozu pouze veltruské větve HOPV zobrazené v Obr. 45 – 48, dochází k výraznému zpomalení šíření kontaminantů, hlavně ve směru k řece Vltavě, protože tam nepůsobí vliv vrtů umístěných právě v tomto směru. Při provozu pouze Veltruské větve HOPV sice 61
dochází k úniku jedné tisíciny původní koncentrace znečištění, ale jak dokládají modelové výsledky, k tomu docházelo i při 100 % čerpání z objektů HOPV. Při uvážení všech faktorů mající vliv na provoz HOPV jako je ochrana podzemních vod, snižování úrovně HPV ve Veltrusech a v neposlední řadě finanční aspekt provozu HOPV, se toto řešení jeví jako ideální, protože úrovně hladin ve Veltrusech jsou vyšší a znečištění se z areálu nešíří a je zachytáváno veltruskou větví HOPV. Ovšem protože areál rafinerie je rozsáhlý a je tam mnoho potenciálních zdrojů znečištění nebylo by vhodné provoz kralupské větve zcela omezit.
Graf 13. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech při čerpání 150 % z veltruské větve HOPV
Název Vypočítaná Pozorovaná studny hodnota hodnota model 25 15.81378 15.84 Vel.150% 17 15.52572 15.67 Vel.150%
Obr.49 Hydroizohypsy při čerpání % z veltruské větve HOPV
150
Tab. 16 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech při čerpání 150 % z veltruské větve HOPV
Při zvýšení čerpaného množství z veltruské větve HOPV o 50 % již dochází ke snížení hladin přibližně na úrovně při čerpání 100 % z celé HOPV.
62
Obr.50 Izolinie šíření kontaminantů čerpání 150 % z veltruské větve HOPV po 1 měsíci
Obr.52 Izolinie šíření kontaminantů čerpání 150 % z veltruské větve HOPV po 12 měsících
Obr.51 Izolinie šíření kontaminantů čerpání 150 % z veltruské větve HOPV po 6 měsících
Obr.53 Izolinie šíření kontaminantů čerpání 150 % z veltruské větve HOPV po 24 měsících
Z modelových výsledků zobrazených na Obr. 50 – 53 je patrné, že při čerpání 150 % z Veltruské větve HOPV dochází k minimálnímu úniku znečištění mimo areál
rafinerie
a
kontaminační
mrak
se
výrazně
nerozšiřuje
do
nekontaminovaných míst. 63
Graf 14. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech při čerpání 100 % z kladrubské a 80 % z veltruské větve HOPV
Název studny
Obr. 54 Hydroizohypsy při čerpání 100 % z kladrubské a 80 % z veltruské větve HOPV
Vypočítaná Pozorovaná hodnota hodnota model Kral. 100% 25 15.8558 15.84 Vel.80% Kral. 100% 17 15.59847 15.67 Vel.80% Tab. 17 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech při čerpání 100 % z kladrubské a 80 % z veltruské větve HOPV
Z výsledků modelu zobrazených na Obr. 54, do kterého byly zadány hodnoty čerpání na 100 % z kralupské a 80 % z veltruské větve HOPV vyplývá že i při čerpání 80 % z veltruské větve HOPV dochází ke snížení hladin v domovních studnách a hodnoty se výrazně neliší od hodnot při celkovém 100 % čerpání jak dokládá graf 14. a Tab. 17.
64
Obr.55 Izolinie šíření kontaminantů čerpání 100 % z kladrubské a 80 % veltruské větve HOPV po1 měsíci
Obr.57 Izolinie šíření kontaminantů čerpání 100 % z kladrubské a 80 % veltruské větve HOPV po 12 měsících
Obr.56 Izolinie šíření kontaminantů čerpání 100 % z kladrubské a 80 % veltruské větve HOPV po 6 měsících
Obr.58 Izolinie šíření kontaminantů čerpání 100 % z kladrubské a 80 % veltruské větve HOPV po 24 měsících
Z modelu šíření znečištění, jehož výsledky jsou zobrazeny na Obr. 55 – 58 vyplývá, že nedochází k výraznému úniku znečištění přes HOPV. Je zde ale patrný vliv většího čerpání z kralupské větve HOPV, protože dochází k přitahování kontaminačního mraku směrem ke kralupské větvi. 65
Graf 15. Porovnání pozorovaných a vypočítaných hodnot piezometrické výšky v domovních studnách ve Veltrusech při čerpání 100 % z veltruské a 80 % z kladrubské větve HOPV
Název studny
Obr.59 Hydroizohypsy při čerpání 100 % z veltruské a 80 % z kladrubské větve HOPV
Vypočítaná hodnota
25
15.88602
17
15.61422
Pozorovaná hodnota
model Kral.80% 15.84 Vel.100% Kral.80% 15.67 Vel.100%
Tab. 18 Vypočítané a naměřené hodnoty piezometrických výšek v domovních studnách ve Veltrusech při čerpání 100 % z veltruské a 80 % z kladrubské větve HOPV
Jedním z možných řešení optimalizace HOPV je zanechání stejného množství čerpání podzemní vody z veltruské větve HOPV, aby se zamezilo šíření znečištění a postupné zmenšování čerpaných množství z kralupské větve HOPV. Z výsledků modelového řešení zobrazeného na Obr. 59 při čerpání 100 % objemu čerpaných množství z veltruské a 80 % z kralupské větve HOPV je zřejmé, že se ještě hladiny v domovních studnách významně neliší což dokládá graf 15. a Tab. 18.
66
Obr.60 Izolinie šíření kontaminantů čerpání 80 % z kladrubské a 100 % veltruské větve HOPV po1 měsíci
Obr.62 Izolinie šíření kontaminantů čerpání 80 % z kladrubské a 100 % veltruské větve HOPV po 12 měsících
Obr.61 Izolinie šíření kontaminantů čerpání 80 % z kladrubské a 100 % veltruské větve HOPV po 6 měsících
Obr.63 Izolinie šíření kontaminantů čerpání 80 % z kladrubské a 100 % veltruské větve HOPV po 24 měsících
Z výsledků modelu šíření kontaminantů zobrazených na Obr. 60 – 63 je zřejmé, že nedochází k rozsáhlému šíření znečištění ani při zmenšení čerpaného množství z kralupské větve HOPV o 20 % proti reálnému stavu.
67
Graf 16. změny hladin v pozorovacích vrtech při zvýšení hladiny podzemní vody
Obr.64 Hydroizohypsy při simulaci zvýšené hladiny podzemní vody
Protože rafinerie leží v blízkosti řeky Vltavy, hrozí areálu postižení záplavami. Při náhlém zvýšení HPV hrozí vyplavení znečištění z nesaturované zóny a zvýšení koncentrace. Při simulaci záplav jsem zvýšil hladinu v řece Vltavě z původních 168.5 m. n. m na 170.5 m. n. m. Jak dokládá 0br. 64 hydroizohypsy více kopírují řeku Vltavu a Miřejovický potok.
68
Obr.65 Izolinie šíření kontaminantů po 1 měsíci při simulaci zvýšeni HPV
Obr.67 Izolinie šíření kontaminantů po 12 měsících při simulaci zvýšeni HPV
Obr.66 Izolinie šíření kontaminantů po 6 měsících při simulaci zvýšeni HPV
Obr.68 Izolinie šíření kontaminantů po 24 měsících při simulaci zvýšeni HPV
Jak je patrné z výsledků modelu zobrazených na Obr. 65 - 64, při zvýšení HPV již provoz HOPV ztrácí na účinnosti a kontaminační mrak přes HOPV uniká směrem k Vltavě. Z tohoto důvodu by bylo nutné při povodňových situacích zvýšit čerpané množství, tak aby HOPV plnila svou funkci.
69
9. Závěr
Hlavní cíle diplomové práce byly: 1) Ovládnutí programu PMWIN Pro a popsání jeho výhod a nevýhod. 2) Posouzení
využití
metod
matematického
modelování
při
eliminaci
ekologických škod na podzemních vodách. 3) Optimalizace provozu hydraulické bariéry na konkrétním případě rafinerie Kaučuk Kralupy a.s.
1) Ovládnutí programu PMWIN Pro a popsání jeho výhod a nevýhod. a) Nevýhody programu: Mezi hlavní nevýhody programu PMWIN pro podle mého názoru patří: •
Problémy související s používáním metody konečných diferencí, která je popsána v kapitole 5.1. Jednou z hlavních nevýhod z MKD vycházející je nemožnost detailně popsat terénní elevace, které mají vliv jak na hydrodynamiku dané oblasti, tak na řešení úloh zabývajícími se šířením znečištění. Při zadání náhlých změn parametrů, terénních skoků nebo náhlých změn velikosti čerpání vody z vrtů dochází k značné nestabilitě výpočtu, nebo výpočet vůbec neproběhne, jak jsem si mohl ověřit při modelování vlivu terénní elevace.
•
Omezení možnosti definování pouze 9000 výpočetních buněk, což vede k problémům při sestavování vícevrstvých modelů.
•
Při sestavování vícevrstvých modelů lze definovat pouze první vrstvu jako vrstvu s volnou hladinou podzemní vody.
•
Jak bylo popsáno výše, při počítání vrstev s volnou hladinou podzemní vody má
model problémy při zadání velkých
piezometrických
výšek
na
okrajových
rozdílů
podmínkách
čili
hodnot velkého
hydraulického gradientu. •
Jednou z největších nevýhod programu je nedokonalý modul na popisování vrtů. Vrty se zadávají do celé buňky, takže čerpaná množství jsou počítána pro celé oblasti buňek, což u modelů popisujících rozsáhlé oblasti může tvořit i několik set metrů čtverečních. Co se týče vlastního zadávání vlastností vrtu, není možné zadávat míru kolmatace vrtu a vliv jeho polohy ve výpočetní buňce.
70
•
Důležitou součástí modelovacích programů jsou nástroje pro vizualizaci dat. Nástroj 2D Visualization je schopen na základě vypočítaných dat vykreslovat izolinie hladin podzemních vod, snížení nebo koncentrace kontaminantů. Nástroj 2D Vizualization ovšem neumožňuje vkládat do map měřítko a legendu. Dále při zadávání parametrů pomocí nástroje Polyline nezobrazuje zadávané hodnoty při ukládání do formátu *.bmp . Nástroj neukládá kontinuální vizualizaci dat do jednoho souboru (např *.avi) takže výsledky vizualizací není možné zobrazovat v jiných programech.
•
Důležitou pomůckou při učení se ovládání a pochopení programu je přikládaný manuál, který je sice přehledný ale nepopisuje v dostatečné míře používané nástroje a na mnoho věcí musí uživatel přijít sám.
b) Výhody programu: Mezi hlavní výhody programu PMWIN pro podle mého názoru patří: •
Možnost použití volně dostupné verze na internetu a tudíž dostupnost pro řešení širokého okruhu ekologických úloh.
•
Mezi hlavní výhody je jeho přehlednost a intuitivnost, která je významná při rychlém zpracování jednotlivých variant.
•
Program se snadno ovládá a práci ulehčuje možnost použití mnoha klávesových zkratek. Kladem programu je možnost zadávání parametrů pomocí polygonů, polyliline nebo přímo do buněk. Hodnoty se mohou přímo kopírovat do jakékoliv vrstvy, což usnadňuje práci při vytváření vícevrstvých modelů.
•
Data se do modelu mohou zadávat rovněž ze souborů vytvořených v jiných programech např. ve formátu ASCII.txt nebo *.Srf. Do těchto typů souborů se také mohou ukládat vypočítané výsledky a je tedy možné výsledky dále použít v dalších programech.
•
Mezi nesporné výhody tohoto programu
patří velké
množství
podporovaných modulů sloužících k počítání šíření znečištění (MT3D, MT3DMS ,MOC3D, RT3D, PHT3D), zadávání vlastností prostředí (řeka, vrt, drén, vodní plocha, evapotranspirace, infiltrace), programů sloužících pro inverzní modelování a kalibraci (MODFLOW – 2000, PEST, UCODE), nebo generování dat (digitizer, field interpolátor, field generator, result extraktor) . 71
2) Posouzení využití metod matematického modelování při eliminaci ekologických škod na podzemních vodách.
I když vytvoření modelu proudění podzemní vody a šíření znečištění jako odraz reálné situace je velice obtížné, výsledky mnou sestavovaného modelu proudění podzemní vody a šíření znečištění pro danou lokalitu dokládají, že matematické modelování je velice silný nástroj při modelování stávajících stavů, optimalizačních úlohách a jako nejdůležitější v úlohách prognózních. Využití správně sestaveného a nakalibrovaného modelu může investorovi ušetřit nemalé finanční prostředky související s čerpáním a ochranou podzemních vod a významně snížit časové nároky na sanaci znečištěných území.
3) Optimalizace provozu hydraulické bariéry na konkrétním případě rafinerie Kaučuk Kralupy a.s.
Výsledky modelu proudění podzemních vod dokazují, že hladina a průběh hladin a tím i proudění podzemní vody jsou zásadně ovlivněny řekou Vltavou, Miřejovickým jezem a přítokem vody z horní kvartérní zvodně. Model šíření znečištění dále ukazuje na možnost šíření znečištění hlavně severním směrem a to k řece Vltavě. Znamená to podle mého názoru, že další provoz hydraulické ochrany podzemních vod je nadále nutný. Z výsledků úloh řešících optimalizaci čerpání vody z objektů HOPV a jeho vliv na šíření znečištění a hladinu podzemních vod v domovních studnách se jako nejlepší řešení ukazuje možnost čerpání vody pouze z veltruské větve HOPV (severní) a tím dosáhnutí zvýšení hladin v domovních studnách v obci Veltrusy a tím i snížení finančních nákladů na provoz vrtů kralupské větve HOPV (jižní větev). Z důvodů malé informační zabezpečenosti modelu a celkem reálného předpokladu, že dochází k úniku znečištění i mimo modelovaná místa bych nedoporučoval úplné odstavení vrtů z kralupské větve HOPV, ale pouze snížení čerpaných množství a tím i omezení nákladů na jejich provoz. Dále doporučuji, aby byly vyhloubeny nové jímací vrty přímo v areálu rafinerie, které zajistí, aby nedocházelo ke zbytečnému postupu kontaminačního mraku proti směru proudění podzemních vod. 72
10. Seznam literatury
BEVEN, K. J., 2001: Rainfall runoff modelling. John Wiley & Sons, 372 s. CÍSLEROVÁ, M., VOGEL, T., 1998: Transportní procesy. Skriptum ČVUT, Praha, 182 s. DINGMAN, L. S., 2002: Physical hydrology. Prentice Hall, 456 s. DOHERTY, J., BREBBER, L., WHYTE P.,1994: PEST - Model-independent parametr estimation. User’s manual, Watermark Computing, Australia, 225 s. HRÁDEK, F. ,KUŘÍK, P., 2002: Hydrologie. Skriptum ČZU, Praha, 280 s. KLINER, K., 2006:Kaučuk Kralupy hydraulická ochrana, zpráva za rok 2006, Karel Kliner vodní zdroje. KOVÁŘ, P., 1990: Využití hydrologických modelů pro určování maximálních průtoků na malých povodích. VŠZ Praha, 136 s. McDONALD,
M.,
HARBAUGH,
A.,
1988:
MODFLOW,
A
modular
threedimensional finite difference ground-water flow model, U. S. Geological Survey, 492 s. MUCHA, I.,ŠESTAKOV, V., 1987:Hydraulika podzemných vod, Bratislava, Alfa, 302 s. Processing modflow Pro, User’s manual, 2002, WebTech 360, 409 s. ŘÍHA, J. et. al., 1997: Matematické modelování hydrodynamických a disperzních jevů. Skriptum VUT, Brno, 185 s. VALENTOVÁ, J., 2001: Hydraulika podzemní vody. Skriptu ČVUT, Praha, 174 s. ZHENG, C., WANG PP.,1999, MT3DMS: A modular three-dimensional multispecies model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems; Documentation and Users Guide, U.S. Army Engineer Research and Development Center, Vicksburg, MS. 146 s. Závěrečná zpráva: Vstupní audit stavu znečištění horninového prostředí a podzemních vod rafinérského komplexu v Kralupech., 1997, DHV Czech Republic, Zpráva o vlivu Kaučuk, a.s., na zdraví, bezpečnost a životní prostředí za rok 2000, Kaučuk a.s., 21 s.
73
Internetové zdroje
Matematické modelování šíření znečištění, 2006, Přednášky k předmětu: Transportní procesy a ochrana vod [online]. Dostupné z www.fle.czu.cz (cit. 2007-04-17) VLNAS, R., 2003: Hodnocení kvantitativního a chemického stavu podzemních vod [online]. Dostupné z www.chmu.cz/hydro/opzv/index.htm. (cit. 2007-03-12) http://www.scisoftware.com http://www.webtech360.com http://www.pmwin.com
Mapové podklady
http://www.mapy.cz http://geoportal.cenia.cz Výpočtová síť modelu, 1 : 5000, KAP spol. s.r.o., 1995 Geological Gross-Sections, 1: 5000, KAP spol. s.r.o., 1995 Iso-Concentration Contours of MTBE in Groundwater, 1: 5000, KAP spol. s.r.o., 1995
74
11. Seznam Příloh
Příloha 1. Geologické profily Příloha 1. 1 Geologický profil směr východ – západ Příloha 1. 2 Geologický profil směr jih – sever Příloha 2. Kalibrace parametrů řeky Příloha 2.1 Cbot= 0.01 m.s-1 Příloha 2.2 Cbot= 0.0001 m.s-1 Příloha 2. 3 Cbot= 0.00000001 m.s-1 Příloha 3. Čerpaná množství z objektů HOPV v [l.s-1] Příloha 4. Zadávané hodnoty čerpaných objemů do objektů HOPV v modelovaných příkladech Příloha 5. 3D zobrazení Příloha 5. 1 Zobrazení dna první vrstvy modelu Příloha5. 2 Hydroizohypsy při 100% čerpání vody z HOPV Příloha 5. 3 Hladina podzemní vody Příloha 5. 4 Izolinie koncentrace znečištění bez HOPV Příloha 6. Velikosti a směry vektorů rychlosti proudění podzemní vody (vypočítaných programem PMATH) Příloha 7. Ukázka vodní bilance modelu (vypočítaná pomocí modulu Water Budget) Příloha 8. Ukázka průběhu odhadu parametrů (pomocí programu PEST)
75
Příloha 1. Geologické profily Příloha 1.1 Geologický profil směr východ – západ
Příloha 1.2 Geologický profil směr jih - sever
76
Příloha 2. Kalibrace parametrů řeky Příloha 2.1
Cbot= 0.01 m.s-1
Příloha 2.2 Cbot= 0.0001 m.s-1
Příloha 2. 3 Cbot= 0.00000001 m.s-1
15
Příloha 3. Čerpaná množství z objektů HOPV v [l.s-1]
Čerpaná množství Q z objektů HOPV [ l.s-1]
8.00 7.00 6.00
Q [ l.s-1]
5.00 4.00 3.00
st1 st2 st3 st4 st5
2.00 1.00 Datum 0.00 29.12.200 17.2.2005 8.4.2005 28.5.2005 17.7.2005 5.9.2005 25.10.200 14.12.200 2.2.2006 4 5 5
16
Čerpaná množství Q z objektů HOPV [ l.s-1] 12.00
10.00
Q [ l.s-1]
8.00 st6 st7 st8 st9 st10
6.00
4.00
2.00
0.00 29.12.2004
17.2.2005
8.4.2005
28.5.2005
17.7.2005
5.9.2005
25.10.2005
14.12.2005
2.2.2006
Datum
Čerpaná množství Q z objektů HOPV [l.s-1] 9.00
8.00
7.00
Q [l.s-1]
6.00 st16 st17 st19 st20 st21
5.00
4.00
3.00
2.00
1.00
0.00 9.12.2004
28.1.2005
19.3.2005
8.5.2005
27.6.2005
16.8.2005
Datum
5.10.2005
24.11.2005
13.1.2006
4.3.2006
Příloha 4. Zadávané hodnoty čerpaných objemů do objektů HOPV v modelovaných situacích
VRT
3
-1
Q [m .s ] 100 %
3
-1
3
-1
3
-1
3
-1
Q [m .s ]
Q [m .s ]
Q [m .s ]
Q [m .s ]
80 %
50 %
20 %
150 %
St17B
-0.00203
-0.00162
-0.00102
-0.00041
-0.00305
St16A
-0.00375
-0.003
-0.00188
-0.00075
-0.00563
St1B
-0.00407
-0.00326
-0.00204
-0.00081
-0.00611
St2A
-0.00403
-0.00322
-0.00202
-0.00081
-0.00605
St3A
-0.00301
-0.00241
-0.00151
-0.0006
-0.00452
St4
-0.00367
-0.00294
-0.00184
-0.00073
-0.00551
St5
-0.00302
-0.00242
-0.00151
-0.0006
-0.00453
St7
-0.00302
-0.00242
-0.00151
-0.0006
-0.00453
St8A
-0.00303
-0.00242
-0.00152
-0.00061
-0.00455
St6B
-0.00418
-0.00334
-0.00209
-0.00084
-0.00627
St10B
-0.00404
-0.00323
-0.00202
-0.00081
-0.00606
St22
-0.0061
-0.00488
-0.00305
-0.00122
-0.00915
St19A
-0.0046
-0.00368
-0.0023
-0.00092
-0.0069
St20
-0.00185
-0.00148
-0.00093
-0.00037
-0.00278
St21
-0.0004
-0.00032
-0.0002
-0.00008
-0.0006
Příloha 5. 3D zobrazení Příloha 5.1 Zobrazení dna první vrstvy modelu
Terénní elevace
Příloha 5.2 Hydroizohypsy při 100% čerpání vody z HOPV
Příloha 5. 3 Hladina podzemní vody
Příloha 5. 4 Izolinie koncentrace znečištění bez HOPV
Příloha 6. Velikosti a směry vektorů rychlosti proudění podzemní vody (vypočítaných programem PMATH)
Příloha 7. Ukázka vodní bilance modelu (vypočítaná pomocí modulu Water Budget) PMWBLF (SUBREGIONAL WATER BUDGET) RUN RECORD FLOWS ARE CONSIDERED "IN" IF THEY ARE ENTERING A SUBREGION THE UNIT OF THE FLOWS IS [L^3/T] TIME STEP
1 OF STRESS PERIOD
1
============================================================= WATER BUDGET OF SUBREGIONS WITHIN EACH INDIVIDUAL LAYER ============================================================= REGION 1 IN LAYER 1 ---------------------------------------------FLOW TERM
IN
STORAGE
OUT
IN-OUT
0.0000000E+00 0.0000000E+00 0.0000000E+00
CONSTANT HEAD
0.0000000E+00 0.0000000E+00 0.0000000E+00
HORIZ. EXCHANGE
3.6698212E-03 0.0000000E+00 3.6698212E-03
EXCHANGE (UPPER) 0.0000000E+00 0.0000000E+00 0.0000000E+00 EXCHANGE (LOWER) 2.0107441E-07 0.0000000E+00 2.0107441E-07 WELLS
0.0000000E+00 3.6700000E-03 -3.6700000E-03
DRAINS
0.0000000E+00 0.0000000E+00 0.0000000E+00
RECHARGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
ET
0.0000000E+00 0.0000000E+00 0.0000000E+00
RIVER LEAKAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
HEAD DEP BOUNDS
0.0000000E+00 0.0000000E+00 0.0000000E+00
STREAM LEAKAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
INTERBED STORAGE 0.0000000E+00 0.0000000E+00 0.0000000E+00 RESERV. LEAKAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
SUM OF THE LAYER 3.6700224E-03 3.6700000E-03 DISCREPANCY [%]
2.2351742E-08
0.00
REGION 1 DOES NOT EXIST IN LAYER 2 -------------------------------------------------------------REGION 2 IN LAYER 1 ---------------------------------------------FLOW TERM STORAGE
IN
OUT
IN-OUT
0.0000000E+00 0.0000000E+00 0.0000000E+00
CONSTANT HEAD
0.0000000E+00 0.0000000E+00 0.0000000E+00
HORIZ. EXCHANGE
1.1241813E-03 1.1242233E-03 -4.2025931E-08
EXCHANGE (UPPER) 0.0000000E+00 0.0000000E+00 0.0000000E+00 EXCHANGE (LOWER) 4.1762767E-08 0.0000000E+00 4.1762767E-08 WELLS
0.0000000E+00 0.0000000E+00 0.0000000E+00
DRAINS
0.0000000E+00 0.0000000E+00 0.0000000E+00
RECHARGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
ET
0.0000000E+00 0.0000000E+00 0.0000000E+00
RIVER LEAKAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
HEAD DEP BOUNDS
0.0000000E+00 0.0000000E+00 0.0000000E+00
STREAM LEAKAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
INTERBED STORAGE 0.0000000E+00 0.0000000E+00 0.0000000E+00 RESERV. LEAKAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
SUM OF THE LAYER 1.1242231E-03 1.1242233E-03 -2.3283064E-10 DISCREPANCY [%]
0.00
WATER BUDGET OF SUBREGIONS OVER THE ENTIRE MODEL
REGION: 1 IN
OUT
IN-OUT
STORAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
CONSTANT HEAD
0.0000000E+00 0.0000000E+00 0.0000000E+00
HORIZ. EXCHANGE
3.6698212E-03 0.0000000E+00 3.6698212E-03
EXCHANGE (UPPER)
0.0000000E+00 0.0000000E+00 0.0000000E+00
EXCHANGE (LOWER)
2.0107441E-07 0.0000000E+00 2.0107441E-07
WELLS
0.0000000E+00 3.6700000E-03 -3.6700000E-03
DRAINS
0.0000000E+00 0.0000000E+00 0.0000000E+00
RECHARGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
ET
0.0000000E+00 0.0000000E+00 0.0000000E+00
RIVER LEAKAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
HEAD DEP BOUNDS
0.0000000E+00 0.0000000E+00 0.0000000E+00
STREAM LEAKAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
INTERBED STORAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
RESERV. LEAKAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
SUM OF REGIO 1)
3.6700224E-03 3.6700000E-03 2.2351742E-08
DISCREPANCY [%]
0.00
REGION: 2 IN
OUT
STORAGE
0.0000000E+00 0.0000000E+00
CONSTANT HEAD
0.0000000E+00 0.0000000E+00
HORIZ. EXCHANGE
1.1241813E-03 1.1242233E-03
IN-OUT 0.0000000E+00 0.0000000E+00 -4.2084139E-08
EXCHANGE (UPPER)
0.0000000E+00 0.0000000E+00 0.0000000E+00
EXCHANGE (LOWER)
4.1762767E-08 0.0000000E+00
WELLS
0.0000000E+00 0.0000000E+00 0.0000000E+00
4.1762767E-08
DRAINS
0.0000000E+00 0.0000000E+00 0.0000000E+00
RECHARGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
ET
0.0000000E+00 0.0000000E+00 0.0000000E+00
RIVER LEAKAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
HEAD DEP BOUNDS
0.0000000E+00 0.0000000E+00 0.0000000E+00
STREAM LEAKAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
INTERBED STORAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
RESERV. LEAKAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
SUM OF REGIO 2)
1.1242230E-03 1.1242233E-03 -3.4924597E-10
DISCREPANCY [%] 0.00 WATER BUDGET OF THE WHOLE MODEL DOMAIN: ============================================================= FLOW TERM STORAGE
IN
OUT
IN-OUT
9.3252018E-02 1.2185055E-01 -2.8598532E-02
CONSTANT HEAD
0.0000000E+00 0.0000000E+00 0.0000000E+00
WELLS
0.0000000E+00 4.7789998E-02 -4.7789998E-02
DRAINS
0.0000000E+00 0.0000000E+00 0.0000000E+00
RECHARGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
ET
0.0000000E+00 0.0000000E+00 0.0000000E+00
RIVER LEAKAGE
5.2386682E-02 2.8179750E-02 2.4206933E-02
HEAD DEP BOUNDS
5.2172564E-02 0.0000000E+00 5.2172564E-02
STREAM LEAKAGE 0.0000000E+00 0.0000000E+00 0.0000000E+00 INTERBED STORAGE 0.0000000E+00 0.0000000E+00 0.0000000E+00 RESERV. LEAKAGE
0.0000000E+00 0.0000000E+00 0.0000000E+00
-------------------------------------------------------------SUM 1.9781126E-01 1.9782031E-01 -9.0450048E-06 DISCREPANCY [%] 0.00
Příloha 8. Ukázka průběhu odhadu parametrů (pomocí programu PEST) PEST RUN RECORD: CASE PESTCTL PEST run mode:Parameter estimation mode Case dimensions:Number of parameters : 2 Number of adjustable parameters : 2 Number of parameter groups : 1 Number of observations : 13 Number of prior estimates : 0 Model command line(s):MODELRUN.BAT Jacobian command line:Na Model interface files:Templates: BCFTPL.DAT GHBTPL.DAT for model input files: BCF.DAT GHB.DAT (Parameter values written using single precision protocol.) (Decimal point omitted where possible to save space.) Instruction files: INSTRUCT.DAT for reading model output files: MODELOUT.DAT PEST-to-model message file:na Derivatives calculation:Param group g1
Increment Increment Increment Forward or Multiplier Method type low bound central (central) (central) relative 1.0000E-02 none always_2 na na
Parameter definitions:Name TransChange Initial Lower Upper formation limit value bound bound hk_1 log factor 8.938709E-04 9.999690E-06 9.999689E-02 ghb_1 log factor 1.781200E-05 1.781200E-07 0.178122 Name Group Scale Offset Model command number hk_1 g1 1.00000 0.00000 1 ghb_1 g1 1.00000 0.00000 1 Prior information:No prior information supplied Observations:Observation name Observation Weight Group 1 16.9600 1.000 no_name 2 16.9600 1.000 no_name 3 16.6200 1.000 no_name 4 16.7100 1.000 no_name 5 16.3800 1.000 no_name 6 16.4600 1.000 no_name
7 8 9 10 11 12 13
16.2100 15.7600 15.7800 15.6100 19.3400 15.8400 15.6700
1.000 1.000 1.000 1.000 1.000 1.000 1.000
no_name no_name no_name no_name no_name no_name no_name
Control settings:Initial lambda : 10.000 Lambda adjustment factor : 2.0000 Sufficient new/old phi ratio per optimisation iteration : 0.30000 Limiting relative phi reduction between lambdas : 1.00000E-02 Maximum trial lambdas per iteration : 8 Maximum factor parameter change (factor-limited changes) : 10.000 Maximum relative parameter change (relative-limited changes) : na Fraction of initial parameter values used in computing change limit for near-zero parameters : 1.00000E-03 Relative phi reduction below which to begin use of central derivatives : na Relative phi reduction indicating convergence : 0.10000E-01 Number of phi values required within this range : 3 Maximum number of consecutive failures to lower phi : 3 Minimal relative parameter change indicating convergence : 0.10000E-01 Number of consecutive iterations with minimal param change : 3 Maximum number of optimisation iterations : 0 OPTIMISATION RECORD INITIAL CONDITIONS: Sum of squared weighted residuals (ie phi) = 2.4352 Current parameter values hk_1 8.938709E-04 ghb_1 1.781200E-05 Optimisation complete: optimisation iteration limit of 0 realized. Total model calls: 1
OPTIMISATION RESULTS Covariance matrix and parameter confidence intervals cannot be determined:No optimisation iterations carried out so Jacobian never calculated. Parameters -----> Parameter Estimated value hk_1 8.938709E-04 ghb_1 1.781200E-05 See file PESTCTL.SEN for parameter sensitivities. Observations -----> Observation Measured Calculated Residual Weight Group value value 1 16.9600 17.0110 -5.104000E-02 1.000 no_name 2 16.9600 17.0513 -9.132000E-02 1.000 no_name 3 16.6200 16.4432 0.176830 1.000 no_name 4 16.7100 16.5257 0.184300 1.000 no_name 5 16.3800 16.0703 0.309680 1.000 no_name 6 16.4600 16.1751 0.284860 1.000 no_name 7 16.2100 15.6538 0.556190 1.000 no_name
8 9 10 11 12 13
15.7600 15.7800 15.6100 19.3400 15.8400 15.6700
15.2641 15.2852 15.1385 20.3091 15.5493 15.3013
0.495930 0.494840 0.471510 -0.969080 0.290670 0.368690
1.000 1.000 1.000 1.000 1.000 1.000
no_name no_name no_name no_name no_name no_name
See file PESTCTL.RES for more details of residuals in graph-ready format. See file PESTCTL.SEO for composite observation sensitivities. Objective function -----> Sum of squared weighted residuals (ie phi) = 2.435 Correlation Coefficient -----> Correlation coefficient = 0.9967 Analysis of residuals -----> All residuals:Number of residuals with non-zero weight = 13 Mean value of non-zero weighted residuals = 0.1940 Maximum weighted residual [observation "7"] = 0.5562 Minimum weighted residual [observation "11"] = -0.9691 Standard variance of weighted residuals = 0.2214 Standard error of weighted residuals = 0.4705 Note: the above variance was obtained by dividing the objective function by the number of system degrees of freedom (ie. number of observations with non-zero weight plus number of prior information articles with non-zero weight minus the number of adjustable parameters.) If the degrees of freedom is negative the divisor becomes the number of observations with non-zero weight plus the number of prior information items with non-zero weight.