Technika TDR w modelowaniu ruchu wody glebowej : (rozprawa habilitacyjna)

TECHNIKA TDR

W MODELOWANIU

RUCHU WODY

GLEBOWEJ

TECHNIKA TDR

W MODELOWANIU

RUCHU WODY

GLEBOWEJ

Opiniodawca prof. dr hab. Henryk Sobczuk

Redaktor merytoryczny

dr hab. in. Krzysztof Pulikowski, prof. nadzw.

Opracowanie redakcyjne mgr Elbieta Winiarska-Grabosz Korekta Janina Szydowska amanie Alina Gebel Projekt okadki mgr in. Stanisaw Rogowski

Rozprawy CCLVII

© Copyright by Uniwersytet Przyrodniczy we Wrocawiu, Wrocaw 2009

ISSN 1897–4732 ISBN 978–83–60574–44–7

WYDAWNICTWO UNIWERSYTETU PRZYRODNICZEGO WE WROCAWIU Redaktor Naczelny – prof. dr hab. Andrzej Kotecki

ul. Sopocka 23, 50–344 Wrocaw, tel. 071 328–12–77 e-mail: wyd@up.wroc.pl

Nakad 100 + 16 egz. Ark wyd. 6,4. Ark. druk. 6,25 Druk i oprawa: Wydawnictwo Tekst Sp. z o.o.

SPIS TRECI

1. WSTP ... 7

2. CEL PRACY ... 9

3. ZASTOSOWANIE TECHNIKI TDR DO WYZNACZANIA WYBRANYCH DANYCH WEJCIOWYCH W RÓWNANIU RICHARDSA ... 11

3.1.WARUNKIPOCZTKOWE ... 14

3.2.WARUNKIBRZEGOWE ... 18

3.3.CZONRÓDOWY ... 22

3.4.PARAMETRYZACJAPRZESTRZENIOBJTEJ MODELOWANIEM ... 24

4. RÓWNANIE RICHARDSA Z RÓNYMI WARIANTAMI DANYCH WEJCIOWYCH ... 27

5. PRAKTYCZNE ZASTOSOWANIA TECHNIKI TDR ... 32

5.1.PRZYKADYWYZNACZANIAWARUNKÓW POCZTKOWYCH ... 32

5.2.SZACOWANIEGÓRNEGOWARUNKUBRZEGOWEGO ... 43

5.2.1. OPAD NETTO ... 43

5.2.2. PAROWANIE Z GLEBY NIEPORONITEJ ... 48

5.3.PRZESTRZENNYROZKADPOBORUWODYPRZEZ ROLINYUWZGLDNIANYJAKOCZONRÓDOWY ... 53

5.4.WYZNACZANIEPRZESTRZENNEJZMIENNOCI PARAMETRÓWGLEBY ... 57

6. WERYFIKACJIA I PORÓWNANIE CZTERECH WARIANTÓW MODELI MATEMATYCZNYCH ... 69

7. PODSUMOWANIE I WNIOSKI ... 84

7

1

WSTP

Racjonalne korzystanie z zasobów naturalnych wymaga rozpoznania praw rzdz-cych w rodowisku przyrodniczym. Umiejtne gospodarowanie wod w glebie jest moliwe, gdy znane s zasady, zgodnie z którymi nastpuje zmiana jej stanu energe-tycznego. Zasady te wyraone s za pomoc opisów matematycznych. Znaczenie sfor-muowa matematycznych i równoczenie pomiarów trafnie opisuj cytowane w pracy Malickiego (1999) sowa Sydenhama, który równie korzysta z myli lorda Kelvina: „Matematyczny opis wiedzy jest probierzem prawdziwoci nauki. Celem nauki jest opisanie wiedzy na podstawie danych wyraonych w kategoriach matematycznych. Ten punkt widzenia odzwierciedla si, midzy innymi, w stwierdzeniu lorda Kelvina: «Czsto powiadam, e kiedy moesz zmierzy to, o czym mówisz i kiedy moesz wyra-zi to liczbami, to ju wiesz cokolwiek o tym, ale jeeli nie moesz zmierzy tego, jee-li nie moesz wyrazi tego jee-liczbami, to twoja wiedza jest uboga i niewystarczajca: to moe by dopiero pocztek wiedzy, uczynie zaledwie may krok w kierunku nauki, bez wzgldu na to, co jest przedmiotem twoich rozwaa. »”

Modelowanie ruchu wody glebowej, której obecno okrela i wyznacza wszystkie wane funkcje gleby, jest wanym narzdziem do przewidywania jej stanu (Sobczuk 1998). Podstaw modelowania zjawisk przyrodniczych – w tym ruchu wody w orod-kach porowatych – stanowi równania róniczkowe. W celu prognozowania stosunków powietrzno-wodnych w glebie korzysta si najczciej z uogólnionego równania Richardsa (Brandyk 1986, Klinge i wsp. 2005, Komisarek, Kozowski 2005, Kowalik 2001, Olszta, Zaradny 2000, Reinhard 1992, Ross 1990) Aparat matematyczny wyko-rzystywany do roz-wizania równania metod rónic skoczonych lub elementów sko-czonych jest do dobrze rozpoznany (Belmans i wsp. 1983, Fipps i wsp. 1986, Kowal-ski 1998, MaciejewKowal-ski 1993, Reinhard 2004, SzulczewKowal-ski 2003, Wosiewicz 1989, Zaradny 1990). Jednak równanie to zawiera zawsze pewne elementy (np. wspóczynni-ki, czony niejednorodne), które charakteryzuj fizyczne waciwoci zjawiska i otocze-nia. S one wyznaczane za pomoc pomiarów, zawsze obarczonych bdem (Sobczyk 1991). Ponadto, w przypadku modelowania ruchu wody glebowej mamy do czynienia z orodkiem powstaym w wyniku procesów przyrodniczych. Dlatego wszystkie gleby charakteryzuj si zmiennoci przestrzenn waciwoci fizycznych (Brandyk i wsp.

8

1993). W zwizku z tym parametry fizyczne, które w tym przypadku reprezentuj mo-delowan przestrze, nie mog by wyraone przez jedn warto (Apul i wsp. 2005, Bykowski i wsp. 2001, Chalfen 1990, Janik 2001, 2008, Kajewski, Kowalski 1990, Li-piec 1983, Peck 1983, Pywaczyk, Kostrzewa 1985, Sobczuk 1998, Usowicz 2001, Walczykowa i wsp. 2005, Warrick, Nielsen 1980). W konsekwencji, do opisu filtracji najlepiej zastosowa teori pól losowych w postaci krigingu (Gnatowski i wsp. 1996, Isaaks, Srivastava 1998, Maciejewski 1998, Stpel 1993, Vaculin i wsp. 1983).

Proces prognozowania stosunków powietrzno-wodnych w glebie, oprócz modelo-wania, wymaga prowadzenia monitoringu wilgotnoci objtociowej. W ostatnich la-tach do tego celu stosowana jest najczciej technika reflektometryczna (TDR – Time Domain Reflectometry) (Malicki 1990, 1996, Malicki, Skierucha 1989, Reeves, Smith 1992, Roth i wsp. 1990, Skierucha 2000, 2005a, 2005b, Sobczuk, Plagge 2007, Topp, Davis 1985, Topp i wsp. 1980). W metodzie TDR bezporednio mierzy si czas propa-gacji impulsu elektromagnetycznego w badanym orodku na odcinku stanowicym podwójn dugo falowodu, który utworzony jest z prtów przewodzcych prd elek-tryczny. Na tej podstawie oblicza si prdko propagacji impulsu zwizan z przeni-kalnoci dielektryczn orodka, a ta jest zalena od wilgotnoci objtociowej badanego materiau. Najnowsze aparaty TDR skonstruowane s tak, aby dane uzyskane z pomia-rów mogy by przekazywane do uytkownika czem radiowym lub za pomoc sieci Internet (Skierucha i wsp. 2004).

Obecnie reflektometria domenowo-czasowa jest uznan na wiecie metod pomiaru wilgotnoci objtociowej gleby. Jednak, jeszcze 18 lat temu, Brandyk w swojej pracy stwierdza e: „Pomiar wilgotnoci aktualnej w warunkach terenowych pomimo istnienia kilku nowoczesnych metod jest cigle kopotliwy, czasochonny oraz wymagajcy dro-giej i skomplikowanej aparatury” (1990). W Polsce miernik TDR powsta w Instytucie Agrofizyki PAN w Lublinie. Do 2003 roku badaniami kierowa prof. Malicki, obecnie doc. Skierucha. Aparaty te stosowane s przez wikszo polskich orodków nauko-wych zajmujcych si oznaczeniami zawartoci wody w glebie. Dobre wyniki osiga si dla pomiarów w piaskach gliniastych lekkich i piaskach sabo gliniastych, w glinach lekkich i lekkich pylastych, a take w glinach cikich, bardzo cikich i cikich pyla-stych (Biniak 2004, yczko i wsp. 2000, Nyc, Pokadek 2004, Pczkowski i wsp. 2001, Pokadek 2001, Somorowska 2004, Szatyowicz i wsp.1994). Badania wykorzystujce technik TDR prowadzone s równie przez orodki badawcze w Europie i USA (Blonquist i wsp. 2005, Crow i wsp. 2005, Driksen, Dasberg 1993, Huisman i wsp. 2001, Ledieu i wsp. 1986, Nobiro 2001, Plagge i wsp. 1990, Thomsen i wsp. 2000, Topp i wsp. 1980, Topp, Reynolds 1998). Metoda reflektometryczna stosowana jest te z powodzeniem do pomiarów wilgotnoci w innych ni gleba orodkach porowatych, np. w materiaach budowlanych (Janik i wsp. 2006, Sobczuk i wsp. 2004).

9

2

CEL PRACY

Zamierzeniem autora bya poprawa dokadnoci prognozowania stosunków po-wietrzno-wodnych w wierzchniej warstwie gleby. Cel ten zrealizowano poprzez zasto-sowanie reflektometrii domenowo-czasowej (TDR) do bardziej precyzyjnego wyzna-czania danych wejciowych w modelach matematycznych opisujcych ruch wody w wierzchniej warstwie glebowej. W równaniu Richardsa wyznaczenie kompletu da-nych wejciowych sprowadza si do okrelenia: warunków pocztkowych, brzegowych oraz parametrów charakteryzujcych modelowan przestrze. W niektórych przypad-kach niezbdne jest równie wyznaczenie czonu ródowego, uwzgldniajcego np. pobór wody przez korzenie rolin. Badania polowe i eksperymenty laboratoryjne opisa-ne w rozdziale 5, mimo e zrónicowaopisa-ne tematycznie, cz dwa aspekty. Po pierwsze, kade z zagadnie dotyczy danych wejciowych w równaniu Richardsa, po drugie – w kadym zastosowano technik TDR. Te dwa punkty wspólne zadecydoway o wybo-rze materiau badawczego.

W niniejszej pracy warunek pocztkowy zosta sformuowany w czterech warian-tach. W pierwszym – po przyjciu pewnych zaoe upraszczajcych – równanie Richardsa sprowadzono do jednego wymiaru i warunek pocztkowy okrelono na pod-stawie pomiarów wilgotnoci aparatem TDR. W dwóch nastpnych wariantach równa-nie to rozpatrzono jako trójwymiarowe. W II – warunek pocztkowy okrelono, podob-nie jak w wariancie I, na podstawie pomiarów wilgotnoci w punktach wzowych. Na-tomiast w III wariancie wartoci wilgotnoci w punktach wzowych odczytano z map, w których interpolacj przeprowadzono metod krigingu. Kriging jest technik stoso-wan do lokalnej estymacji, w której tylko dane w pobliu obszaru estymacji s uwzgldniane w procesie szacowania. Technika ta stanowi metod interpolacyjn red-niej waonej ruchomej (Namysowska-Wilczyska 2006). W czwartym wariancie rów-nanie ponownie potraktowano jako jednowymiarowe, ale warunek pocztkowy w punktach wzowych realizowano jako zmienn losow. W kadym z wariantów technik TDR wykorzystano do wyznaczenia wilgotnoci objtociowej. Zastosowano j równie do szacowania warunków brzegowych – w szczególnoci do ustalenia gór-nego warunku brzegowego II rodzaju. Podano sposób szacowania opadu netto oraz pa-rowania z nieporonitej powierzchni gleby. Ponadto technik TDR wykorzystano do

10

wyznaczania czonu ródowego w równaniu Richardsa, co szczegóowo opisano w rozdziale 5.3. Dokadno zaproponowanych poniej zastosowa techniki TDR prze-analizowano na podstawie bada laboratoryjnych i polowych. Wyznaczanie parametrów charakteryzujcych modelowan przestrze i stosowanych w równaniu Richardsa, po-dobnie jak warunek pocztkowy, przeprowadzono w czterech wariantach. W wariancie I zaoono, e parametry s jednakowe w kadym punkcie modelowanej przestrzeni a ich wartoci, dla danego typu gleby, okrelono na podstawie danych literaturowych (Genuchten van 1991). Ustalenia musz by poprzedzone badaniami pozwalajcymi oznaczy skad granulometryczny materiau glebowego. Drugi sposób (wariant II) wy-korzystuje technik TDR i polega na przeprowadzeniu tzw. kalibracji modelu przepy-wu wody dla równania Richardsa. Pozwolio to na wyznaczenie wartoci parametrów glebowych w wielu pionowych przekrojach modelowanej przestrzeni. Kalibracj mode-lu w niniejszej pracy nazywano równie procedur identyfikacyjn. Trzeci i czwarty sposób to, podobnie jak w przypadku warunku pocztkowego, podanie wartoci para-metrów w postaci map (wariant III) lub w postaci zmiennej losowej (wariant IV).

W rozdziale 6 przedstawiono przykady zastosowania modelu matematycznego do prognozowania stosunków powietrzno-wodnych w wierzchniej warstwie gleby uytku kowego. W tym celu równanie Richardsa zostao rozwizane z danymi wejciowymi zaproponowanymi w niniejszej pracy. Powstay w ten sposób cztery warianty modelu pozwalajcego prognozowa stosunki powietrzno-wodne w wierzchniej warstwie gleby. Przydatno modeli przeanalizowano i zweryfikowano, wykorzystujc niezaleny mate-ria uzyskany na podstawie bada polowych na uytku kowym pooonym w miej-scowoci Brenna w województwie lskim.

11

3

ZASTOSOWANIE TECHNIKI TDR

DO WYZNACZANIA WYBRANYCH

DANYCH WEJCIOWYCH

W RÓWNANIU RICHARDSA

Pomiar wilgotnoci objtociowej gleby technik TDR odróniaj od tradycyjnej metody grawimetrycznej trzy cechy. Po pierwsze, pomiar charakteryzuje maa praco-chonno. Umoliwia to przeprowadzenie oznacze w wielu punktach, w krótkim cza-sie, i w konsekwencji rozpoznanie przestrzennej zmiennoci wilgotnoci w modelowa-nej przestrzeni. Dziki temu mona precyzyjnie wyznaczy warunek pocztkowy – np. w postaci map wilgotnoci lub te moe on by realizowany jako zmienna losowa. Po drugie, pomiar wykonywany technik TDR jest praktycznie nieinwazyjny. W zwizku z tym, oznaczenia wilgotnoci mog by prowadzone wielokrotnie w tym samym punk-cie. Po trzecie, konstrukcja sterujca prac miernika, którego budow opisano m.in. w pracy Wilczka i Skieruchy (Wilczek, Skierucha 2007) (chodzi tu o zastosowanie dwupoziomowego przeczania, multipleksowania czujników), umoliwia cigy moni-toring mierzonej wielkoci – np. co jedn minut w wielu punktach jednoczenie. Ce-chy te pozwalaj na budow bilansu wodnego w elementach przestrzeni powstaych po

jej dyskretyzacji (Vi,j,k) (rys. 1). Jest to pomocne, przy pewnych zaoeniach, w

oblicza-niu zmian objtoci wody w tych elementach dla dowolnie krótkich kroków czasowych. Na tej podstawie wyznaczono warunek brzegowy II rodzaju oraz czon ródowy. Ponadto, specyfika pomiaru wilgotnoci gleby aparatami TDR pozwala zastosowa procedury identyfikacyjne, których efektem jest wyznaczenie przestrzennego zrónico-wania parametrów modelowanego orodka.

12

k+1

Rys. 1. Dyskretyzacja modelowanej przestrzeni

Vi,j,k – element przestrzeni powstay w wyniku dyskretyzacji; i,j,k – cakowita wysoko

cinienia w elemencie Vi,j,k; Ti,j,k – wilgotno w elemencie Vi,j,k; i ,j ,k – indeksy

odlego-ci w kierunkach x, y, z; L, M, N, – liczba elementów w kierunkach x, y, z Fig. 1. Discretisation of modelled space

Vi,j,k – element of space formed as a result of discretisation; i,j,k – total level of pressure

in element; Vi,j,k, Ti,j,k – moisture in element Vi,j,k; i ,j ,k – indexes of distance in directions

x, y, z; L, M, N, – number of elements in directions x, y, z

W niniejszej pracy w modelu prognozujcym stosunki powietrzno-wodne zastoso-wano równanie Richardsa, które w oryginalnej postaci opisuje ruch wody jedynie w strefie nienasyconej (Richards 1931). Dla przestrzeni trójwymiarowej, przedstawio-nej na rysunku 1, przyjmuje ono, z dodatkowym czonem ródowym, nastpujc po-sta (Brandyk 1990, Reinhard 2001):

S z ) h ( K z y ) h ( K y x ) h ( K x t ) h ( C  »¼ º «¬ ª w ) w w w  » ¼ º « ¬ ª w ) w w w  »¼ º «¬ ª w ) w w w w ) w ˜ , (1)

gdzie: C(h) – róniczkowa pojemno wodna, dh d

C T , cm-1, K(h) – przewodno hydrauliczna, cm min-1,

)– cakowita wysoko cinienia, ) hz, cm H2O,

h – wysoko cinienia, cm H2O,

T – wilgotno, m3

m-3,

x, y, z – wspórzdne przestrzeni objtej modelowaniem, cm, t – czas, min,

S – czon ródowy, min-1.

Vi,j,k k = 1,2,3,…,N z y )i,j,k i = 1,2,3,…,L _x j = 1,2,3,…,M

kolumna Vi, j, column Vi,j

k-1 k Ti, j, k

13

Natomiast rónicow posta tego równania z czonem ródowym (schemat jawny) zapisano nastpujco (Reinhard T., Reinhard A. 2005):





 

V , Q z K z K y K y K x K x K t C k , j , i t k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i k , j , i 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 1  '              ˜ '  ' )  ) ˜  ' )  ) ˜  ' )  ) ˜  ' )  ) ˜  ' )  ) ˜  ' )  ) ˜ )  ) ˜ ' W W W W W W W W W W W W W W W W W W W W W (2) gdzie: W k , j , i

C – róniczkowa pojemno wodna dla elementu przestrzeni powstaym po jej dys-

kretyzacji (w dalszej czci pracy nazwany elementem Vi,j,k), cm-1,

W k , j , i

) – cakowita wysoko cinienia dla elementu Vi,,j,k, cm H2O,

W 1  k , j , i

K – przewodno hydrauliczna dla elementu Vi,j,k, cm min-1,

't – krok czasowy, min,

'x, 'y, 'z – krok przestrzenny w kierunku x, y, z, cm,

t k , j , i Q'

' – zmiana objtoci wody w elemencie Vi,,j,k w czasie 't, cm3 min-1,

Vi,,j,k – objto elementu ze wskanikami i,j,k, cm3,

 – indeks czasu, –,

i, j, k – indeksy w kierunkach x, y, z, –.

Do wyznaczenia krzywej przewodnoci hydraulicznej zastosowano formu (Genuchten van 1991): 2 ] ) | | ( 1 [ } ] ) ( 1 [ ) ( 1 { ) ( 2    ˜   m n m n mn s h h h K h K D D D , (3)

i zaleno odpowiadajc krzywej pF wyznaczono ze wzoru:

m n r s r h h ] |) | ( 1 [ ) ( D T T T T    , (4) gdzie:

m, n,  – wspóczynniki zalene od typu gleby, m =1 – n–1, –,

Tr – zawarto wody residualnej, m3m-3,

Ts – wilgotno w strefie penego nasycenia, m3m-3,

Ks – wspóczynnik filtracji, cm min-1,

h – wysoko cinienia, cm H2O,

14

Zalenoci (2), (3) i (4) wykorzystano do budowy modelu matematycznego pozwa-lajcego prognozowa stosunki powietrzno-wodne. Do rozwizania równania (2) dla kolejnych chwil czasowych wykorzystano procedury obliczeniowe napisane przez auto-ra w progauto-ramie Excel oauto-raz procedury w progauto-ramie Matlab napisane przez dr. in. ukasza Balbusa z Politechniki Wrocawskiej z Instytutu Matematyki i Informatyki. Poprawno procedur sprawdzono poprzez porównanie wyników symulacji z wynikami symulacji zawartymi w pracy prof. Reinharda (2004).

3.1. WARUNKI POCZTKOWE

Opisane zostan teraz cztery warianty formuowania warunków pocztkowych. W pierwszym przyjto, e równanie (2) jest jednowymiarowe. W drugim i trzecim – model potraktowany jest jako trójwymiarowy, za w czwartym wariancie warunek po-cztkowy realizowano jako zmienn losow, a równanie ponownie potraktowano jako jednowymiarowe.

Zaómy na pocztek, e modelowan przestrzeni jest równomiernie poronity uytek kowy o takim samym ksztacie jak przestrze przedstawiona wczeniej na rysunku 1. Jednoznaczno rozwizania rónicowej postaci równania Richardsa (zale-no (2)) wymaga okrelenia warunku pocztkowego. W wariancie I, najbardziej upraszczajcym zagadnienie, przyjto, e ruch wody odbywa si wycznie w kierunku pionowym oraz warunki ruchu s jednakowe w kadym przekroju. W zwizku z tym, w modelu caa przestrze moe by reprezentowana wycznie przez wybran kolumn

Vi,j, któr przedstawiono na rysunku 1. Zaoenie to jest zgodne z rzeczywistoci, gdy

teren jest paski, zasilanie dolnej warstwy jest jednakowe w kadym elemencie Vi,j,N,

jednakowe s równie warunki otoczenia wpywajce na ruch wody w górnej warstwie przestrzeni. Po spenieniu tych zaoe wilgotno objtociowa w kadej warstwie

osobno moe by reprezentowana wycznie przez jedn warto (Ti,j,k dla staego k jest

jednakowa). Równanie Richardsa rozpatrywane jest wtedy wycznie w kierunku osi z.

W wariancie I do wzoru (2) jako warunek pocztkowy naley wprowadzi wartoci

potencjau ) w chwili pocztkowej, przy czym wskaniki i,j mog by tu pominite. _kpo

Wartoci te uzyskano na podstawie wzoru (5). W wariancie II jako warunek

pocztko-wy podawana jest warto potencjau )_ipo_,j,k w chwili pocztkowej (po) w kadym

ele-mencie Vi,j,k,. Wielkoci )_ipo_,j,k obliczono z zalenoci:

k , j , i po k , j , i po k , j , i h z ) , (5) gdzie: po k , j , i

) – potencja cakowity wody dla elementu Vi,j,k w chwili pocztkowej, cm H2O,

po k , j , i

h – wysoko cinienia w elemencie Vi,j,k w chwili pocztkowej, cm H2O,

15

Rys. 2. Wybrana para punktów w k-tej warstwie, dla której obliczona jest semiwariancja

i,j,k – wilgotno w elemencie Vi,j,k; i+li, j+lj, k – wilgotno w punkcie odlegym o l;

pozostae oznaczenia jak w rys. 1

Fig. 2. Selected pair of points in k-th layer, for which semivariance is calculated

i,j,k – moisture in element Vi,j,k; i+li, j+lj, k – moisture at a point distant by l; other

desig-nations as in Fig. 1

Wielko z we wzorze (5) jest zalena od gbokoci pooenia rozpatrywanego

elementu Vi,j,k.Wielko hi,j,k obliczono po odpowiednim przeksztaceniu zalenoci (4).

W zwizku z tym, wymagane jest zmierzenie wilgotnoci objtociowej T_i,j,kw

ka-dym elemencie Vi,j,k. Wtedy zmniejszanie kroków przestrzennych 'x, 'y wie si

z przeprowadzeniem wielu pracochonnych pomiarów. Przykadowo, chcc modelowa

ruch wody na obszarze o wymiarach 20 m x 20 m z krokiem 'x = 0,2 m i 'y = 0,2 m,

naleaoby tylko w jednej warstwie zmierzy wilgotno w 100 x 100 = 10 000 punk-tach. Ponadto, stosujc ten wariant, nie ma moliwoci wyeliminowania tych oznacze

0 t k , j , i

T aparatem TDR, które wynikaj np. z niejednorodnoci orodka (Skierucha i wsp.

2004, Wasilewski i wsp. 2005).

Warunek pocztkowy w wariancie III formuowany jest podobnie jak w wariancie

II, z t rónic, e wielkoci T_it_,j0_,k odczytywane s z map wilgotnoci zbudowanych dla

kadej k-tej (k = 1,2,3,...,N) warstwy (rys. 2). Mapy mog by zbudowane na podstawie

mniejszej liczby punktów pomiarowych. Ponadto, procedury interpolacyjne eliminuj bdy pomiarów oraz pozwalaj na odczytanie interpolowanej wartoci wilgotnoci w kadym punkcie modelowanego obszaru – równie w tym, w którym nie dokonywa-no pomiarów. Wymaga to jednak przeprowadzenia analizy obszarowej zmiendokonywa-noci

wil-gotnoci objtociowej w kadej k-tej warstwie. Stosuje si w tym celu metody

geosta-k = 1,2,3,…,N z x y j = 1,2,3,…,M i = 1,2,3,…,L Ti,j,k Ti+li, j+lj, k lj,k li,k lk

16

tystyczne (Cahil i wsp. 1999, Goldstein, Skrzypek 2004, Grenhalz i wsp. 1988, Grego i wsp. 2006, Mohanty i wsp. 2000, Nielsen 1973, Nyberg 1996, Usowicz B., Usowicz . 1999). W pierwszej kolejnoci budowane s wariogramy, których parametry pozwa-laj na optymalizacj oblicze interpolacyjnych metod krigingu.

Wariogramy dla jednej paszczyzny (np. warstwy k) definiuje si jako wariancj

wilgotnoci w funkcji odlegoci porównywanych ze sob punktów (Dbowska, Za-wadzki 2005, Kossowski, Usowicz 2000, Mc Bratney, Webster 1986, Usowicz 1999, Wang i wsp. 2001):





2 ) ( 1 , , ) ( 2 1 ) (

¦

 _ ˜ k k l m p k l p k p k k k l m l T T J , (6) gdzie:

Jk(lk) – wariogram dla k-tej warstwy, (m3m-3)2,

m(lk) – liczba porównywanych par punktów odlegych od siebie o dystans lk, (lk – rys. 2),

p,k=i,j,k – wilgotno w elemencie Vi,j,k, m3m-3,

p+lk,k = i+li, j+l,j k – wilgotno w elemencie odlegym od elementu Vi,j,k o dugo lk,

(rys. 2), m3m-3.

Podczas budowy semiwariogramów naley zwraca uwag na ich stabilno zalen od liczby porównywanych par punktów (Janik, Miciorak 2007). W niektórych przy-padkach wartoci semiwariogramów obliczone na podstawie zalenoci (6) s takie, e dobranie do nich modelu matematycznego, wykorzystywanego podczas interpolacji, jest niemoliwe. Wtedy naley podj prób przeprowadzenia analizy semiwariogra-mów standaryzowanych obliczonych z zalenoci (Keim i wsp. 2005, Usowicz i wsp. 2004):

 



k

  

k k k k k k ks l l l l G G J J ˜ 0 ) ( , (7) gdzie: ) (k ks l

J – semiwariogram standaryzowany zbudowany dla k-tej warstwy, (–)2,

 

k

k l

J – semiwariogram empiryczny zbudowany dla k-tej warstwy, (m3m-3)2,



k 0



k l

G – odchylenie standardowe wartoci zmiennej losowej, dla k-tej warstwy,

w punkcie zaczepienia wektora lk, m3m-3,

 

k

k l

G – odchylenie standardowe wartoci zmiennej losowej, dla k-tej warstwy

w punkcie oddalonym od punktu zaczepienia wektora o odlego lk, m3m-3.

Zbudowane semiwariogramy pozwalaj na zastosowanie fraktali. Warto wymiaru fraktalnego okrela, na ile badana cecha jest zdeterminowana, a na ile jej rozkad ma charakter losowy (Anderson i wsp. 1998, Gimenez i wsp. 1997, Sawiski i wsp. 2002, Usowicz 1999). Semiwariogramy (klasyczne bd standaryzowane) wykorzystano na-stpnie podczas interpolacji do budowy map. Na kocowym etapie, w kadej warstwie oddzielnie, przeprowadzono analiz zgodnoci danych uzyskanych z map i danych zmierzonych. Zastosowano wymiennie nastpujce miary: maksymaln bezwzgldn

17

rónic pomidzy wartociami zmierzonymi i obliczonymi (R|max|) (zaleno (8)),

red-ni warto bezwzgldn rónic wartoci zmierzonych i obliczonych (SWB) (zaleno

(9)), pierwiastek redniej kwadratowej rónic wartoci zmierzonych i obliczonych (R.M.S.) (zaleno (10)). Ponadto, w niektórych przypadkach, maksymaln rónic

pomidzy wartociami zmierzonymi i obliczonymi Pmax (zaleno (11)), minimaln

rónic tych samych wielkoci Pmin (zaleno (12)) oraz rozstp Rs i wspóczynnik

zmiennoci Z. Poniej podano definicje wybranych miar dla k = const. (Barnet 1982):

^

obl

`

k j i pom k j i j i k R _{, , , ,} , max max T T , (8)

¦¦

 ˜ L i M j obl k , j , i pom k , j , i k WB M L S 1 1 1 T T , (9) M L S M R L i M j obl k j i pom k j i k ˜ 

¦¦

1 2 1 , , , , ) ( . .

T

T

, (10)

^

obl

`

k j i pom k j i j i k P _{,, , ,} , max maxT T , (11)

^

obl

`

k j i pom k j i j i k P _{, , ,,} , min minT T , (12) gdzie: pom k , j , i

T – zmierzona warto wilgotnoci w elemencie Vi,j,k, m3m-3,

obl k , j , i

T – odczytana z mapy warto wilgotnoci dla elementu Vi,j,k, m3m-3,

k

max

R – maksymalna bezwzgldna rónica pomidzy wartociami zmierzonymi i odczy-

tanymi z map, m3m-3,

k

WB

S – rednia warto bezwzgldna rónic wartoci zmierzonych i obliczonych, m3m-3,

R.M.Sk – redni bd kwadratowy, m3m-3,

k max

P – maksymalna rónica pomidzy wartociami zmierzonymi i obliczonymi, m3m-3,

k

P

_min – minimalna rónica pomidzy wartociami zmierzonymi i obliczonymi, m3m-3,

L – liczba elementów w kierunku osi x, –, M – liczba elementów w kierunku osi y, –.

Ostatecznie, przy sformuowaniu warunku zgodnie z wariantem III modelowana

przestrze moe by dyskretyzowana na dowolnie mae objtoci Vi,j,k i rónicowa

posta równania Richardsa (zaleno (2)) moe by rozwizana z dowolnie maym krokiem przestrzennym. Jedynym ograniczeniem jest dugo prowadzonych oblicze przez komputer.

18

Równanie Richardsa, jak napisano powyej, zawiera elementy (np. wspóczynniki, czony ródowe, równie warunki pocztkowe i brzegowe) wyznaczane na podstawie pomiarów. Na skutek bdów pomiarowych nie powinny by one wyraane przez jedn warto. W zwizku z tym, bardziej realistycznym modelem ruchu wody w glebie b-dzie model wykorzystujcy stochastyczne równanie róniczkowe, to znaczy równanie, które zawiera elementy losowe (Janik 2007b, Maciejewski i wsp. 1994, Russo 1986a, 1986b, 1994, Sobczyk 1991). Rozróniane s równania z losowymi warunkami poczt-kowymi, równania z losowym czonem niejednorodnym i równania z losowymi wspó-czynnikami. W niniejszej pracy przeanalizowano równanie z losowym warunkiem pocztkowym i losowymi parametrami.

Obecnie rozpatrzony zostanie przykad warunku pocztkowego realizowanego jako zmienna losowa. Bdzie to IV wariant sformuowania warunku pocztkowego. W opi-sanym powyej wariancie I przyjto, e ruch wody odbywa si wycznie w kierunku pionowym. Jak zaznaczono, wymaga to przyjcia zaoenia, e wilgotno

objtocio-wa oznaczana w wielu punktach k-tej warstwy jest jednakowa. W warunkach

natural-nych zaoenie to nie jest spenione. Wynika to ze zrónicowania fizycznatural-nych waciwo-ci gleby w punktach, które s pooone od siebie w niewielkiej odlegowaciwo-ci (Janik 2005a, Kostrzewa, Pulikowski 1993, Rigby, Porporato 2006, Webster, Oliwer 1990).

W zwizku z powyszym, podanie cile okrelonej wartoci wilgotnoci i,j,k nie jest

reprezentatywne dla caej warstwy, co moe by powodem niezadowalajcych wyników

symulacji. Fakt zrónicowania wilgotnoci w k-tej warstwie moe by uwzgldniony

poprzez potraktowanie i,j,k jako zmiennej losowej. Równanie Richardsa potraktowane

jest wtedy jako jedna z klas stochastycznych równa róniczkowych – w tym przypadku jest to równanie z losowym warunkiem pocztkowym. Sposób uwzgldnienia losowego charakteru parametrów wynikajcego z waciwoci warstwy opisano w rozdziale 3.4.

3.2. WARUNKI BRZEGOWE

Rozpatrzone zostan teraz przykady zastosowania techniki TDR do wyznaczania górnego warunku brzegowego. Na rysunku 3 przedstawiono fragment paskiego uytku kowego, w którym naley wyznaczy górny warunek brzegowy II rodzaju – tzw. wa-runek Neumana. Wawa-runek ten wystpuje wtedy, gdy na górnym brzegu modelowanej przestrzeni znane s wartoci pochodnej normalnej do brzegu (Kowalski 1998). Ozna-cza to konieczno okrelenia natenia przepywu wody przez górn powierzchni, na

który wpywaj: opad netto Pi,j,1, ewaporacja Ei,j,1 oraz transpiracja z wierzchniej

war-stwy Ti,j,1. W rozpatrywanym tu przypadku przyjto, e dolna warstwa jest

nieprzepusz-czalna. Zmiana objtoci wody 'Qi,j,1 w górnej warstwie moe by wywoana

wielko-ciami Pi,j,1, Ei,j,1, Ti,j,1 oraz dopywem wody z warstwy poniej 21

, , j , i K lub wypywem do warstwy poniej K_i_,,_j2_,1.

19

Rys. 3. Czynniki wpywajce na górny warunek brzegowy II rodzaju

Pi,j,1 – opad netto; Ti,j,1 – transpiracja rolin z górnej warstwy; Ei,j,1 – ewaporacja z

po-wierzchni gleby; 2 1 , , j , i

K – objto wody, która dopyna do elementu Vi,j,1 z elementu

Vi,j,2; 21 , , j , i

K – objto wody, która odpyna z elementu Vi,j,1 do elementu Vi,j,2; Ki k,,j,k1 –

objto wody, która dopyna do ele-mentu Vi,j,k z elementu Vi,j,k-1; Ki k,,j,k1 – objto

wody, która dopyna do elementu Vi,j,k z elementu Vi,j,k+1; K_i k_,,_j,k1– objto wody, która

odpyna z elementu Vi,j,k do elementu Vi,j,k-1; K_i k_,,_j,k1 – objto wody, która odpyna

z elementu Vi,j,k do elementu Vi,j,k+1; Ki N,,j,N1 – objto wody, która dopyna do N-tego

elementu z elementu N-1;  N, 1 N , j , i

K – objto wody, która odpyna z N-tego elementu do

elementu N-1

Fig. 3. Factors that affect the upper boundary condition of the 2nd kind

Pi,j,1 – net precipitation; Ti,j,1 – transpiration of plants from upper horizon; Ei,j,1 –

evapora-tion from soil surface; 2

1 , , j , i

K – volume of water that migrated to element Vi,j,1 from

ele-ment Vi,j,2, K_i_,,_j2_,1 – volume of water that migrated from element Vi,j,1 to element Vi,j,2;

1   k, k , j , i

K – volume of water that migrated to element Vi,j,k from element Vi,j,k-1, K_i k_,,_j,k1 –

vol-ume of water that migrated to element Vi,j,k from element Vi,j,k+1; Ki k,,j,k1– volume of water

that migrated from element Vi,j,k to element Vi,j,k-1; Ki k,,j,k1 – volume of water that migrated

from element Vi,j,k to element Vi,j,k+1; K_i N_,,_j,N1 – volume of water that migrated to N-th

ele-ment from eleele-ment N-1;  N, 1

N , j , i

K – volume of water that migrated from N-th element to

element N-1 z x y Ti,j,1 Pi,j,1 Ei,j,1 K+,2 i,j,1 K-,2 i,j,1 K+,k+1 i,j,k K+,k-1 i,j,k K-,k+1 i,j,,k K-,k-1 i,j,k K+, N-1 i,j,N K -,N-1 i,j,N

20

Przyjmijmy dodatkowo, e Pi,j,1, Ei,j,1 oraz Ti,j,1, a take Ki,j,2,1 i Ki,,j2,1 s jednakowe

dla kadego wskanika i, j. Dlatego w dalszej czci tego rozdziau wskaniki te zostay

pominite. Równanie bilansu wodnego dla elementu V1 zapisano nastpujco:

1 2 1 1 1 1 2 1 1 1 1˜'    ˜'      ˜'  'Q t P K , t T E K , t , (13) gdzie:

'Q1 – zmiana objtoci wody w elemencie V1, cm3,

P1 – opad netto, cm3min-1,

T1 – transpiracja rolin z górnego elementu, cm3min-1,

E1 – ewaporacja z powierzchni (skadniki P1 oraz E1 wystpuj wymiennie), cm3min-1,

2 1

,

K – objto wody, która dopyna do elementu V1 z elementu V2, cm3,

2 1

,

K – objto wody, która odpyna z elementu V1 do elementu V2, cm3,

t

'

– krok czasowy, dla którego wyznaczana jest wielko 'Q1, min.

Dla k >1 (poza warstw doln) mona zapisa:





k k , k k , k k , k k , k k t K K K K t T Q ˜'    ˜'  ' 1  1  1  1  1 1 , (14) gdzie:

Tk – transpiracja rolin z elementu Vk, cm3min-1,

pozostae oznaczenia jak na rysunku 3 i w zalenoci (13)

oraz dla k = N, mona zapisa:

N N , N N t K t T Q ˜' ˜'  ' 1  1 1 _lub N N , N N t K T Q ˜'  ' 1  1 _{, (15)} gdzie:

TN – transpiracja z elementu VN, cm3min-1,

pozostae oznaczenia jak na rysunku 3 i w zalenoci (13).

Do wyznaczenia wielkoci Qk zastosowano technik TDR (Janik 2004). Na

rysun-ku 4 zaprezentowano dowolny element gleby Vk, w której znajduje si faza staa,

powie-trze i woda. W lewej czci rysunku przedstawiono stan w chwili pocztkowej tpo, za

w prawej, po czasie 't, w chwili kocowej tko. Jeeli w chwili pocztkowej wilgotno

objtociowa w dowolnym elemencie wynosi T a w chwili kocowej _ipo T , to zgodnie _kko

z rysunku 4 mona wyznaczy zmian objtoci wody w elemencie Vk w czasie t wg

zalenoci: 1  ' _{ ˜'} 'Q_kt V_k(T_kko T_kpo) t , (16) gdzie: t k Q'

21

k

V – objto elementu Vk,cm3,

po k

T – pocztkowa wilgotno objtociowa w elemencie Vk, m3m-3,

ko k

T – kocowa wilgotno objtociowa w elemencie Vk, m3m-3,

t

' – krok czasowy, min.

W powyszych rozwaaniach naley dodatkowo przyj, e wilgotno w kadym

punkcie elementu Vk, w danym kroku czasowym, jest jednakowa. Wtedy tylko wielko

t k

Q'

' wyznaczona zostanie poprawnie.

Rys. 4. Zmiana wilgotnoci objtociowej jako informacja o zmianie objtoci wody w elemencie Vk

tpo – chwila pocztkowa; tko – chwila kocowa; t – krok czasowy; Vk – objto

elemen-tu; Vgk – objto fazy gazowej; Vwk – objto wody; Vsk – objto fazy staej; k –

wil-gotno objtociowa; po – indeks dla chwili pocztkowej; ko – indeks dla chwili kocowej

Fig. 4. Change in volumetric moisture content as information on change of water volume in element Vk

tpo – initial moment; tko – final moment; t – time step; Vk – volume of element; Vgk –

volume of gaseous phase; Vwk – volume of water; Vsk – volume of solid phase; k –

volu-metric moisture content; po – index for initial moment; ko – index for final moment

W literaturze mona znale wiele metod pozwalajcych wyznaczy opad netto, transpiracj rolin oraz ewaporacj (Bac, Rojek 1990, Banach i wsp. 1983, Bry K., Bry T. 2001, Jaworski 1990, Kowanetz 2004, Musia 2001, Musia, Rojek 1990a). Ich zastosowanie wymaga jednak okrelenia, podczas bada laboratoryjnych bd polo-wych, wielu empirycznych wspóczynników, z których kady obarczony jest bdem

(Zaradny 1990). Poniej przedstawiono metod pozwalajc szacowa wielkoci P1, TK,

E1 bez koniecznoci wyznaczania jakichkolwiek empirycznych wspóczynników (Janik

2004, 2006). Jest to metoda porednia wykorzystujca ide rozwizania zadania od-wrotnego. Polega ono na pomiarze wielkoci, która jest skutkiem ocenianego zjawiska.

Przykadowo, chcc w czasie 't wyznaczy opad netto P1, dokonywany jest pomiar

wilgotnoci objtociowej k w elementach Vk w chwili pocztkowej i kocowej. Biorc

pod uwag specyfik pomiaru aparatem TDR, mona to uczyni wielokrotnie w tym samym punkcie. Powinien by przy tym speniony warunek selektywnoci, tzn. skutek zjawiska (zmiany wilgotnoci objtociowej) powinien by wywoany wycznie

Vwkok Vgkok Vskok Vw pok Vs pok Vg pok z Vgkok Vw po k z Vwkok Vspok = Vs ko k Vpok Tpo k zTkok Vg pok tpo + _'t = tko Vkok

22

opadem netto P1. Warunek selektywnoci musi by speniony równie wtedy, gdy

chcemy obliczy ewaporacj lub transpiracj. W rzeczywistoci, gdy gleba jest poro-nita, warunki te nie s spenione. W zwizku z tym, do zastosowania metody

pozwa-lajcej obliczy oddzielnie wielkoci P1, T1 i E1 musz by poczynione pewne zaoenia.

Jeeli

z równania (13) chcemy wyznaczy wielko P1, to naley zagwarantowa, e TK = 0

dla kadego k od 1 do N. W wikszoci przypadków podczas opadu

¦

N k k T 1 << P1,

w zwizku z tym w przyblieniu mona przyj, e TK | 0.Ponadto, z definicji opadu

netto wynika, e gdy wystpuje skadnik P1, to skadnik E1 naley pomin, pomimo to

e wystpuje ewaporacja. Ostatecznie jedyn niewiadom równania (13) pozostaje

wielko P1. Jeeli natomiast chcemy obliczy wielko E1, naley zagwarantowa, e

w danym okresie P1 = 0 i gleba jest nieporonita (T1 = 0). Podobnie, jeeli wyznaczamy

T1, ustalamy, e P1 = 0 i E1 = 0 (brak parowania z powierzchni gleby). Zaoenie o braku

parowania z powierzchni gleby mona speni, np. eliminujc kontakt gleby z powie-trzem poprzez zastosowanie nieprzepuszczalnej warstwy. Naley podkreli, e

wielko-ci TK zostan obliczone poprawnie jedynie wtedy, gdy przyjmiemy, e ruch wody

w obrbie masy korzeniowej wynika gównie z poboru wody przez rolin (woda

po-midzy warstwami nie przemieszcza si w czasie t). Jeeli warunek ten nie zostaby

speniony, to wielko

¦

N k k T 1

i tak zostaaby wyznaczona poprawnie. Podsumowujc powysze rozwaania, stwierdzamy, e w pewnych przypadkach znajc jedynie

dy-namik wilgotnoci k w elementach Vk, mona wyznaczy ewaporacj E1,

transpira-cj TK oraz opad netto P1. Wielkoci te s niezbdne do podania warunku brzegowego

II rodzaju.

3.3.

CZON

RÓDOWY

Zmiana wilgotnoci w modelowanej przestrzeni moe nastpi w wyniku np. reakcji chemicznych przebiegajcych z wydzieleniem lub wizaniem wody, ale w pracy si tego nie rozwaa. Poza tym zmiana wilgotnoci moe zaj na skutek dziaania innego mechanizmu ni opisywany równaniem Richardsa (np. poprzez stosowanie nawodnie wgbnych lub transport wody przez korzenie rolin (Jeznach 1996)). W modelu mate-matycznym najlepiej uwzgldni to poprzez czon ródowy w rónicowej postaci rów-nania Richardsa – równanie (2). Wymagane jest wtedy okrelenie objtoci podawanej

wody Qti,j,k do kadego wewntrznego elementu Vi,j,k (rys. 3). Wielko Qti,j,k

w przypadku wgbnych systemów nawadniajcych jest atwa do wyznaczenia, ponie-wa to urzdzenia techniczne decyduj o miejscu i iloci podawanej wody. Jeeli nato-miast czon ródowy ma reprezentowa pobór wody przez roliny, to zagadnienie jest zoone (Janik 2005b). Pojawia si wtedy konieczno wyznaczenia intensywnoci po-bieranej wody – najlepiej w kadym elemencie dyskretyzowanej przestrzeni, w której

23

wystpuje masa korzeniowa rolin. Intensywno poboru wody jest zalena od gstoci masy korzeniowej, stadium rozwojowego roliny, zmiennych w cigu doby warunków atmosferycznych oraz stopnia zagszczenia gleby (Lipiec i wsp. 1998, Mioduszewski 2006, Nosalewicz, Lipiec 2002). W literaturze mona znale wiele formu matema-tycznych do wyznaczania intensywnoci poboru wody przez roliny, tzw. transpiracji aktualnej (Musia, Rojek 1990a, 1990b, 1999, Rojek 1987, Szulczewski 1998, yromski 1990). Zoono procesu sprawia, e jego ilociowy opis wymaga wyznaczenia danych wejciowych. Przykadowo, w modelu zaproponowanym przez Fedesa i wsp. (1978) tylko do wyznaczenia ewapotranspiracji potencjalnej naley poda: strumie energii pochodzcy od radiacji netto, strumie ciepa odprowadzany z powierzchni do atmosfe-ry w wyniku przewodzenia i konwekcji, strumie ciepa zuyty na parowanie, ciepo utajone parowania, strumie pary oraz strumie ciepa skierowany do gleby. Ponadto, w modelach tych masa korzeniowa roliny potraktowana jest jako cao, w zwizku

z tym nie ma moliwoci podania przestrzennego rozkadu funkcji ródowej

–

podob-nie jak w metodach ewaporometrycznych. Istpodob-niej rówpodob-nie metody pozwalajce szaco-wa przestrzenne zrónicowanie poboru wody przez rolin, wykorzystujc mudne makroskopowe pomiary gstoci masy korzeniowej roliny. W niniejszym rozdziale opisano metod zastosowania techniki TDR do wyznaczenia aktualnej intensywnoci poboru wody przez korzenie rolin. Zalet tej metody jest to, e nie ma koniecznoci opisu adnych empirycznych wspóczynników. Na rysunku 5 przedstawiono dwa spo-soby uwzgldniania w modelach poboru wody przez korzenie rolin. Po stronie lewej,

gdy masa korzeniowa znajduje si wycznie w brzegowych elementach Vi,j,1,

stosowa-ny jest najczciej omawiastosowa-ny w rozdziale 3.2 górstosowa-ny warunek brzegowy II rodzaju. Gdy jednak korzenie rolin znajduj si poza warstw brzegow, to naley to uwzgldni poprzez wprowadzenie czonu ródowego. Kolorem zielonym zaznaczono te elementy poza brzegami, w których wystpuje pobór wody przez rolin. Jeeli ponownie tak jak

w rozdziale 3.3 zaoymy, e w czasie 't woda nie przemieszcza si midzy

elementa-mi Vi,j,k , to korzystajc z zalenoci (13) rozszerzonej o indeksy i,j, mona obliczy dla

czasu 't zmian objtoci wody 'Q_i'_,t_j,k. Gdy opad netto Pi,j,1 = 0 i ewaporacja Ei,j,1 = 0,

to speniony jest warunek selektywnoci i mona wtedy z zalenoci (13) i (16) obliczy

24

Rys. 5. Sposoby uwzgldniania w modelu matematycznym poboru wody przez korzenie rolin x – punkty brzegowe; – punkty, w których naley uwzgldni czon ródowy; –

po-zostae punkty; t k , j , i Q'

' – zmiana objtoci wody w elemencie Vi,j,k wywoana poborem

wody przez rolin

Fig. 5. Methods of representing water uptake by plant roots in mathematical model

x – boundary points; – points at which core element should be taken into account; –

other points; t k , j , i Q'

' – change in water volume in element Vi,j,k caused by water uptake by

plant

3.4.

PARAMETRYZACJA

PRZESTRZENI

OBJTEJ

MODELOWANIEM

Oprócz znajomoci warunków pocztkowych i brzegowych niezbdny jest opis przestrzeni objtej modelowaniem (Czamara 1998, Kowalik 1973, Olszta 1981,

Pywa-czyk iwsp. 1992, Russo 1998, Walczak i wsp. 1998, Wosiewicz 1986). W przypadku

równania Richardsa wspóczynnikami opisujcymi gleb s parametry w zalenociach

podanych przez van Genuchtena (zalenoci (3) i (4)). Do parametrów, które naley

wyznaczy, nale n, r, s,  i Ks. Ich wartoci charakteryzuj si przestrzennym

zró-nicowaniem (Olszta, Zaradny 1990, Si i wsp. 1999, Simùnek i wsp. 2000, Sawiski 2003, Vogej i wsp. 2000, Wosiewicz 1986, Zaradny 1987, Zhu, Mohanty 2002).

Para-metry temog by wyznaczane z tabel na podstawie rozkadu granulometrycznego

(Ge-nuchten van i wsp. 1991, Porbska i wsp. 2006, Sawiski 2003). Stosowane s w tym celu równie sieci neuronowe (Lamorski, Walczak 2002), a take tzw. zagadnie-nia odwrotne (Szulczewski 1990, 2003). W niniejszej pracy parametry te wyznaczono

poprzezkalibracj jednowymiarowego modelu przepywu wody opartego na równaniu

Richardsa przeprowadzon w kolumnach glebowych (Hills i wsp. 1989a, 1989b, Janik górna warstwa brzegowa,

upper boundary layer

dolna warstwa brzegowa, lower boundary layer

k

z j

x

pd roliny, plant stem

y y x j k z i Qt i,j,k Qt i,j,k warunek brzegowy, boundary condition czon ródowy, sink term

dolna warstwa brzegowa, lower boundary layer

25

i wsp. 2006). Polega to na takim doborze parametrów, aby wartoci wilgotnoci zmie-rzone aparatem TDR w warstwach kolumny glebowej i obliczone na podstawie symula-cji róniy si jak najmniej (byy najlepiej dopasowane). Ruch wody podczas ekspery-mentu wywoano podsikiem. Podczas symulacji komputerowych wykorzystano róni-cow posta równania Richardsa (zaleno (2)), sprowadzon do jednego wymiaru. Jako miar dopasowania wartoci zmierzonych i obliczonych przyjto sum odchyle

przecitnych (Bc) wyznaczon ze wzoru:

¦

  1 2 2 1 N k k c B N B , (17) gdzie:

Bc – suma odchyle przecitnych w warstwach od k = 2 do k = N – 1, m3m-3,

N – liczba warstw, na które podzielono kad kolumn, –,

k

B – odchylenie przecitne, w k-tej przestrzeni ustalonej kolumny, m3m-3.

Poniewa w trakcie prowadzonych symulacji zastosowano warunek brzegowy I

ro-dzaju, to wielko Bk wyznaczono wycznie w elementach, które le wewntrz

mode-lowanego obszaru. Odchylenie przecitne (Bk) dla k-tego elementu obliczano ze wzoru:

¦

W  W W W W T T N pom k obl k k N B 1 1 _{, (18)} gdzie:

Bk – odchylenie przecitne dla k-tej warstwy (k-tego elementu) kolumny, m3m-3,

NW – ilo porównywanych par wilgotnoci w k-tej warstwie kolumny podczas trwania

eksperymentu, –,

 – indeks kolejnych chwil czasowych, –,

obl k

W

T – wilgotno obliczona w k-tej warstwie kolumny w chwili , m3m-3,

pom k

W

T – wilgotno pomierzona w k-tej warstwie kolumny w chwili , m3m-3.

Wzór ten moe by uyty tyko wtedy, gdy czas mierzony jest staym krokiem.

Mi-nimaln warto bdu Bc otrzymano w sposób opisany w pracy dotyczcej

zastosowa-nia równazastosowa-nia Richardsa do opisu ruchu wody w betonie autoklawizowanym (Janik

i wsp. 2006). W pierwszej kolejnoci obliczono warto Bc dla parametrów n, r, s, 

i Ks, wybranych na podstawie tablic podanych przez van Genuchtena (1991), po wcze-

niejszym oznaczeniu rozkadu granulometrycznego czci mineralnych gleby.

Nastp-nie obliczono warto Bc, zmniejszajc i zwikszajc wartoci pierwszego parametru n.

Obliczenia prowadzono do chwili, gdy wielko Bc osigna warto minimaln. Na

tym etapie wielko n uznano jako ustalon. Nastpnie rozwizano rónicow posta równania Richardsa (zaleno (2)), zwikszajc i zmniejszajc wartoci nastpnego

parametru r i obliczano warto bdu Bc. Obliczenia ponownie prowadzono do chwili,

26

ustalony. Analogiczne obliczenia kontynuowano, zwikszajc i zmniejszajc wartoci nastpnych parametrów. Po przeprowadzeniu symulacji dla ostatniego zmienianego

parametru Ks zakoczono pierwsz ptl oblicze. Obliczenia w nastpnych ptlach

wykonywano w taki sam sposób, rozpoczynajc minimalizacj bdu Bc od zmian

usta-lonego w ptli pierwszego parametru n. Obliczenia prowadzono do chwili, w której

zmiana któregokolwiek z optymalizowanych parametrów powodowaa wycznie

zwikszenie wartoci Bc. Otrzymane w ten sposób wartoci parametrów gwarantuj

uzyskanie minimalnej wartoci odchylenia przecitnego Bk obliczanego z zalenoci

(17). Nie ma natomiast pewnoci, e bd one dokadnie odzwierciedla waciwoci opisywanej gleby.

W niniejszej pracy wspóczynniki n, r, s, , Ks, wyznaczono, a nastpnie

wprowa-dzono do równania Richardsa w czterech wariantach – podobnie jak warunki pocztkowe. W pierwszym wariancie, w którym równanie Richardsa sprowadzono do jednego wy-miaru, przestrze reprezentowana jest przez cile okrelone wartoci parametrów re-prezentujcych ca modelowan przestrze. Wartoci parametrów ustalano na podsta-wie tablic podanych przez van Genuchtena po wczeniejszym ustaleniu typu gleby w punkcie uznanym jako reprezentatywny. W drugim, dziki przeprowadzeniu procedur

identyfikacyjnych w wielu kolumnach Vi,j reprezentujcych przestrzenie powstae

w wyniku dyskretyzacji, róne wartoci parametrów wprowadzono bezporednio do trójwymiarowego równania Richardsa. W trzecim wariancie w pierwszej kolejnoci

przeprowadzono analiz geostatystyczn parametrów n, r, s, , Ks w sposób opisany

w rozdziale 3.1. Zbudowano, a nastpnie metod krigingu uporzdkowano mapy. Osta-tecznie, do równania Richardsa wprowadzono wartoci parametrów odczytane z map. W czwartym wariancie równanie znowu – jak w przypadku warunku pocztkowego – roz-patrzono jako jednowymiarowe, za wspóczynniki potraktowano jako zmienne losowe.

27

4

RÓWNANIE RICHARDSA Z RÓNYMI

WARIANTAMI DANYCH WEJCIOWYCH

W niniejszym rozdziale opisano róne sposoby prognozowania stosunków po-wietrzno-wodnych w wierzchniej warstwie gleby. Syntetycznie scharakteryzowano je w tabeli 1. W I wariancie, posiadajcym najwicej zaoe upraszczajcych, przyjto, e ca przestrze przedstawion na rysunku 1 reprezentuje wycznie jeden przekrój. Dlatego, model jest jednowymiarowy. Warunek pocztkowy stanowi wartoci

cako-witej wysokoci cinienia ) w chwili pocztkowej (indeksy i, j pomijamy). Górny _kpo

warunek brzegowy sformuowano jako warunek I rodzaju (zagadnienie Dirichleta). Wystpuje on wtedy, gdy w górnym brzegu dane s wartoci cakowitej wysokoci

ci-nienia ) (zaleno (2)) dla kadego kroku czasowego. Obliczono je analogicznie W1

jak wielkoci )_ipo_,j,k ze wzoru (5); z t rónic, e zamiast

h

_ipo_,j,1 wstawiono do wzoru

W 1

h (h₁W – wysoko cinienia w elemencie Vk dla kadego kroku czasowego 't).

Wiel-ko h1W wyznaczono po przeksztaceniu zalenoci (4) i znajc wilgotno

objtocio-w T₁W. Analogicznie ustalono dolny warunek brzegowy. W wariancie I przyjto

ponad-to, e gleba jest jednorodna w caej przestrzeni. Dlatego, parametry glebowe wystpuj-ce w równaniach van Genuchtena przyjmuj cile okrelone wartoci, wyznaczone na podstawie danych literaturowych. Wynikiem symulacji s cile okrelone wartoci wilgotnoci w przekroju uznanym jako reprezentatywny.

W II wariancie przyjto, e model jest trójwymiarowy. Warunek pocztkowy

sta-nowi wartoci )_ipo_,j,k dla kadego elementu Vi,j,k , (II wariant warunku pocztkowego –

rozdzia 3.1). Warunek brzegowy opisano w poprzednim akapicie, z tym e indeksy i, j zachowano. Konieczna jest wic znajomo dynamiki wilgotnoci podczas trwania

sy-mulacji w elementach Vi,j,k w górnej i dolnej warstwie. Parametry gleby modelowanej

przestrzeni w rzeczywistoci róni si w kierunkach poziomym i pionowym. Jednak zaproponowany w pracy sposób identyfikacji (rozdzia 3.4) pozwala na zrónicowanie ich wartoci jedynie w kierunku poziomym. W punktach, w których ze wzgldów

tech-28

nicznych nie udao si przeprowadzi identyfikacji, przyjto redni warto ze wszyst-kich wyznaczonych wartoci (II wariant parametryzacji przestrzeni zakada, e korela-cja waciwoci fizycznych gleby w ssiednich obszarach nie jest znana).

Wynik symulacji komputerowej stanowi wartoci wilgotnoci w kadym elemencie

Vi,j,k objtej modelowaniem. Wariant III równania Richardsa róni si tym od II, e do

formuowania warunku pocztkowego oraz opisu parametrów gleby wykorzystuje si metody geostatystyczne w celu budowy odpowiednich map. Opis wariantu IV zawarto w tabeli 1.

W czterach rozpatrywanych powyej wariantach do rozwizania równania zastoso-wano metod rónic skoczonych. Rónicow posta ze schematem jawnym przedsta-wiono w rozdziale 3 (zaleno (2)). Metoda rozwiza dla wariantów I, II i III, gdy warunki pocztkowe i parametry gleby przedstawione s w sposób zdeterminowany, nie wymaga oddzielnego komentarza. Szerszy opis potrzebny jest do zaprezentowania spo-sobu rozwizania równania Richardsa, gdy warunek pocztkowy realizowany jest jako zmienna losowa i uwzgldniono losowy charakter parametrów. Rozwizania analitycz-ne takiej klasy równa stochastycznych znaanalitycz-ne s jedynie dla najprostszych przypadków. Mona je znale w pracy prof. Sobczyka w rozdziale powiconym rozwizaniom analitycznym stochastycznych równa róniczkowych (Sobczyk 1991). Najczciej jednak stosuje si metody przyblione bazujce na rozwizaniach numerycznych (Sob-czyk 1987). W niniejszej pracy do rozwizania równania Richardsa z warunkiem po-cztkowym i parametrami gleby realizowanymi jako zmienna losowa zaproponowano sposób wykorzystujcy procedury losujce napisane w programie Matlab (Janik 2007b). Najpierw przeliczmy wszystkie moliwe przypadki sformuowania warunku pocztko-wego. Jeeli ponownie rozpatrzymy przestrze przedstawion na rysunku 1 (rozdz. 3),

w której wilgotno jest zmierzona w L˜M punktach w kadej warstwie, to warunek

pocztkowy moe by podany na (L˜M)N sposobów (L, M, N – rys. 1). Przykadowo, dla

L = 10, M = 10 i N = 5 powstaje (L˜M)N = 1010 zestawów warunków pocztkowych.* Chcc przeanalizowa wszystkie moliwe przypadki, naleaoby tyle razy rozwiza równanie Richardsa. Cho schematy jawne zastosowane w tej pracy cechuj si krótkim czasem oblicze (dla jednego zestawu ok. 0,02 s), to obliczenia musiayby trwa i tak ponad 6 lat. Dodatkowo, ilo moliwych rozwiza ulega zwikszeniu, gdy parametry gleby podane s równie w postaci zmiennej losowej. W zwizku z powyszym, zapro-ponowano nastpujcy sposób postpowania, który zaprezentowano równie w pracy autora (Janik 2007b), gdzie rozpatrzono wycznie losowy warunek pocztkowy. W pierwszej kolejnoci wybrano w sposób losowy ma liczb zestawów warunków

pocztkowych i parametrów opisujcych gleb (np. nL = 5). Taki sposób losowania

spowodowa, e wartoci wilgotnoci w warunku pocztkowym oraz wartoci parame-trów opisujcych gleb nie pochodziy z jednego przekroju. Dla kadego wylosowanego zestawu rozwizano rónicow posta równania Richardsa.

Tabela 1 Table 1

Warianty modeli do prognozowania warunków powietrzno-wodnych

w czynnej warstwie gleby

Variants of models for prediction of air-water conditions in active layer of soil

Komentarz Comments

Warunek pocztkowy

Initial condition

Warunek brzegowy Boundary conditions Parametry w równaniu Parameters in equation

Wynik symulacji

Result of simulation

Wariant I Variant I

Przyjmuje si, e jeden przekrój reprezentuje przestrze, model roz- patrywany wycznie w kierunku pionowym It is assumed that one section represents 3D space, the model is considered only in the vertical direction Okrelany na podstawie wilgotnoci w przekroju reprezentatywnym Determined on the basis of moisture values in repre- sentative section

Górny i dolny I rodzaju, Upper and lower of the 1st kind

Wartoci parametrów n ,  r , s ,  i K s ustalane na

pod-stawie danych literaturo- wych bd za pomoc pro- cedury identyfikacyjnej Values of parameters

n ,  r , s ,  and K s are determined

on the basis of literature data or by means of identifi- cation procedure Warto wilgotnoci w przekroju reprezenta- tywnym Value of moisture in representative section

Wariant II Variant II Model przestrzenny Spatial (3D) model Okrelany na podstawie wartoci wilgotnoci zmie- rzonych w kadym

punkcie

przestrzeni powstaej na skutek dyskretyzacji obszaru Determined on the basis of moisture value at every point of space formed through discretisation of the area Górny i dolny I rodzaju Upper and lower of the 1st kind

Obliczone bezporednio na podstawie procedury identy- fikacyjnej Calculated directly on the basis of identification pro- cedure Warto wilgotnoci w kadym punkcie przestrzeni powstaej w wyniku dyskretyzacji Value of moisture at every point of space formed through discre- tisation

Tabela 1 cd. _{Table 1 cont.}

Wariant III Variant III Model przestrzenny Spatial (3D) model

Okrelany na podstawie wartoci wilgotnoci w kadym punkcie odczyta- nych z map wilgotnoci Determined on the basis of moisture values at every point read from moisture distribution maps Górny i dolny I rodzaju Upper and lower of the 1st kind Wartoci parametrów od- czytane z map Values of parameters read from maps

Mapy wilgotnoci sporzdzone na pod- stawie danych otrzyma- nych z symulacji Maps of moisture dis- tribution created on the basis of data from simulations

Wariant IV Variant IV Model jednowymiarowy Mono-dimensional (1D) model Realizowany jako zmienna losowa wilgotnoci Realized as a random variable of moisture Górny i dolny I rodzaju Upper and lower of the 1st kind W postaci zmiennej losowej As a random variable

Wilgotno w przekroju reprezentujcym obszar podana jako zmienna losowa Moisture values in section representing the area, given as a random variable

31

Symulacj prowadzono dla okrelonych warunków brzegowych i okrelonego kroku

czasowego. Czas symulacji okrelono jako T1. W wyniku symulacji komputerowych,

w kadej warstwie, otrzymywano nL = 5 zestawów wartoci wilgotnoci k, które

przed-stawiono w postaci histogramów wilgotnoci. Nastpnie prowadzono analogiczne

sy-mulacje, z tym e za kadym razem zwikszano liczb nL. Symulacje kontynuowano do

momentu, w którym ksztat uzyskanych histogramów w kadej warstwie nie zmienia

si jednoczenie ze wzrostem liczby nL. Jako kryterium stabilnoci ksztatu

histogra-mów przyjto wielko EnL, któr dla ustalonego k wyznaczono ze wzoru (Janik 2007b):

 

T

 

T T T T d f f E_{nL Nmax nL} max min 

³

, (19) gdzie:

 

T max N

f – gsto rozkadu wilgotnoci dla Nmax, –,

 

T

nL

f – estymator funkcji

f

 

T

dla nL losowa, –,

Nmax – graniczna liczba losowa, –,

T – wilgotno objtociowa, m3

m-3,

Tmax – maksymalna wilgotno objtociowa zmierzona w k-tej warstwie, m3m-3,

Tmin – minimalna wilgotno objtociowa zmierzona w k-tej warstwie, m3m-3.

W zalenoci (19) pojawia si wielko Nmax. Jest to przyjta liczba losowa, która

gwarantuje, e wielko EnL bdzie staa jednoczenie ze wzrostem Nmax. Liczb Nmax

ustalono na podstawie eksperymentu numerycznego. Ostatecznie, liczb losowa nL

naley uzna za wystarczajc, gdy EnL nie bdzie zmieniao swojej wartoci. Na

ostat-nim etapie przeprowadzono ocen dobroci dopasowania histogramów danych uzyska-nych z pomiarów do histogramów uzyskauzyska-nych na podstawie dauzyska-nych zmierzouzyska-nych. Jako kryterium dobroci dopasowania danych empirycznych do danych uzyskanych z

symula-cji przyjto wielko E, któr zdefiniowano w nastpujcy sposób:

 

T

 

T T T T W f f d E _{m pom} max min 

³

, (20) gdzie: W

E – miara dopasowania histogramu uzyskanego na podstawie oblicze do histogra-

mów uzyskanych na podstawie symulacji, –,

 

T

m

f – estymator funkcji f

 

T uzyskany na podstawie symulacji, –,

 

T

pom

f – estymator funkcji f

 

T uzyskany na podstawie pomiarów, –,

pozostae oznaczenia jak w zalenoci (19).

Ostatecznie, wyniki uzyskane z symulacji przeprowadzonej zgodnie z czterema wa-riantami porównano z wynikami uzyskanymi na podstawie bada polowych wykona-nych na uytku kowym.

32

5

PRAKTYCZNE ZASTOSOWANIA

TECHNIKI TDR

5.1.

PRZYKADY

WYZNACZANIA

WARUNKÓW

POCZTKOWYCH

Poniej przedstawiono przykady zastosowania techniki TDR i metod geostaty-stycznych do formuowania rónych wariantów warunków pocztkowych w równaniu Richardsa. W tym celu przeprowadzono badania polowe w miejscowoci Brenna w wo-jewództwie lskim. Pooenie obiektu badawczego uwarunkowano spenieniem zaoe opisanych w rozdziale 3. Na rysunku 6 przedstawiono usytuowanie poletka dowiadczalnego stanowicego uytek kowy wraz z zaznaczonymi przekrojami,

w których zmierzono wilgotno objtocioww 100 punktach na gbokoci z1 = 2,5 cm

od powierzchni terenu, nastpnie w 100 punktach na gbokoci z2 = 7,5 cm, z3 = 12,5 cm,

z4 = 17,5 cm oraz z5 = 22,5 cm. (W dalszej czci pracy punkty na gbokoci z1 = 2,5 cm

nazywano „warstwa z1 = 2,5 cm” itd.). Odlego krawdzi poletka od rowu wynosi 150

m i od utwardzanej drogi okoo 75 m. Gwarantuje to, e ani rów, ani utwardzona droga nie miay wpywu na wilgotno w jego wntrzu. Do pomiaru zastosowano rczny apa-rat TDR z czujnikami polowymi FP/m. W kadym punkcie wilgotno mierzono

jedno-krotnie. Pomiar przeprowadzono w dniu 30 sierpnia w 2006 r. w godz. od 1420 do 1535.

Dodatkowo, w trakcie prowadzenia pomiarów w rodkowym przekroju obszaru (prze-krój nr 56 na rys. 6) prowadzono cigy monitoring wilgotnoci, co 1 minut na

gbo-kociach od z1 do z5 (przy braku zmian wilgotnoci odczyty wystarczyo prowadzi z

duszym krokiem czasowym, np. 't = 15 min). Stwierdzono, e czynnik czasu

w trakcie ustalania warunku pocztkowego, tj. w cigu 1 godz. i 15 min nie wpywa na dokadno pomiaru, poniewa rónica pomidzy maksymaln i minimaln

wilgotno-ci wyniosa jedynie 0,003 m3m-3.

W wariancie I zaoono, e jeden przekrój reprezentuje ca modelowan przestrze przedstawion na rysunku 6. Przyjto, e jest to przekrój nr 56 o wspórzdnych

(10 m; 10 m). Zmierzone wilgotnoci wyniosy: Tz561 0,297m

3

m-3, Tz562 0,288m

3

33

246 0 56 3 , z T m3m-3, T_z56₄ 0,235m3m-3, T_z56₅ 0,215m3m-3. Zaoenie reprezentatywnoci

tak sformuowanego warunku pocztkowego musi generowa bdy podczas modelowania.

Rys. 6. Usytuowania poletka dowiadczalnego w Brennej

a – poletko dowiadczalne; b, c – zabudowania gospodarcze; i, j – indeksy punktów po-miarowych w kierunkach x i y

Fig. 6. Location of experimental plot in Brenna

a – experimental plot; b, c – farm buildings; i, j – indexes of measurement points in direc-tions x and y

Wystarczy porówna wartoci wilgotnoci w przekroju nr 69 o wspórzdnych

(16 m; 12 m), gdzie: Tz691 0,274m 3 m-3, Tz692 0,278m 3 m-3, Tz693 0,266m 3 m-3, 221 0 69 4 , z T m3m-3, Tz69₅ 0,197m 3 m-3 bd w przekroju 54 o wspórzdnych (6 m; 10 m)

– tam wartoci wilgotnoci wyniosy Tz541 0,334m

3 m-3, T_z54₂ 0,262m3m-3, 285 0 54 3 , z T m3m-3, T54_z4 0,263m3m-3, Tz545 0,269m 3 m-3. Na rysunku 7 przedstawiono

inne moliwe przypadki sformuowania warunku pocztkowego. Jest to klasyczne podanie warunku pocztkowego, gdy równanie dotyczy przestrzeni jednowymiarowej (wariant I formuowania warunku pocztkowego).

szczegó a, detail a b c 75 m 15 0 m Brenna Skoczów a 2 4 6 8 18 16 14 10 12 0 y [m] i = 1,2,3, … ,L 8 6 j = 1, 2, 3, … , M i, j i+2, j i+1, j+3 1 2 3 21 22 23 65 66 67 68 69 91 92 93 97 98 99 100 2 4 10 12 14 16 18 x [m] 55