Oscilace hladiny pozorované ve vrtu a jejich matematické modelování

Oscilace hladiny pozorované ve vrtu a jejich matematické modelování Jiří Mls1 , Tomáš Ondovčin1 1 Universita

Karlova v Praze, Přírodovědecká fakulta, Ústav hydrogeologie, inženýrské geologie a užité geofyziky, Praha – [email protected]

Abstrakt V článku jsou uvedeny dva možné modely mechanismu, který přenáší slapové síly působící na geologické těleso na změny napětí a tlakové výšky ve zkoumaném aquiferu. Následně je modifikována Biotova teorie poroelasticity a z jejích rovnic je odvozena řídící rovnice procesu, jejímž jedním parametrem je vektor slapového zrychlení. Úloha je řešena numericky a její výsledky jsou porovnány se zaznamenanými oscilacemi hladiny ve vrtu.

Klíčová slova: Biotova teorie poroelasticity, Darcyův zákon, oscilace hladiny, slapové zrychlení, slapy v podzemní vodě 1 ÚVOD Zatímco slapové jevy na mořském pobřeží a v podzemní vodě komunikující s mořem jsou známy od nepamněti, jejich pozorování v podzemní vodě ve vnitrozemí jsou relativně nového data. Prvním známým záznamem je (Klönne, 1880). Pozdější výzkum ukázal, že podstata vnitrozemských slapových jevů v podzemní vodě je principiálně odlišná od podstaty přímořských slapových jevů a souvisí spíše s deformací horninových těles provázenou změnou tlakové výšky než s působením slapových sil na podzemní vodu. Podle Rojstaczera a Agnewa (1989) jsou aquifery s velkou pórovitostí citlivější na na změny atmosferického tlaku, zatímco vysoká citlivost na deformace horniny se projevuje u aquiferů s nížší pórovitostí. U obou pak tato citlivost roste s klesající stlačitelností pevné fáze aquiferu. Obrázek 1: Hladina ve vrtech VS-20 a V-34 (Krásný et al. 2002)

Bodvarsson (1970) přidal k Darcyovu zákonu inerciální člen a formuloval rovnice deformací indukovaných tlakových fluktuací v aquiferech. Ve svém článku zdůrzňovat nezbytnost, aby vrt v

němž mají být fluktuace zaznamenány byl napojen na aquifer s dostatečnou napjatostí hladiny. Slapové fluktuace hladiny v různých vrtech zaznamenávají také Zeumann et al. (2008). Ve své práci dospívají ke zjištění, že hlubší vrty odrážejí zemské slapy s větší citlivostí a zřetelněji. Podobnou problematikou se zabývají také Ritzi et al. (1991), kteří do svého výzkumu navíc zahrnují také vliv atmosferického tlaku. Bernard a Delay (2008) porovnávají zaznamenané oscilace s oscilacemi slapových sil a využívají spektrální analýzu. Rojstaczer and Riley (1990) dospívají k závěru, že vrty, ve kterých nalézáme slapové oscilace hladiny, indikují lokality, kde napětí pevné fáze neodpovídá působícím slapovým silám. V takových místech tlak v kapalné fázi přebírá část napětí pevné fáze. Tento článek presentuje práci zaměřenou na oscilace hladiny podzemní vody zjištěné v Polické pánvi, které tam byly zjištěny v několika vrtech; Obr. 1 ukazuje dva takové záznamy. Počáteční předpoklad, že oscilace jsou vyvolány periodickým čerpáním z vrtu VS-13 byl potvrzen ve všech případech až na jednu výjimku. Oscilace ve vrtu pokračovaly beze změny i po zastavení čerpání a prokázaly nezávislost jak na Obrázek 2: Tektonika v okolí vrtu V-34 (Krásný et al. 2002)

čerpání tak na atmosférickém tlaku. Vzhledem k frekvenci oscilací, Krásný et al. (2002) došli k závěru, že podstatou pozorovaného jevu jsou slapové oscilace tlakové výšky podzemní vody. Ve své práci vycházíme ze znalosti specifických geologických podmínek polohy vrtu V-34. Navrhujeme dva modely mechanismu vyvolávajícího změny tlakové výšky v příslušném aquiferu a pomocí modifikované Biotovy teorie poroelasticity odvozujeme matematický model zkoumaného procesu.

2 GEOLOGICKÉ PODMÍNKY Polická pánev sestává z vrstev druhohorních pískovců, slínovců, prachovců a jílů, které byly během třetihor rozlámány na několik bloků. Na Obr. 2 je dobře vidět, že vrt V-34 se nachází v malém bloku pánve odděleném od ostatních ze tří stran velkými zlomy – Polickým zlomem, Skalským zlomem a zlomem rovnoběžným se Skalským zlomem. Vrt V-34 je hydraulicky propojen se zhruba 10m mocnou vrstvou silně rozpukaných rohovců, která je navzdory své malé mocnosti významným kolektorem. Na Obr. 2 vidíme také linii podélného řezu severní částí Polické pánve označenou AA’. Tento řez je zobrazen na Obr. 3 a je na něm vidět vzájemnou vertikální polohu jednotlivých částí druhohorních vrstev. Podrobné informace o geologii a hydrogeologii polické pánve lze nalézt v článku (Krásný et al. 2002). Obrázek 3: Podélný řez severní částí Polické pánve (Krásný et al. 2002)

2 MATEMATICKÝ MODEL 2.1 Rovnice poroelasticity Základem matematického modelu je předpoklad, že změny tlakové výšky jsou vyvolány změnami napětí v určité oblasti porosího prostředí. Protože pozorované fluktuace jsou blízké periodickým změnám, je na místě považovat vyvolanné deformace za elestické, viz (Ondovčin, Mls 2011) a (Ondovčin et al. 2012). Z tohoto důvodu je vhodným nástrojem pro vytvoření modelu Biotova teorie poroelasticity, (Biot 1941). Tato teorie vychází z Biotových předpokladů, které vyjadřují tensor τ napětí v pevné fázi porosního prostředí a objemový zlomek θ kapalné fáze (v našem případě vlhkost) ve tvaru   ∂uj ∂ui 2σµ 2 µ (1 + σ) + +  δi,j − p δi,j , i = 1, 2, 3 , (1) τi,j = ∂xj ∂xi 1 − 2σ 3 H (1 − 2 σ) a

 θ = θ0 + α  +

1 α − R H

 p,

kde x = (x1 , x2 , x3 ) jsou prostorové souřadnice, t je čas, u = (u1 , u2 , u3 ) je vektor posunutí hmotných bodů pevné fáze, µ je modul pružnosti ve smyku pevné fáze, σ je Poissonovo číslo

pevné fáze, p je tlak kapalné fáze, θ0 je vlhkost při nulovém tlaku a napětí,  = div u, H, R a α jsou konstanty, které zavedl Biot; 1/H je mírou objemové deformace porosního prostředí při jednotkové změně tlaku kapalné fáze, 1/R vyjadřuje změnu objemového zlomku vody vztaženou k jednotkové změně tlaku kapalné fáze. Koeficient α je definován vztahem α=

2 µ (1 + σ) . 3 H (1 − 2 σ)

Výsledná teorie se dá po určitých úpravách shrnout do následujících čtyř rovnic, jimiž je vektorová rovnice rovnováhy působících sil a rovnice kontinuity proudění podzemní vody, µ ∆ui + a

µ ∂ ∂p −α = 0, 1 − 2 σ ∂xi ∂xi

∂ α + ∂t



1 α − R H



i = 1, 2, 3 ,

∂p k = ∆p , ∂t η

(2)

(3)

kde k je propustnost aquiferu vystaveného působícímu napětí a η je dynamická viskosita kapalné fáze. Pro nasycené prostředí uvádí Biot (1941) tento vztah vztah α→1 a

HR → ∞ pro θ → n , H − αR

(4)

kde n je pórovitost. Vzhledem k tomu, že registrované oscilace hladiny odpovídají tlakovým změnám v rohovcovém aquiferu, bude oblastí našeho zájmu jeho část vymezená výše uvedenými zlomy. Na základě známých geologických podmínek lze považovat zlomy za nepropustné; tento předpoklad byl rovněž potvrzen provedenou zkouškou. Navíc můžeme předpokládat zanedbatelné změny hydraulické výšky ve směrech vertikálním a rovnoběžným se Skalským zlomem. Budeme vycházet z následujícího koncepčního modelu. Model definujeme v časoprostorové oblasti Ω=X ×T , kde X je jednorozměrný interval X = (0, L) orientovaný vodorovně kolmo na Skalský zlom. Jeho délka L je tedy vzdálenost mezi Skalským zlomem a zlomem s ním rovnoběžným. Časový interval T = (tB , tE )odpovídá délce záznamu oscilací, viz Obr. 1. V důsledku volby oblasti X budou případné vertikální přetoky mezi ní a nadložním nebo podložním aquiferem vyjadřovány jako zdrojové členy. Rovnice kontinuity bude mít tvar ∂θ ∂v + = q, ∂t ∂x kde θ je objemový zlomek vody (v uvažovaném případě pórovitost, protože se jedná o nasycené prostředí) a q je zdrojový člen, udávající objem vody, která se objeví v jednotkovém objemu prostředí za jednotku času. Zdrojovým členem je tedy hustota toku přetékající z nadloží dělená mocností aquiferu,   b k pb − p q= g ρw + , (5) ηZ Zb kde b k je propustnost nadložní polopropustné vrstvy, pb tlak kapalné fáze na horní hranici polopropustné vrstvy, Z je mocnost aquiferu a Zb je mocnost polopropustné vrstvy. Protože v jednodimensionálním případě je x1 = x ,

u1 = u ,

τ1,1 = τ ,

bude napětí v pevné fázi podle vztahu (1) vyjádřeno ve tvaru τ=

2 µ(1 − σ) ∂u − αp 1 − 2 σ ∂x

(6)

a rovnice rovnováhy působících sil (2) ve tvaru 2 µ(1 − σ) ∂ 2 u ∂p −α = 0. 2 1 − 2 σ ∂x ∂x

(7)

Po zavedení zdrojového členu a při vyjádření pro jednodimensionální problém se rovnice kontinuity (3) modifikuje na tvar   ∂2u 1 α ∂p k ∂2p +q. (8) α + − = ∂t ∂x R H ∂t η ∂x2 Poronáním rovnic (6) a (7) vidíme, že napětí v pevné fázi nezávisí na prostorové souřadnici, je tedy pouze funkcí času, τ = τ (t) . To znamená, že ke znalosti napětí v oblasti Ω stačí znát jeho hodnotu jako funkci času v jednom bodě intervalu X . Označme tedy T (t) časový průběh napětí v hraničním bodě x = 0 a T 0 (t) jeho derivaci podle času. Derivováním obou stran vztahu (6) podle času a jednoduchou úpravou dostáváme α (1 − 2 σ) ∂p 1 − 2σ ∂2u = + T0 . ∂t ∂x 2 µ(1 − σ) ∂t 2 µ(1 − σ) Tato rovnice nám umožňuje vyloučit druhou smíšenou derivaci z rovnice (8). Se současným využitím vztahů (4) a (5) pak dostáváme výslednou řídící rovnici zkoumaného procesu   ∂p 2 µ k (1 − σ) ∂ 2 p 2µb k (1 − σ) pb − p = + g ρw + − T0 . (9) ∂t η (1 − 2 σ) ∂x2 η Z (1 − 2 σ) Zb Jedná se o parabolickou parciální diferenciální rovnici druhého řádu s neznámou funkcí p, definovanou v oblasti Ω. Všechny parametry vystupující v této rovnici byly řádně definovány a až na derivaci T 0 napětí v bodě x = 0 je zřejmé, jak mohou být určeny. 2.2 Dva modely vzniku slapových jevů Napětí T (t) = τ (0, t) definované v intervalu T je dáno působením okolního horninového prostředí na zkoumaný blok x ∈ X v bodě x = 0. Z řídící rovice je zřejmé, že jedině toto napětí může vyvolávat oscilace tlaku p kapalné fáze v oblasti X . Existuje nějaký mechanismus, který časově proměnné slapové síly působící na sousední horninové prostředí transformuje na napětí T . V tomto odstavci uvedeme dva možné modely takového mechanismu. V obvyklých úlohách hydrauliky podzemní vody je zrychlení vyvolávající síly působící na vodu v pórech dáno vztahem g = (0, 0, −g), když souřadnicový systém je orientován tak, že osa x3 směřuje svisle vzhůru. V našem případě musíme ke gravitačnímu zrychlení vyvolávanému hmotností Země přidat slapové zrychlení. Omezíme se přitom na dvě tělesa s nejsilnějším účinkem, tj. na Měsíc a Slunce. Výsledné působící zrychlení tedy bude součtem g + γ, kde γ = γ(t) = (γ1 (t), γ2 (t), γ3 (t)) je slapové zrychlení vyvolávané v daném místě a čase uvednými dvěma tělesy. Vzhledem ke geologickým podmínkám v dané lokalitě lze navrhnout následující dva modely. První model určuje napětí v bodě x = 0 jako důsledek smykových deformací geologických bloků v bodech x = 0 a x = L se zkoumanou oblastí X . V důsledku hlubšího uložení jižního sousedního

bloku (x = 0) oproti severnímu (x = L) je magnituda smykové deformace v bodě x = 0 větší než v bodě x = L, orientace je stejná, napětí T vzniká rozdílem velokostí, viz Obr. 4. Působí-li horizontální složka slapového zrychlení ve směru osy x, je horninový blok X stlačován a funkce T nabývá záporných hodnot. Předpokládáme, že existuje konstanta MN , která vyjadřuje souhrnný účinek aktivně působící hmotnosti geologických těles zároveň s příslušným pákovým efektem. Za tohoto předpokladu je napětí T dáno vztahem Obrázek 4: Model č. 1 – napětí v oblasti X vytváří horizontální složka slapového zrychlení

T (t) = −MN (γ1 (t) ν1 + γ2 (t) ν2 ) , kde ν = (ν1 , ν2 ) je jednotkový vektor vnější normály oblasti X v bodě x = 0, tj. vektor kolmý na Skalský zlom orientovaný proti směru osy x. Druhý model určuje napětí v bodě x = 0 jako reakci dvěma sousedními zlomy přerušené vrstvy na vertikální složku slapového zrychlení. Svisle vzhůru orientovaná síla působící na obě křídla vrstev rozvírá prostor mezi nimi, tj. roztahuje oblast X , a tím v ní způsobuje kladné napětí; tento Obrázek 5: Model č. 2 – napětí v oblasti X vytváří vertikální složka slapového zrychlení

mechanismus je znázorněn na Obr. 5. Zavedením konstanty MN stejně jako výše dostáváme pro napětí T vyjádření T (t) = MN γ3 (t) , kde γ3 je, jak už bylo definováno, vertikální složka slapového napětí.

2.3 Formulace úloh a jejich řešení Problém na který se zaměřujeme je ověření, zda oscilace ve vrtu V-34, viz Obr. 1, mohou či nemohou být slapového původu. Určením mechanismu jímž vzniká napětí T v bodě x = 0, tedy přijetím jednoho z modelů uvedených v předchozím článku, máme definovány všechny parametry vystupující v řídící rovnici (9). Hranicí oblasti X , ve které úlohu řešíme jsou body x = 0 a x = L. Obě tyto části hranice odpovídají zlomům o nichž víme, že jsou v úrovni zkoumaného aquiferu nepropustné. Předepisujeme zde tedy nulový průtok hranicí; vzhledem k tomu, že neznámou funkcí je tlak p, bude okrajová podmínka vyjádřena ve tvaru ∂p (0, t) = 0 , ∂x

∂p (L, t) = 0 , ∂x

pro t ∈ T .

Počáteční podmínka je určena zjištěnou hodnotou pB tlaku v počátečním čase tB a může navíc zohledňovat fázi p(x, tB ) = pB , pro x ∈ X . b vystupujících v úloze byly určeny a pubHodnoty hydrogeologických parametrů k, b k, pb, Z a Z, likovány v práci (Krásný et al. 2002). Parametry µ a σ byly určeny pro daný druh křemence. Je prakticky nemožné určit explicitně konstantu MN , protože je závislá na neznámých částech sousedních geologických těles, jejich tuhosti a pákovém efektu. Tato konstanta nemá vliv na frekvenci působících slapových sil, ale je přímo úměrná jejich amplitudě. Dá se snadno určit kalibrací na měřenou amplitudu oscilací hladiny ve vrtu. Je to jediný parametr modelu, který byl kalibrován. Všechny ostatní byly určeny nezávisle na použitých modelech. Vektor slapového zrychlení γ byl pro daný čas a polohu určen výpočtem pomocí numerického kódu SP_SM_01, (Ondovčin 2007) využívajícího katalog slapového potenciálu (Hartmann a Wenzel 1995). Označíme-li p p , a h= , h= g ρw g ρw bude neznámou funkcí v řídící rovnici úlohy tlaková výška, což přímo odpovídá měřeným výškám hladiny ve vrtu; řídící rovnice bude ∂h 2 µ k (1 − σ) ∂ 2 p 2µb k (1 − σ) 2µb k (1 − σ) b b T0 = − h + h + Z − , b (1 − 2 σ) ∂t η (1 − 2 σ) ∂x2 η Z Zb (1 − 2 σ) g ρw ηZZ

(10)

počáteční podmínka bude h(x, tB ) = hB ,

pro x ∈ X .

(11)

a okrajová podmínka bude ∂h (0, t) = 0 , ∂x

∂h (L, t) = 0 , ∂x

pro t ∈ T .

(12)

Všimněme si nyní toho, že počáteční podmínka předepisuje neznámé funkci konstantní hodnotu nezávislou na x, že okrajová podmínka je homogenní a že také parametry úlohy nejsou závislé na prostorové proměnné. Vyplývá odtud, že každá funkce h daná vztahem h(x, t) = w(t) pro (x, t) ∈ X × T ,

(13)

je řešením úlohy (10), (11), (12), jestliže funkce w je v intervalu T řešením počáteční úlohy w0 = −

2µb k (1 − σ) 2µb k (1 − σ) b b T0 w+ h+Z − , b (1 − 2 σ) b (1 − 2 σ) g ρw ηZZ ηZZ w(tB ) = hB .

(14) (15)

Přirozeně vzniká otázka, jaká je množina řešení úlohy (10), (11), (12), která nejsou řešením úlohy (14), (15). Vzhledem k tomu, že počátešně okrajová úloha s parabolickou rovnicí tupu rovnice (10) a Neumannovou počáteční podmínkou je, za splnění některých dalších podmínek hladkosti, které v našem případě splněny jsou, jednoznačně řešitelná, viz (Vladimirov 1971), můžeme výše uvedený výrok zesílit: funkce h(x, t) je řešením úlohy (10), (11), (12) právě když může být vyjádřena ve tvaru (13), kde funkce w je řešením počáteční úlohy (14), (15). Do rovnice (14) můžeme zavést ještě tuto velmi užitečnou substituci: y(t) = w(t) +

T (t) . g ρw

Tímto způsobem převedeme úlohu (14), (15) na úlohu y0 = −

2µb k (1 − σ) 2µb k (1 − σ) 2µb k (1 − σ) b b y− T+ h+Z , b (1 − 2 σ) b (1 − 2 σ) ηZZ η Z Zb g ρw (1 − 2 σ) ηZZ y = hB +

T (tB ) . g ρw

(16)

(17)

Vzhledem k tomu jak je dána funkce T , je nutné řešit jak úlohu (14), (15) tak i úlohu (16), (17) numericky. Obě úlohy jsou si velmi podobné, ale důležitá výhoda řešení úlohy (16), (17) vůči řešení úlohy (14), (15) spočívá v tom, že nemusíme numericky derivovat funkci T . 3 VÝSLEDKY A DISKUSE Na Obr. 6 jsou vyneseny výsledky, které pro zkoumané místo a zkoumané období dává první model vzniku slapových oscilací. Pro srovnání je v obrázku zároveň vynesen záznam získaný ve vrtu V-34, známý už z Obr. 1. Při řešení úlohy byly všechny parametry modelů, s výjimkou konObrázek 6: Model č. 1 – srovnání výsledků s hodnotami měřenými ve vrtu

stanty MN , která nemá vliv na oscilatorický charakter výsledku, určeny nezávisle na naměřených výškách hladiny. Ze srovnání vypočtených hodnot s hotnotami měřenými vyplývá, že oscilace zjištěné ve vrtu V-34 jsou slapového původu. Totéž lze usuzovat i z výsledků získaných z druhého modelu, které jsou vyneseny na Obr. 7, opět pro srovnání společně se záznamem výšky hladiny ve vrtu V-34.

Kromě oscilací s amplitudou zhruba 10 – 15 mm je ze záznamu výšky hladiny ve vrtu patrný plynulý neoscilatorický vzestup hladiny probíhající od polovine druhého dne do konce pátého dne. Vzestup hladiny je následně vystřídán jen o trochu menším, rovněž neoscilatorickým, poklesem. Obrázek 7: Model č. 2 – srovnání výsledků s hodnotami měřenými ve vrtu

Tento vývoj je vysvětlován nárůstem hydraulické výšky ve zkoumané oblasti přetokem z aquiferu CV přes aquitard A/C, viz Obr. 3. I když výsledky potvrzují slapový původ oscilací, nedá se říci, že by některý z modelů přesně popisoval mechanismus vzniku oscilací tlakové výšky ve zkoumané oblasti. Z počáteční části zkoumaného období, zhruba první tři dny, lze usuzovat, že druhý model mnohem lépe vykresluje pozorovaný průběh oscilací. Je zde ale velký fázový rozdíl, který naopak mizí ve druhé části zkoumaného období. Dostáváme se k otázkám: je možné, aby významná změna (značně převyšující amplitudu oscilací) tlakových poměrů způsobila fázový posun reakce aquiferu na oscilace slapového zrychlení? Byl nárůst tlakové výšky skutečně důsledkem jejího nárůstu v sousedním aquiferu, nebo má jinou, např. seismickou, příčinu? Na tyto otázky nemohou předložené výsledky dát odpověď, jejich objasnění by vyžedovalo další výzkum. LITERATURA BERNARD, S., DELAY, F. Determination of porosity and storage capacity of a calcareous aquifer (France) by correlation and spectral analyses of time series. Hydrogeol. J. 2008, XVI. Nr. 7, pp. 1299–1309. BIOT, M. A. General theory of three-dimensional consolidation. J. Appl. Phys. 1941, XII. pp. 155–164. BODVARSSON, G. Confined fluids as strain meters. J. Geophys. Res. 1970, LXXV. Nr. 14, pp. 2711–2718. HARTMANN, T., WENZEL, H.-G. The HW95 tidal potential catalogue. Geophys. Res. Lett. 1995, XXII. pp. 3353–3556. KLÖNNE, F. W. Die periodischen Schwankungen des Wasserspiegels in den inundirten Kohlenschächten von Dux in der Periode vom 8. April bis 15. September 1879. Akad. Wiss. Wien Sitzungsber. 1880, LXXXI. pp. 114-116. KRÁSNÝ, J., BUCHTELE, J., ČECH, S., HRKAL, Z., JAKEŠ, P., KOBR, M., MLS, J., ŠANTRŮČEK, J., ŠILAR, J., VALEČKA, J. Hydrogeology of the Police Basin: optimisation of groundwater development and protection. J. Geol. Sci. 2002, XXII. pp. 5–100.

ONDOVČIN, T. Popis numerického kódu SPZ_SM_01. Přírodovědecká fakulta UK, Praha: 2007, Manuskript. ONDOVČIN, T., MLS, J. Relation between phase density and component concentration in groundwater flow modelling. Journal of Hydrology and Hydromechanics. 2007, LV. Nr. 4, pp. 236-245. ONDOVOČIN, T., MLS, J., HERRMANN, L. Mathematical modeling of tidal effects in groundwater. Transport in Porous Media. 2012, VC. Nr. 2, pp. 483-495. RITZI, R. W., SOROOSHIAN, S., HSIEH, P. A. The estimation of fluid-flowproperties from the response ofwater levels in wells to the combined atmospheric and Earth tide forces. Water Resour. Res. 1991, XXVII. Nr. 5, pp. 883–893. ROJSTACZER, S., AGNEW, D. C. The influence of formation material properties on the response of water levels in wells to earth tides and atmospheric loading Journal of Geophysical Research-Solid Earth and Planets. 1989, XCIV. Nr. B9, pp. 12403-12411. ROJSTACZER, S., RILEY, F.S. Response of the water level in a well to Earth tides and atmospheric loading under unconfined conditions. Water Resour. Res. 1990, XXVI. pp. 1803–1817. Vladimirov, V. S. Rovnice matematické fysiky. Moskva: Nauka, 1941 (v ruštině). ZEUMANN, S., WEISE, A., JAHR, T. Tidal and non-tidal signals in groundwater boreholes in the KTB area. Germany. J. Geodyn. 2009, IIL. pp. 115–119.