Modelowanie spływu powierzchniowego w obszarach miejskich z wykorzystaniem GIS

Leszek Ziółkowski

Modelowanie spływu powierzchniowego w obszarach miejskich z wykorzystaniem GIS

praca magisterska

wykonana pod kierunkiem: dr hab. Jacka Kozaka, prof. UJ

Praca magisterska złożona zgodnie z wymogami programu studiów magisterskich UNIGIS (Master of Science) w zakresie "Geographical

Information Science & Systems"

Uniwersytet Jagielloński w Krakowie, Uniwersytet Parisa Lodrona w Salzburgu

2013

2 Oświadczam, że wszystkie źródła wykorzystane w pracy zostały wymienione zgodnie z regułami cytowania. Niniejszą pracę wykonałem samodzielnie. Praca ta nie była i nie będzie składana jako praca dyplomowa w innym miejscu.

Data ………... Podpis………

3 Składam serdeczne podziękowania mojemu promotorowi Panu dr hab. Jackowi Kozakowi, prof. UJ, za pomoc i cenne uwagi udzielone w czasie pisania niniejszej pracy.

4 Spis treści

1 Wstęp ... 5

1.1 Cel pracy ... 6

2 Charakterystyka obszaru badań ... 7

2.1 Ukształtowanie terenu ... 8

2.2 Klimat ... 9

2.3 Użytkowanie ziemi ... 10

2.4 Stosunki wodne i sieć kanalizacji deszczowej ... 11

3 Dane wejściowe ... 13

3.1 Numeryczny model terenu ... 13

3.2 Dane o użytkowaniu ziemi ... 13

3.3 Kanalizacja deszczowa ... 14

4 Metody ... 15

4.1 Uproszczony dwuwymiarowy model spływu powierzchniowego ... 16

4.1.1 Założenia teoretyczne ... 16

4.1.2 Dane wejściowe do modelu ... 17

4.1.3 Warunki brzegowe modelu ... 19

4.1.4 Struktura modelu ... 21

4.2 Symulacja spływu powierzchniowego w programie ANUGA ... 21

4.2.1 Założenia teoretyczne ... 21

4.2.2 Części składowe modelu ... 22

4.3 Porównanie modeli ... 23

5 Wyniki ... 25

5.1 Wyniki symulacji spływu powierzchniowego w modelu uproszczonym w programie PCRaster ... 28

5.2 Wyniki symulacji spływu powierzchniowego w modelu ANUGA ... 32

5.3 Porównanie wyników symulacji z obydwu modeli ... 34

6 Dyskusja ... 38

7 Wnioski ... 41

8 Literatura ... 42

Spis rycin ... 45

Spis tabel ... 45

5 1 Wstęp

Powodzie miejskie stanowią jedno z większych i najbardziej kosztownych zagrożeń naturalnych. Ze względu na dużą gęstość zaludnienia oraz koncentrację zabudowy i infrastruktury technicznej, powodzie miejskie powodują większe szkody niż powodzie na obszarach wiejskich. Duże szkody i zniszczenia są często efektem nagłych powodzi lokalnych typu „flash flood”, które charakteryzuje duża objętość wody i krótki czas trwania. Przykładem tego typu powodzi jest zdarzenie, które miało miejsce w Gdańsku w lipcu 2001 r. W wyniku nawalnego opadu, trwającego 4 godziny, podtopiona została duża część miasta. Szkody oszacowano na około 150 mln złotych, a liczba poszkodowanych osób wyniosła 5000 (Majewski 2008). Procesy urbanizacji sprzyjają zwiększaniu się powierzchni nieprzepuszczalnych w zlewni. Intensywne opady deszczu na terenach z dużym udziałem powierzchni nieprzepuszczalnych mogą doprowadzić do gwałtownych spływów powierzchniowych, powodując lokalne podtopienia terenów najniżej położonych. Szacuje się, że urbanizacja (zwiększenie powierzchni obszarów uszczelnionych) przyczynia się nawet do sześciokrotnego zwiększenia odpływu (Zevenbergen i in. 2011, Hapuarachchi i in. 2011).

Zazwyczaj systemy kanalizacji deszczowej projektowane są na opady występujące raz na 10 lat (prawdopodobieństwo 10%), ponieważ budowa i utrzymanie kanalizacji deszczowej mogącej przeprowadzić wody o mniejszym prawdopodobieństwie wystąpienia byłoby zbyt kosztowne (Schmitt i in. 2004). W związku z tym, w wyniku bardzo intensywnych opadów, często dochodzi do sytuacji, kiedy przekroczona zostaje przepustowość kanalizacji i nadmiar wody wydostaje się na powierzchnię terenu tworząc spływ powierzchniowy. Analizowanie i badanie pracy kanalizacji deszczowej w odpowiedzi na zdefiniowany opad jest obecnie powszechnie przeprowadzane z wykorzystaniem modeli hydrodynamicznych. Najbardziej powszechną grupę stanowią modele jedno­wymiarowe (1D), których prekursorem jest model EPA SWMM stworzony przez Agencję Ochrony Środowiska Stanów Zjednoczonych (U.S.

Environmental Protection Agency). Jednakże za pomocą tych modeli nie jest możliwa symulacja spływu powierzchniowego – nadmiar wody magazynowany jest w wirtualnym zbiorniku nad przeciążoną studzienką i powraca do kanalizacji deszczowej w momencie, gdy znowu będzie to możliwe (Mark i in. 2004). W ostatnich latach silny nacisk kładziono na połączenie i integrację modeli kanalizacji deszczowej oraz modeli spływu powierzchniowego (Mark i in. 2004). W wyniku tych badań

6 powstała grupa modeli 1D/1D (np. MIKE Urban, Infoworks CS) – jedno­wymiarowy model kanalizacji oraz spływu powierzchniowego oraz 1D/2D – jedno­wymiarowy model kanalizacji z dwu­wymiarowym modelem spływu powierzchniowego (np. MIKE MOUSE, Sobek Urban, SIPSON/UIM). Porównanie poszczególnych modeli nie było do tej pory szeroko podejmowane w literaturze (Lhomme i in. 2006, Leandro i in.

2009). Dodatkowo większość modeli zintegrowanych dostępna jest jako rozwiązania komercyjne.

Każde z wyżej zaproponowanych rozwiązań ma swoje wady i zalety. Prostsze modele wymagają mniej pracy i są mniej kosztowne. Dodatkowo, nie zawsze najbardziej zaawansowane rozwiązanie jest najbardziej użyteczne, dlatego złożoność modelu powinna być zawsze dostosowana do celów danego projektu (Lan­Anh i in. 2012).

Wiele z programów do modelowania powodzi wykorzystuje narzędzia i bazuje na systemach informacji geograficznej (GIS;np. Bates, De Roo 2000, Lhomme i in.

2004, Lhomme i in. 2006, Boonya­Aroonnet i in. 2007).

1.1 Cel pracy

Celem niniejszej pracy jest przetestowanie i porównanie wybranych metod symulacji spływu powierzchniowego wody wydostającej się z przeciążonej kanalizacji deszczowej. Metody te zostały opracowane w środowisku GIS dla wybranej zlewni w Kielcach. Wybór obszaru, dla którego przeprowadzone zostały obliczenia podyktowany był dostępnością danych dotyczących przebiegu sieci kanalizacji deszczowej.

Porównane zostały dwie metody symulacji spływu powierzchniowego. W pierwszej z metod obliczenia dokonywane były pomiędzy pikselami rastra przy założeniu ruchu jednostajnie ustalonego. Do tego celu wykorzystane zostały funkcje i narzędzia w oprogramowaniu PCRaster. Druga metoda obliczeń opierała się na modelu hydrodynamicznym programu ANUGA, w którym rzeźba terenu odwzorowana jest w postaci siatki nieregularnych trójkątów.

Otrzymane wyniki zostały porównane w środowisku GIS. Ostatecznie dokonana została ocena przydatności zaproponowanych metod do symulacji spływu powierzchniowego w mieście.

7 2 Charakterystyka obszaru badań

Opracowanie dotyczy fragmentu miasta Kielce odwadnianego przez system kanalizacji deszczowej do rzeki Silnicy. Teren ten obejmuje zlewnie kolektora wód opadowych Si­8 o powierzchni 1,64 km², w skład której wchodzą centralne rejony Kielc wraz z osiedlami Wielkopole i Czarnockiego (ryc. 2.1).

Zgodnie z regionalizacją fizycznogeograficzną według Kondrackiego (2000) obszar objęty analizą w całości leży w mezoregionie Góry Świętokrzyskie (makroregion Wyżyna Kielecka, podprowincja Wyżyna Małopolska, prowincja Wyżyny Polskie).

Ryc. 2.1 Obszar badań na tle mapy topograficznej

8 2.1 Ukształtowanie terenu

Miasto Kielce położone jest w Niecce Łagowskiej, zwanej także Padołem Kielecko­

Łagowskim. Obniżenie to osiąga szerokość 9 km i charakteryzuje się występowaniem lekko falistych lub płaskich wzniesień, z których wystają ostańcowe pagóry z wapieni skalistych ­ Kadzielnia, Karczówka, góra Cmentarna, Słoneczna i Wietrznia. Dno Niecki Łagowskiej rozcięte zostało przez doliny poprzeczne (m.in.. Bobrzy i Silnicy) i podłużne wypełnione miąższymi osadami żwirów, piasków, iłów i mułków plejstoceńskich (Klimaszewski 1972).

Rzeźba obszaru badań została silnie przekształcona przez działalność człowieka.

Charakterystyczne jest występowanie antropogenicznych form rzeźby takich jak:

podcięcia i nasypy drogowe, starasowane zbocza i stoki. W wielu miejscach teren został zniwelowany oraz nastąpiły zmiany nachyleń. Zmiany te wraz z silną urbanizacją miały wpływ na zmianę naturalnych kierunków spływu wody.

Południową granice zlewni stanowi ostańcowy pagór wapienny wchodzący w skład większego zespołu o nazwie Wietrznia. Zlokalizowany jest tam również najwyższy punkt obszaru badań o wysokości 315,8 m n.p.m. Najniższy punkt (259,5 m n.p.m.) położony jest w okolicy ujścia kolektora kanalizacji deszczowej do Silnicy. Stoki wapiennych pagórów w południowej części zlewni charakteryzują się największymi nachyleniami terenu. Lokalnie duże spadki występują na skarpach nasypów drogowych (ryc. 2.2).

9 Ryc. 2.2 Ukształtowanie terenu badanej zlewni

2.2 Klimat

Kielce leżą w obrębie regionu klimatycznego wyżyn śląsko­małopolskich, którego część stanowią Góry Świętokrzyskie. Średnia roczna suma opadu z wielolecia 1954­2003 dla stacji Kielce­Suków (położona około 8 km od centrum miasta) wynosi 634,2 mm. Najwyższa roczna suma opadów w tym czasie zanotowana została w 2001 r.

i wyniosła 948,1 mm. Najwyższe miesięczne sumy opadów występują w lipcu, natomiast najniższe w lutym. Średnia roczna temperatura dla tego wielolecia wyniosła 7,3°C.

Zurbanizowana część miasta charakteryzuje się występowaniem podwyższonej temperatury w stosunku do terenów sąsiednich, co związane jest z występowaniem tzw.

miejskiej wyspy ciepła. Efektem tego jest pojawianie się wyraźnych maksimów w dobowym rozkładzie opadów w godzinach popołudniowych i wieczornych.

10 Wzmożone procesy konwekcyjne odpowiadają również za wzrost częstości i intensywności opadów nawalnych w centralnej części miasta (Ciupa 2009).

2.3 Użytkowanie ziemi

Obszar badań, ze względu na centralne położenie w mieście jest silnie zurbanizowany.

W użytkowaniu dominuje zabudowa zwarta (44%) oraz tereny przemysłowe, handlowe lub usług publicznych wraz z przyległymi parkingami (16%), które przeplatają się z zabudową luźną o dużym zgęszczeniu (8%). Południową część obszaru charakteryzującą się największymi nachyleniami terenu, zajmują ogródki działkowe i obiekty infrastruktury sportowej (13%) oraz tereny rolne i obszary semi­naturalne (5%). Dużą część terenu zlewni stanowią drogi i tereny z nimi związane (12%).

Z południa na północ zlewnię przecina dwujezdniowa ul. Źródłowa, natomiast wzdłuż północnej granicy zlewni przebiegają dwujezdniowe ul. IX Wieków Kielc i ul. Sandomierska. Na obszarze badanej zlewni brak jest obszarów leśnych, a miejskie tereny zielone zajmują 1% powierzchni zlewni (ryc. 2.3).

11 Ryc. 2.3 Użytkowanie ziemi w zlewni Si­8, na podstawie European Urban Atlas (2013) 2.4 Stosunki wodne i sieć kanalizacji deszczowej

Obszar objęty analizą obejmuje fragment zlewni rzeki Silnicy, będącej lewobrzeżnym dopływem Bobrzy. Rzeka ta z kolei uchodzi do Czarnej Nidy. Badana zlewnia jest typową zlewnią miejską o dużym stopniu uszczelnienia, co powoduje zwiększony i przyspieszony odpływ wód deszczowych do rzeki Silnicy. W związku z dużym stopniem zagospodarowania terenu, brak jest tutaj naturalnych cieków. Sieć hydrograficzna jest reprezentowana przez nieliczne rowy przechodzące często w kanały zamknięte sieci odprowadzania wód deszczowych.

Zlewnia kolektora wód opadowych Si­8 odpowiada za 4,8% z całości wód opadowych prowadzonych korytem Silnicy dla deszczu o prawdopodobieństwie wystąpienia 1%

(Gmina Kielce 2012). Długość sieci kanalizacji deszczowej w tej zlewni wynosi około 14 km, w związku z czym zagęszczenie sieci wynosi 8,5 km/km². Liczba studzienek

12 kanalizacji deszczowej w badanej zlewni wynosi 428 (ryc. 2.4).

Ryc. 2.4 Sieć kanalizacji deszczowej w zlewni SI­8 (Gmina Kielce 2010)

13 3 Dane wejściowe

Modelowanie spływu powierzchniowego wymaga dokładnych danych wysokościowych tj. numerycznego modelu terenu pochodzącego ze skaningu laserowego. Dodatkowo bardzo istotną informacją jest typ pokrycia terenu, od którego zależy współczynnik szorstkości uwzględniany w obliczeniach symulacyjnych. W celu zlokalizowania miejsc wypływu wody z przeciążonej kanalizacji deszczowej niezbędne było pozyskanie informacji o lokalizacji studzienek oraz przebiegu sieci kanalizacji deszczowej.

3.1 Numeryczny model terenu

Z Centralnego Ośrodka Dokumentacji Geodezyjnej i Kartograficznej pozyskano dwa arkusze numerycznego modelu terenu (NMT) o rozdzielczości 1 m obejmujące swym zasięgiem obszar badań. NMT powstał w wyniku interpolacji chmury punktów pochodzącej z lotniczego skaningu laserowego. Powyższe dane wysokościowe zostały utworzone w ramach projektu pn. „Informatyczny System Osłony Kraju przed nadzwyczajnymi zagrożeniami” (ISOK 2013). Pomiary dla wybranych obszarów miejskich (w tym Kielc) wykonane zostały z gęstością równą 12 punktów na m². Średni błąd wysokości NMT pochodzących z projektu ISOK wynosi 0,15 m dla powierzchni odkrytych i utwardzonych oraz 0,25­0,30 m dla powierzchni zalesionych (GUGIK 2013).

3.2 Dane o użytkowaniu ziemi

Dane dotyczące pokrycia terenu pozyskano z bazy danych European Urban Atlas udostępnianej przez Europejską Agencję Środowiska (European… 2013). Dane te dotyczą roku referencyjnego 2006 i obejmują miasta zamieszkałe przez więcej niż 100000 mieszkańców. Do klasyfikacji wykorzystano między innymi zobrazowania satelitarne o rozdzielczości przestrzennej nie mniejszej niż 2,5 m. Minimalna jednostka kartowania dla terenów antropogenicznych wynosiła 0,25 ha, a dla pozostałych obszarów ­ 1 ha. Dodatkowo dla terenów antropogenicznych dane dostępne są do czwartego poziomu szczegółowości według hierarchicznego podziału klas użytkowania ziemi zaproponowanego przez Europejską Agencję Środowiska. Dla pozostałych klas szczegółowość danych kończy się na pierwszym poziomie (tab. 3.1).

14 Tab. 3.1 Hierarchiczny podział na klasy użytkowania ziemi w bazie danych European

Urban Atlas (European… 2011)

Kod klasy Urban

Atlas Nazwa klasy

1 Tereny antropogeniczne

1.1 Zabudowa miejska

1.1.1 Zabudowa zwarta (poziom uszczelnienia >80%) 1.1.2 Zabudowa luźna (poziom uszczelnienia 10%­80%)

1.1.2.1 Zabudowa luźna o dużym zagęszczeniu (stopień uszczelnienia 50%­80%) 1.1.2.2 Zabudowa luźna o średnim zagęszczeniu (stopień uszczelnienia 30%­50%) 1.1.2.3 Zabudowa luźna o małym zagęszczeniu (stopień uszczelnienia 10%­30%) 1.1.2.4 Zabudowa luźna o bardzo małym zagęszczeniu (stopień uszczelnienia <10%) 1.1.3 Pojedyncze obiekty zabudowy

1.2 Tereny przemysłowe, handlowe, usług publicznych, wojskowe i komunikacyjne 1.2.1 Tereny przemysłowe, handlowe, usług publicznych i wojskowe

1.2.2 Tereny komunikacyjne i związane z komunikacją drogową i kolejową 1.2.2.1 Szybkie drogi tranzytowe i tereny z nimi związane

1.2.2.2 Pozostałe drogi i tereny z nimi związane 1.2.2.3 Koleje i tereny kolejowe

1.2.3 Porty

1.2.4 Lotniska

1.3 Kopalnie, wyrobiska i budowy

1.3.1 Miejsca eksploatacji odkrywkowej, wysypiska i zwałowiska

1.3.3 Budowy

1.3.4 Tereny bez określonego użytkowania

1.4 Zantropogenizowane tereny zielone nie użytkowane rolniczo 1.4.1 Miejskie tereny zielone

1.4.2 Tereny sportowe i wypoczynkowe

2 Tereny rolne, obszary roślinności semi­naturalnej, obszary podmokłe

3 Lasy

5 Wody

3.3 Kanalizacja deszczowa

Dane dotyczące przebiegu sieci kanalizacji deszczowej pozyskano z opracowań udostępnianych przez Urząd Miasta Kielce: Koncepcji zagospodarowania wód deszczowych dla Miasta Kielce (2010) oraz Wykonania koncepcji ochrony przeciwpowodziowej w zlewni Bobrzy, Silnicy, Sufragańca i Lubrzanki na obszarze miasta Kielce (2012).

15 4 Metody

W opracowaniu wykorzystano dwie metody symulacji spływu powierzchniowego bazujące na numerycznym modelu terenu pochodzącym z lotniczego skaningu laserowego. Pierwsza metoda polega na uproszczonym dwuwymiarowym modelowaniu spływu powierzchniowego, w którym przepływ wody pomiędzy pikselami wyrażony jest równaniem Manninga (ruch jednostajny ustalony). Metoda ta stosowana była między innymi do modelowania przepływu na terasach zalewowych (Bates, De Roo 2000), a także symulacji powodzi miejskich (Hunter i in. 2008). Druga z metod bazuje na dwuwymiarowym modelu hydrodynamicznym dostępnym w darmowym oprogramowaniu – ANUGA. Program ten do obliczenia przepływu, wykorzystuje równanie płytkiej wody, które rozwiązywane jest za pomocą metody objętości skończonych (Geoscience … 2010).

Numeryczny model terenu, wykorzystany do obliczeń, został zmodyfikowany poprzez dodanie informacji o położeniu i wysokości budynków. Wszelkie zabudowania stanowią barierę dla przepływu wody oraz mogą powodować zmianę kierunku jej płynięcia. Dane dotyczące położenia budynków pozyskano z bazy danych OpenStreetMap. Są to wolne dane udostępniane na podstawie licencji Open Database License. Wysokość wszystkich budynków przyjęto jako 10 m, gdyż dokładna wartość nie jest istotna, a jedynie fakt stanowienia przez nie bariery.

Do przeprowadzenia obliczeń hydraulicznych niezbędne jest określenie wartości współczynnika szorstkości dla całego obszaru objętego modelowaniem. Do tego celu wykorzystano informację o użytkowaniu ziemi zawartą w bazie danych European Urban Atlas. Poszczególnym klasom użytkowania nadano odpowiednie wartości współczynnika szorstkości Manninga n (tab. 4.1).

Tab. 4.1 Wartości współczynnika szorstkości Manninga dla poszczególnych klas użytkowania ziemi (Chow 1959)

Nazwa klasy

Wartość współczynnika Manninga n [s/m^1/3]

Zabudowa zwarta (poziom uszczelnienia >80%) 0,15

Zabudowa luźna o dużym zagęszczeniu (stopień uszczelnienia 50%­80%) 0,15 Tereny przemysłowe, handlowe, usług publicznych i wojskowe 0,067

Pozostałe drogi i tereny z nimi związane 0,025

Miejsca eksploatacji odkrywkowej, wysypiska i zwałowiska 0,05

16

Nazwa klasy

Wartość współczynnika Manninga n [s/m^1/3]

Tereny bez określonego użytkowania 0,04

Miejskie tereny zielone 0,025

Tereny sportowe i wypoczynkowe 0,56

Tereny rolne, obszary roślinności semi­naturalnej, obszary podmokłe 0,5

4.1 Uproszczony dwuwymiarowy model spływu powierzchniowego 4.1.1 Założenia teoretyczne

Prostym sposobem symulacji spływu powierzchniowego pomiędzy pikselami numerycznego modelu terenu jest zastosowanie równania ciągłości w odniesieniu do objętości wody jaką magazynuje każda komórka rastra. Zmiana objętości wody w pikselu w danej jednostce czasu jest równa ilości dopływającej i odpływającej z niego wody w jednostce czasu (Bates, De Roo 2000):

= + + + ,

gdzie:

V – objętość piksela [m³], t – czas [s],

Qd, Qg, Ql, Qp – dopływ (dodatni lub ujemny) z sąsiednich pikseli [m³/s].

Natomiast przepływ pomiędzy dwoma sąsiednimi pikselami i i j może zostać obliczony w uproszczony sposób za pomocą równania Manninga, które wyraża się wzorem (Puzyrewski, Sawicki 1998):

Q_ij=A_ij⋅R_ij^2/3⋅S_ij^1/2 n , gdzie:

Qij – przepływ między pikselem i i j [m³/s],

Aij – powierzchnia przekroju na granicy pikseli [m²], Rij – promień hydrauliczny [m],

Sij – spadek zwierciadła wody pomiędzy pikselami i i j,

17 n – współczynnik szorstkości Manninga [s/m^1/3].

Promień hydrauliczny jest to stosunek powierzchni przekroju czynnego do obwodu zwilżonego.

4.1.2 Dane wejściowe do modelu

Powyższy model zbudowano z wykorzystaniem oprogramowania PCRaster pozwalającego na przeprowadzanie symulacji zjawisk, których parametry zmieniają się zarówno w czasie jak i przestrzeni. Program ten daje możliwość tworzenia własnych modeli z wykorzystaniem wbudowanego języka skryptowego lub języka Python.

Obliczenia bazują na serii map rastrowych stanowiących dane wejściowe do modelu oraz kombinacji operatorów i funkcji. PCraster jest programem typu open source.

Uproszczony, dwuwymiarowy model spływu powierzchniowego bazował na dwóch zdefiniowanych w programie PCRaster funkcjach ­ dynamicwaveh oraz dynamicwaveq.

Funkcje te, na podstawie równań opisanych powyżej, obliczają stan wody oraz wartość przepływu w poszczególnych pikselach NMT w każdym kroku obliczeniowym. Każda komórka rastra traktowana była jako przybliżenie koryta cieku, aby możliwe było obliczenie przepływu wody między pikselami. W związku z powyższym, jako dane wejściowe do modelu musiały zostać utworzone następujące mapy rastrowe:

 rzędna dna koryta – wysokość terenu z NMT,

 kierunek spływu wody,

 współczynnik szorstkości,

 długość koryta,

 szerokość koryta przy dnie,

 głębokość koryta,

 kształt koryta.

Każda wejściowa mapa rastrowa musiała posiadać taki sam zasięg oraz rozdzielczość, aby każdemu pikselowi danej mapy odpowiadał dokładnie jeden piksel innej mapy.

Do obliczeń przyjęto rastry o rozdzielczości równej 1 m, dostosowanej do rozdzielczości NMT, aby uniknąć generalizacji niektórych form ukształtowania terenu mających wpływ na kierunek spływu wody.

18 Kierunki spływu wody wyznaczone zostały na bazie NMT z naniesionymi budynkami.

Dla każdego piksela określony został kierunek spływu wody do sąsiedniego piksela po linii największego spadku. Do tego celu wykorzystano funkcję lddcreate programu PCRaster, która nadaje komórkom rastra wartości od 0 do 8 w zależności od kierunku płynięcia wody (ryc. 4.1). Wartość 5 przyjmuje piksel wyznaczający lokalne zagłębienie bezodpływowe, do którego spływają wody z sąsiednich komórek. Duża część obszarów bezodpływowych spowodowana jest błędami w NMT powstałymi na etapie jego tworzenia. Tego typu lokalne zagłębienia powinny zostać usunięte. Program PCRaster daje możliwość decydowania, które zagłębienia zostaną usunięte na podstawie jego wielkości. Powierzchnia obszaru bezodpływowego obejmuje wszystkie piksele zlewni tego zagłębienia, których wysokość jest mniejsza od wysokości całkowitego wypełnienia tego zagłębienia (ryc. 4.2). Tak zdefiniowany obszar bezodpływowy może zostać scharakteryzowany przez 4 parametry:

 głębokość (outflow depth),

 pojemność (core volume),

 powierzchnię (core area),

 wysokość opadu powodującego zapełnienie zagłębienia.

Dla każdego piksela o wartości równej 5 program oblicza powyższe parametry. Jeżeli wartość któregokolwiek parametru będzie większa niż określona wartość progowa, to takie zagłębienie nie zostanie usunięte (PCRaster 2011). Do celów modelowania spływu powierzchniowego usunięte zostały wszystkie obszary bezodpływowe o głębokości mniejszej niż 0,15 m.

Ryc. 4.1 Wartości pikseli dla poszczególnych kierunków płynięcia wody

19 Ryc. 4.2 Schemat obszaru bezodpływowego oraz jego parametrów (źródło: PCRaster

2011)

Długość koryta, dla każdego piksela wyznaczono na podstawie mapy kierunków spływu. Dla czterech głównych kierunków długość koryta określono zgodnie z wielkością piksela, czyli jako 1 m. Dla pozostałych kierunków wartość tą ustalono na 1,44 m, czyli długość przekątnej kwadratu o boku 1 m. Program daje możliwość ustawienia długości koryta większej niż wielkość piksela, co pozwala na uwzględnienie m.in. meandrów przy analizach dokonywanych na rastrach o małej rozdzielczości przestrzennej. W obrębie jednego piksela może znajdować się wtedy meandrujący fragment rzeki, którego długość jest większa niż wielkość piksela.

Dla całego obszaru przyjęto założenie, iż piksele odzwierciedlają płytkie, prostokątne koryta między którymi odbywa się przepływ wody na całej ich szerokości. W związku z czym, szerokość koryta przy dnie ustalono na 1 m, a głębokość koryta – 0,1 m.

4.1.3 Warunki brzegowe modelu

Kolejnym etapem pracy było określenie warunków brzegowych modelu, tzn. wyznaczenie miejsc i sposobu, w jakim woda dostaje się na obszar uwzględniony w modelu i w jakim obszar ten opuszcza. Bazując na opracowaniu wykonanym przez Gminę Kielce (2010), w zlewni kolektora Si­8 wybrano studzienki podtopione i wylewające dla deszczu o prawdopodobieństwie wystąpienia 20% (raz na 5 lat).

Dla obu przypadków przyjęto stały wypływ wody o natężeniu odpowiednio 0,1 i 0,2 m³/s przez 60 minut trwania symulacji. Wartości te zostały ustalone arbitralnie, na podstawie analizy wyników modelowania przepływu wody w kanalizacji deszczowej z wspomnianego opracowania. Takie podejście wynika również z faktu, iż celem pracy jest porównanie wyników modelowania spływu powierzchniowego dla dwóch modeli,

20 a nie określenie zasięgu zalewu dla konkretnego zdarzenia. Miejsce, w którym woda mogła opuścić model, wyznaczono w korycie rzeki Silnicy w sąsiedztwie zamknięcia zlewni Si­8. Tak zdefiniowany warunek brzegowy pozwalał na swobodny wypływ wody z modelu o stałym natężeniu równym 0,5 m³/s.

Powodzie miejskie powodowane są przez krótkotrwałe, lecz bardzo intensywne opady nawalne. Z tego względu czas symulacji ograniczony został do 90 minut. Do 60 minut, w których następowało rozlewanie się wody ze studzienek kanalizacji deszczowej, dodano 30 minut trwania symulacji, w których natężenie wypływu ze studzienek wyniosło 0 m³/s. Pozwoliło to spływającej wodzie na dotarcie do lokalnych zagłębień terenu lub na opuszczenie modelu. W tym czasie nie zezwolono pozostałej na powierzchni wodzie na powrót do systemu kanalizacji deszczowej.

Bardzo istotnym parametrem obliczeń hydrodynamicznych jest krok czasowy, z jakim mają być one przeprowadzane. Czas jest jedną ze zmiennych, które mają decydujący wpływ na ilość wody przepływającej między komórkami obliczeniowymi, a co za tym idzie na stabilność modelu. Krok czasowy obliczeń powinien być na tyle mały, aby w trakcie jego trwania woda nie przemieściła się dalej niż o jeden piksel.

Na sprawdzenie tego warunku pozwala bezwymiarowy parametr zwany liczbą Couranta (Courant i in. 1967):

= ∆

∆ ≤ 1, gdzie:

U – prędkość maksymalna fali [m/s], Δt – krok czasowy obliczeń [s],

Δx – wielkość komórki obliczeniowej [m].

Z powyższego równania wynika, iż im mniejsza komórka obliczeniowa tym mniejszy powinien być krok czasowy obliczeń, aby spełnić warunki stabilności modelu hydrodynamicznego i osiągnąć wiarygodne wyniki. Niewielki rozmiar piksela, w przyjętym modelu, wymusił zastosowanie niewielkiego kroku czasowego obliczeń wynoszącego 0,1 s. Oba te czynnik mają decydujący wpływ na czas potrzebny do przeprowadzenia symulacji. Wyniki obliczeń – stan wody oraz wartość przepływu – były zapisywane co 100 krok czasowy, czyli co 10 s.

21 4.1.4 Struktura modelu

Program PCRaster daje możliwość wykorzystania wyników obliczeń z danego kroku czasowego jako dane wejściowe do obliczeń w następnym kroku. W przypadku stworzonego modelu spływu powierzchniowego, obliczony stan wody dla każdego piksela stanowił w danym kroku obliczeniowym dane wejściowe do obliczeń w następnym kroku. Dla każdego kroku czasowego aktualizowano również informację o kierunku spływu wody, obliczanego na podstawie NMT, do którego dodana została wartość stanu wody w danym pikselu. W ten sposób osiągnięto efekt, w którym spływająca po powierzchni woda może zmieniać kierunek płynięcia w zależności od aktualnego stanu wody w danej komórce oraz komórkach sąsiednich. Schemat zbudowanego modelu przedstawiony został na ryc. 4.3.

Ryc. 4.3 Schemat obliczeń w uproszczonym modelu spływu powierzchniowego 4.2 Symulacja spływu powierzchniowego w programie ANUGA

4.2.1 Założenia teoretyczne

Program ANUGA jest dwuwymiarowym modelem hydrodynamicznym, który bazuje na równaniu płytkiej wody rozwiązywanym metodą objętości skończonych. Obliczenia dokonywane są na siatce nieregularnych trójkątów (ryc. 4.4). Równanie to przyjmuje postać:

22 ℎ+ + = ,

gdzie:

h – głębokość wody,

q_x – przepływ w kierunku x, q_y – przepływ w kierunku y, q_m – dopływ/odpływ źródłowy.

Obliczenia głębokości wody oraz przepływu dokonywane są w każdej komórce obliczeniowej – centroidzie trójkąta. Szczegółowy opis powyższego równania, z podziałem na równania zachowania masy i pędu, można znaleźć m.in. w opracowaniu Szymkiewicz (2000) oraz w przewodniku do programu ANUGA (Geoscience … 2010).

Ryc. 4.4 Siatka trójkątów w metodzie objętości skończonych (źródło: Roberts i in.

2008)

4.2.2 Części składowe modelu

Dane wejściowe do niniejszego modelu stanowiły: NMT oraz mapa współczynnika szorstkości Manninga n. Wprowadzenie tych danych do programu ANUGA odbywało się poprzez import pliku tekstowego zawierającego następujące informacje o punktach:

 współrzędną x,

 współrzędną y,

 współrzędną z,

 wartość współczynnika szorstkości.

W pierwszej kolejności dokonano konwersji rastrowego NMT do wektorowej warstwy punktowej. Następnie dla każdego punktu odczytano wartość współczynnika szorstkości. Tak przygotowane dane zapisano do pliku tekstowego o odpowiedniej

23 strukturze.

Na podstawie danych z poszczególnych punktów w programie ANUGA stworzony został NMT w postaci siatki trójkątów z zapisaną informacją o wartości współczynnika szorstkości dla każdego punktu. Użytkownik ma możliwość zmiany parametrów generowanej siatki poprzez określenie maksymalnej wielkości trójkąta oraz minimalnego dopuszczalnego kąta. Pierwszy z tych parametrów ustawiono na 1 m, a dla drugiego pozostawiono wartość domyślną (28°).

W niniejszym modelu przyjęto takie same założenia co do natężenia i czasu trwania wypływu wody z systemu kanalizacji deszczowej jak w uproszczonym modelu spływu powierzchniowego. Dla studzienek podtapianych i wylewających przyjęto stały wypływ wody o natężeniu odpowiednio 0,1 i 0,2 m³/s przez 60 minut trwania symulacji.

Warunki brzegowe zdefiniowane zostały na odpowiednich krawędziach trójkątów, które odpowiadały lokalizacji studzienki.

Krok czasowy obliczeń w programie ANUGA obliczany jest automatycznie na podstawie wielkości trójkąta i maksymalnej prędkości przepływu wody między trójkątami. Użytkownik ustala tylko czas trwania symulacji oraz interwał z jakim zapisywane mają być wyniki. Parametry te ustawiono odpowiednio na 5400 i 10 s.

Program ANUGA nie posiada interfejsu graficznego umożliwiającego budowanie modelu hydraulicznego. Cały ten proces odbywa się poprzez skrypty języka Python wykorzystujące odpowiednie moduły i funkcje ANUGA. Program umożliwia przeglądanie wyników w postaci animacji.

4.3 Porównanie modeli

W celu porównania wyników symulacji spływu powierzchniowego z obu modeli pozyskano informację o:

 maksymalnym stanie wody w każdej komórce obliczeniowej,

 maksymalnej wartości przepływu.

Wynikowe rastry przekonwertowano do warstw wektorowych (wielobok) przedstawiających maksymalny zasięg zalewu. Następnie dokonano analizy powierzchni zalewu w poszczególnych obszarach zlewni Si­8.

Do porównania głębokości zalewu oraz wartości przepływu wykorzystano narzędzia

24 algebry map (odjęcie dwóch map rastrowych). Natomiast do wyznaczenia obszarów które zostały zalane w obu modelach, a które tylko w jednym z modeli, wykorzystano funkcję Combine programu ArcGIS. Narzędzie to nadaje pikselom unikalne wartości dla każdej unikalnej kombinacji wartości pikseli z rastrów wejściowych (ryc. 4.5).

Aby możliwe było zastosowanie funkcji Combine, rastry wynikowe przeklasyfikowano w ten sposób, aby każdy zalany piksel otrzymał wartość 1, a pozostałe ­ 0.

Ryc. 4.5 Schemat działania funkcji Combine (ESRI 2012)

25 5 Wyniki

W obrębie zlewni Si­8 znajdowało się 948 budynków, które zostały wprowadzone do NMT (ryc. 5.1). Drugą bardzo istotną informacją niezbędną do przeprowadzenia obliczeń hydraulicznych jest mapa współczynnika szorstkości. W zlewni Si­8 przeważają obszary o wysokiej wartości współczynnika szorstkości Manninga n – 0,15 s/m^1/3 – odpowiadające zwartej zabudowie oraz zabudowie luźnej o dużym zagęszczeniu. Ponad trzykrotnie mniejszą powierzchnię zajmują obszary przemysłowe i handlowe o szorstkości wynoszącej 0,067 s/m^1/3. Tereny dróg stanowiące najbardziej korzystne warunki dla przepływu wody (0,025 s/m^1/3) obejmują powierzchnię około 0,2 km² (ryc. 5.2).

Ryc. 5.1 NMT z wprowadzoną informacją o wysokości budynków

26 Ryc. 5.2 Mapa wartości współczynnika szorstkości Manninga n dla obszaru objętego

modelowaniem

Na podstawie opracowania Gminy Kielce (2010) z obszaru zlewni Si­8 wybrano 48 podtopionych i 28 wylewających studzienek (ryc. 5.3). Większość z nich skupia się w trzech obszarach. Pierwszy, to okolice ul. Bohaterów Warszawy między ul. Szymanowskiego i Szczecińską, gdzie zlokalizowanych jest 9 studzienek podtopionych i 4 wylewające. Wzdłuż ul. Źródłowej, w jej najbliższym sąsiedztwie, znajdują się 23 studzienki, z czego 18 to studzienki podtopione. Obszar z największą liczbą studzienek wylewających (9) położony jest wzdłuż ul. Kościuszki i ul. Bodzentyńskiej. Pozostałe zlokalizowane są w kilku miejscach, tworząc grupy od 2 do 5 studzienek. Hydrogramy wypływu wody ze studzienek podtopionych i wylewających przedstawiono na ryc. 5.4. Zgodnie z przyjętymi założeniami, maksymalne natężenie przepływu wynosi dla nich odpowiednio 0,1 i 0,2 m³/s.

27 Ryc. 5.3 Podtapiane i wylewające studzienki na obszarze zlewni Si­8 (Gmina Kielce

2010)

Ryc. 5.4 Hydrogramy wypływu wody ze studzienek

0 0,05 0,1 0,15 0,2 0,25

0 1000 2000 3000 4000 5000

Q [m3/s]

t [s]

podtopione wylewające

28 5.1 Wyniki symulacji spływu powierzchniowego w modelu uproszczonym

w programie PCRaster

Podczas każdego kroku obliczeniowego aktualizowana była informacja o kierunkach spływu. W pierwszym kroku czasowym mapa kierunków spływu stworzona została na podstawie NMT, a w każdym kolejnym powstawała na bazie NMT, do którego dodano wartości rzędnej wody z kroku poprzedniego. Mapa uzyskana na bazie NMT (w pierwszym kroku czasowym) zawierała 1243 obszary bezodpływowe – piksele o wartości 5 (ryc. 5.5).

Ryc. 5.5 Fragment mapy rastrowej przedstawiającej kierunki spływu

Ponieważ wyniki symulacji uproszczonego modelu spływu powierzchniowego były zapisywane co 10 s otrzymano w sumie 540 map rastrowych przedstawiających stan wody oraz tyle samo pokazujących wartości przepływu wody między komórkami. Daje to możliwość przeglądania wyników obliczeń w postaci animacji oraz obserwowania w jakiej kolejności zalewane zostają poszczególne tereny. Jednakże do celów porównawczych pomiędzy dwoma modelami posłużono się mapami rastrowymi pokazującymi maksymalne wartości stanu wody oraz przepływu dla obliczonej symulacji.

Maksymalna powierzchnia zalanego obszaru wyniosła 102890 m² (ryc. 5.6).

W większości miejsc głębokość zalewu wynosiła od 0,01 do 0,5 m (90% całego

29 zalanego obszaru). Maksymalna głębokość dochodziła do 3,15 m w lokalnych zagłębieniach terenu. Głębokość wody większa niż 1 m wystąpiła na niecałych 4%

powierzchni zalanego obszaru. W kilku pikselach zaobserwowano bardzo duże głębokości wody, wynoszące od kilku do kilkunastu metrów. Tak odstające wyniki, niemające odzwierciedlenia w morfologii terenu, są efektem niestabilności modelu.

Głębokość zalewu dla tych komórek została dostosowana do wartości z komórek sąsiednich. Średnia maksymalna głębokość zalewu dla całego obszaru wyniosła około 0,24 m. Wylewająca się z kanalizacji deszczowej woda zalała w większości sąsiednie tereny ulic. Podtopionych zostało 108 budynków, w tym niektóre zostały otoczone przez wodę ze wszystkich stron.

Ryc. 5.6 Wyniki symulacji modelu uproszczonego – maksymalna głębokość zalewu Duże wartości głębokości zalewu występowały w bliskim sąsiedztwie miejsc wypływu wody z sieci kanalizacji deszczowej (ryc. 5.7). Nie we wszystkich przypadkach

30 pokrywało się to z występowaniem lokalnych zagłębień terenu. Największy zasięg i głębokości zalewu wystąpiły w miejscach koncentracji studzienek podtapianych i wylewających – skrzyżowanie ul. Bodzentyńskiej i ul. Kościuszki, a także wzdłuż ul. Źródłowej.

Ryc. 5.7 Głębokość zalewu w miejscu największej koncentracji studzienek – model uproszczony

Maksymalne wartości przepływu w uproszczonym modelu spływu powierzchniowego wahały się od 0 do prawie 54 m³/s (ryc. 5.8). Prawie 60% zalanego obszaru charakteryzowało się przepływem maksymalnym mniejszym niż 0,5 m³/s. Z drugiej strony 39% pikseli otrzymało wartość większą niż 1 m³/s. Na bardzo niewielkiej części zalewanych terenów występowały przepływy maksymalne o wartościach od 0,5­1 m³/s.

31 Ryc. 5.8 Wyniki symulacji modelu uproszczonego – maksymalna wartość przepływu Bardzo wysokie wartości przepływu występowały na odcinkach łączących zalane obszary w postaci cienkiej strugi wody (ryc. 5.9). Istotny jest fakt, iż na tych odcinkach (o szerokości jednego piksela) nie została zarejestrowana głębokość wody na wynikowej mapie maksymalnych głębokości. Po dokładniejszej analizie danych okazało się, iż maksymalne głębokości wody w tych miejscach osiągają wartości od 1,0·10^­6 do 1,0·10^­16. Zbadanie wyników symulacji krok po kroku, pozwoliło na zaobserwowanie licznych niestabilności w uproszczonym modelu spływu powierzchniowego. Co któryś krok czasowy następowały skokowe wzrosty wartości przepływu i następowała wymiana wody między poszczególnymi obszarami zalanymi.

32 Ryc. 5.9 Przykładowe miejsce występowania bardzo wysokich wartości przepływu

w modelu uproszczonym

5.2 Wyniki symulacji spływu powierzchniowego w modelu ANUGA

Z powstałego pliku wynikowego wygenerowana informację o maksymalnej głębokości wody oraz maksymalnych wartościach przepływu. Dane te wyeksportowano do plików rastrowych.

Maksymalna powierzchnia zalanego obszaru wyniosła 150124 m² (ryc. 5.10).

Większość terenów (86%) podtopiona została do głębokości 0,5 m. Maksymalna głębokość dochodziła do 3,4 m w lokalnych zagłębieniach terenu przy ul. Warszawskiej oraz ul. Bohaterów Warszawy. Głębokość wody większa niż 1 m wystąpiła na niecałych 4% powierzchni zalanego obszaru. Średnia maksymalna głębokość zalewu dla całego obszaru wyniosła około 0,29 m. Wylewająca się z kanalizacji deszczowej woda spływała ulicami łącząc się w jedną strugę ostatecznie wpadającą do rzeki Silnicy.

Podtopionych zostało 120 budynków, w tym niektóre zostały otoczone przez wodę ze wszystkich stron.

Najgłębsze i najbardziej rozległe rozlewisko powstało u zbiegu ul. Bodzentyńskiej

33 i ul. Kościuszki, a także dalej na wschód w stronę rzeki Silnicy. Były to miejsca, w których gromadziła się nie tylko woda wydostająca się z sąsiednich studzienek, ale również woda spływająca z miejsc położonych w wyższych partiach zlewni.

Ryc. 5.10 Wyniki symulacji modelu ANUGA – maksymalna głębokość zalewu

Maksymalne wartości przepływu w uproszczonym modelu spływu powierzchniowego wahały się głównie w przedziale od 0 do 0,1 m³/s (ryc. 5.8). Nie całe 98% zalanego obszaru charakteryzowało się przepływem maksymalnym mniejszym niż 0,1 m³/s.

Znikoma była ilość komórek (372), w których zaobserwowano przepływ większy niż 0,2 m³/s. Maksymalna wartość przepływu wyniosła 1,16 m³/s, jednakże był to jedyny piksel, w którym przepływ przekroczył 1 m³/s.

Występowanie najwyższych maksymalnych wartości przepływu w większości zaobserwowane zostało wzdłuż ulic, czyli w miejscach charakteryzujących się najmniejszą szorstkością.

34 Ryc. 5.11 Wyniki symulacji modelu ANUGA – maksymalna wartość przepływu

5.3 Porównanie wyników symulacji z obydwu modeli

Obszary zalane zarówno w modelu uproszczonym jaki i w modelu ANUGA zajmują powierzchnie 63279 m², co stanowi odpowiedni 62% i 42% całości zalanych terenów w tych modelach. Powierzchnia zlewni, która została podtopiona w jednym bądź drugim modelu wyniosła 176553 m² (ryc. 5.12).

Obszary, które zostały zalane tylko w wyniku symulacji w programie ANUGA, pokrywają w dużej mierze odcinki ulic i przyległe do nich tereny, którymi woda spływała w kierunku ujścia do rzeki Silnicy. Teren osiedla i parku zlokalizowanego przy ul. Bohaterów Warszawy również został podtopiony tylko w tym modelu. Całkowita powierzchnia omawianych powyżej terenów wyniosła 73663 m².

Wyniki symulacji w modelu zbudowanym na bazie programu PCRaster w największym stopniu odbiegają od wyników uzyskanych w modelu ANUGA w miejscach

35 zlokalizowanych w sąsiedztwie podtapianych i wylewających studzienek.

Ma to odzwierciedlenie zarówno w zasięgu jak i głębokości zalewu. W modelu uproszczonym woda w dużo większym stopniu rozlewała się na tereny sąsiadujące ze studzienkami. Całkowita powierzchnia obszarów podtopionych tylko w tym modelu wyniosła 39611 m².

Ryc. 5.12 Porównanie zasięgu zalewanych obszarów w obu modelach

Średnia różnica maksymalnej głębokości zalewu pomiędzy obydwoma modelami wyniosła 0,28 m. Głębokości uzyskane z programu ANUGA były większe na 7%

obszaru wspólnie podtopionego przez modele, przy średniej różnicy równej 0,03 m (ryc. 5.13). Zalew o zbliżonej głębokości dla obu symulacji (różnica do 0,25 m) zanotowano na 66% powierzchni terenu. Różnica głębokości wody wynosząca ponad 1 m dotyczyła 6% uzyskanych wyników. Maksymalne wartości różnicy głębokości dochodziły odpowiednio do 0,64 m i 3,14 m dla obszarów o głębszym zalewie

36 w programie ANUGA niż w modelu uproszczonym i w sytuacji odwrotnej. Średnia wartość różnicy uzyskanych głębokości wyniosła 0,28 m. Największe różnice między wynikami obu symulacji zaobserwowano w sąsiedztwie studzienek oraz lokalnych zagłębień terenu.

Ryc. 5.13 Różnica głębokości wody pomiędzy wynikami uzyskanymi z obu modeli Wartości maksymalnych przepływów w 17% wspólnie zalanej powierzchni były wyższe w wynikach uzyskanych z programu ANUGA (ryc. 5.14). Średnia różnica wyniosła dla tego obszaru tylko 0,02 m³/s. Dla porównania ta sama wartość dla całego obszaru wyniosła 2,33 m³/s. Niewielkie różnice w maksymalnych wartościach przepływu między wynikami z obu modeli (do 0,1 m³/s) zaobserwowano na 33%

wspólnie zalanej powierzchni. Największe różnice między modelami zanotowano w bezpośrednim sąsiedztwie studzienek oraz w miejscach występowania dużych niestabilności w modelu uproszczonym (wąskie strugi wody łączące zalewane obszary).

37 Maksymalne wartości przepływów uzyskane z tego modelu były w ponad 50%

przypadków większe o ponad 1 m³/s niż te uzyskane z modelu ANUGA.

Ryc. 5.14 Różnica maksymalnych wartości przepływów pomiędzy wynikami uzyskanymi z obu modeli

38 6 Dyskusja

Przeprowadzone analizy wykazały, iż wyniki przeprowadzonych obliczeń symulacji spływu powierzchniowego w obu modelach znacznie się od siebie różnią. Zasięg zalewu uzyskany z uproszczonego modelu był znacznie mniejszy (o około 30%) od tego uzyskanego w modelu ANUGA. Dodatkowo w modelu uproszczonym, w przeciwieństwie do dwuwymiarowego modelu hydrodynamicznego, spływająca woda nie połączyła się w jedną strugę płynącą do najniższego punktu, czyli rzeki Silnicy. Jednakże na większości obszaru (93%), który został zalany w obu symulacjach, większe głębokości występowały w modelu bazującym na programie PCRaster.

Szczególnie wyraźne różnice występowały w sąsiedztwie studzienek, gdzie zasięg i głębokość zalewu były większe w przypadku modelu uproszczonego. Jeszcze większe różnice można było zaobserwować pomiędzy uzyskanymi wartościami przepływów, gdzie bardzo wyraźnie ujawniły się niestabilności występujące w uproszczonym modelu.

Uproszczony model spływu powierzchniowego utworzony został na bazie programu PCRaster, który był również podstawą dla modelu LISPFLOOD­FP (Bates, De Roo 2000). LISPFLOOD­FP w dużej mierze wykorzystywany był do modelowania przepływu wody w rzekach oraz na terasach zalewowych. Obliczenia te dokonywane były zazwyczaj na rastrach o średniej bądź dużej rozdzielczości (Bates, De Roo 2000, Horit, Bates 2001, Aronica i in. 2002). W badaniach tych, zasięg wygenerowanej strefy zalewowej był zgodny w 70­90% z danymi rzeczywistymi. Model ten testowany był również do symulacji powodzi miejskich (Hunter i in. 2008), gdzie uzyskano wyniki porównywalne z innymi modelami dwuwymiarowymi. Duże rozbieżności między wynikami z obu porównywanych w opracowaniu modeli wynikają w większości z niestabilności jakie występowały w modelu uproszczonym. W celu jego udoskonalenia wprowadzony powinien zostać mechanizm, kontrolujący bilans masy wody, która znajduje się w modelu w każdym kroku obliczeniowym. Pozwoliłoby to na uniknięcie sytuacji, w której w wyniku błędów obliczeniowych z danego piksela odpływa więcej wody niż jest w nim zawarte. Mechanizm taki zastosowany został w wspomnianym wyżej modelu LISPFLOOD­FP. Znaczenie tego mechanizmu do redukcji niestabilności w tego typu modelach było już wielokrotnie odnotowywane (Bates, De Roo 2000, Hunter i in. 2008). Bates i DeRoo (2000) zaobserwowali również zwiększanie się błędów związanych z bilansem masy wody wraz ze zmniejszaniem

39 wielkości piksela i zwiększaniem się ilości komórek obliczeniowych. Poprawę stabilności mogłoby przynieść więc zwiększenie rozdzielczości NMT użytego do modelowania, a także wygładzenie w miejscach, gdzie następuje nagła zmiana wysokości z piksela na piksel. Modele opisujące ruch wody w sposób uproszczony mogą dawać niewiarygodne rezultaty w miejscach, gdzie występują duże zmiany wysokości pomiędzy pikselami oraz nagłe zmiany prędkości przepływu (Mark i in.

2004). W przypadku programu ANUGA nie zaobserwowano znaczących niestabilności.

Jest to w dużej mierze wynikiem zastosowania bardziej złożonego mechanizmu obliczeniowego. Dużą rolę odgrywa również mechanizm automatycznego ustalania kroku czasowego obliczeń na podstawie wielkości trójkąta oraz prędkości przepływu wody.

Dynamiczne modele symulacyjne środowiska przyrodniczego powinny być kalibrowane na rzeczywistych danych w celu odpowiedniego dopasowania parametrów modelu i uzyskania wiarygodnych wyników. Wymagane jest również przeprowadzenie weryfikacji modelu na podstawie innych danych historycznych niż tych użytych w procesie kalibracji (Longley i in. 2008). Problem braku danych pomiarowych do kalibracji i weryfikacji modelu powodzi miejskich był bardzo często poruszany podczas analizy i oceny jego wyników (Mark i in. 2004, Aronica, Lanza 2005, Hunter i in.

2008). Właściwie nie możliwe jest pozyskanie danych pomiarowych o natężeniu przepływu wody spływającej powierzchniowo. Z powyższych względów najczęściej do kalibracji wykorzystuje się informację o zasięgu i głębokości zalewu oraz wartości natężenia opadu dla tego zdarzenia (Mark i in. 2004). Problemem pozostaje jednak lokalny charakter zjawiska, gdzie dane opadowe pozyskane ze stacji oddalonej nawet o kilka kilometrów od miejsca zdarzenia mogą okazać się nieprzydatne.

Brak danych pomiarowych nie pozwala na jednoznaczne stwierdzenie, który model lepiej opisuje zjawisko spływu powierzchniowego w mieście. Problem ten był wielokrotnie zauważany przy próbach porównania modeli hydraulicznych (Lhomme i in. 2006). Zaobserwowane w uproszczonym modelu znaczące niestabilności powodują, iż wyniki w nim uzyskane obarczone są dużym błędem. Z powyższych względów model hydrodynamiczny ANUGA w większym stopniu pozwala na zbadanie zjawiska spływu powierzchniowego.

Należy podkreślić, że modelowanie spływu powierzchniowego w mieście jest utrudnione ze względu na brak naturalnych ścieżek spływu oraz skomplikowaną

40 topografię. Żaden model nie jest w stanie dokładnie odzwierciedlić rzeczywistości i uwzględnić wszystkich czynników wpływających na dane zjawisko, jednakże stanowi niezwykle przydatne (choć uproszczone) narzędzie pozwalające na prognozowanie skutków powodzi miejskich. Zrozumienie i poznanie źródeł niepewności stosowanych modeli jest więc niezmiernie ważne, zwłaszcza jeżeli na ich podstawie podejmowane są decyzje wpływające na funkcjonowanie miasta (Hapuarachchi i in. 2011).

41 7 Wnioski

Wysoka rozdzielczość przestrzenna danych wymusza wysoką rozdzielczość czasową prowadzonych obliczeń. Wykorzystanie NMT o wysokiej rozdzielczości przestrzennej do modelowania hydraulicznego powoduje konieczność przeprowadzenia obliczeń z bardzo małym krokiem czasowym. Mimo użycia bardzo małego kroku czasowego obliczeń (0,1 s) nie udało się osiągnąć zadowalającej stabilności w uproszczonym modelu spływu powierzchniowego, co prawdopodobnie spowodowane było brakiem mechanizmu kontrolującego bilans masy wody, która znajdowała się w modelu w każdym kroku obliczeniowym.

Brak danych pomiarowych uniemożliwia jednoznaczne stwierdzenie, który model daje rezultaty bardziej zbliżone do rzeczywistości. Jednakże zaobserwowane niestabilności w uproszczonym modelu skłaniają do stwierdzenia, iż do symulacji spływu powierzchniowego na terenach miejskich bardziej odpowiedni jest model ANUGA.

Dodatkowo, brak danych pomiarowych nie pozwala na przeprowadzenie procesu kalibracji i weryfikacji modelu. Jest to jeden z głównych powodów, dla którego wyniki takich symulacji powinny być traktowane z dużą ostrożnością.

Do opisu spływu powierzchniowego w terenie o skomplikowanej topografii jaka występuje w mieści wymagane jest użycie modelu hydrodynamicznego, który pozwala na obliczenia w warunkach szybko zmieniającego się poziomu wody oraz nagłych zmian prędkości przepływu.

42 8 Literatura

1. Aronica G, Bates PD, Horritt MS (2002) Assessing the uncertainty in distributed model predictions using observed binary pattern information within GLUE.

Hydrological Processes 16: 2001­2016.

2. Aronica G, Lanza LG (2005) Drainage efficiency in urban areas: a case study.

Hydrological Processes 19: 1105­1119.

3. Bates PD, De Roo APJ, (2000) A simple raster­based model for flood inundation simulation. Journal of Hydrology 236: 54­77.

4. Boonya­Aroonnet S, Maksimović Č, Prodanović D, Djordjević S (2007) Urban pluvial flooding: Development of GIS based path­ways model for surface flooding and interface with surcharged sewer model. Proc., 6th NOVATECH International Conf., Workshop I, Graie, Lyon, France, 481–488.

5. Chen AS, Hsu MH, Chen TS, Chang TJ (2005) An integrated inundation model for highly developed urban areas. Water Science and Technology 51: 221–229.

6. Chow VT (1959) Open channel hydraulics. McGraw­Hill Co., New York.

7. Ciupa T (2009) Wpływ zagospodarowania terenu na odpływ i transport fluwialny w małych zlewniach na przykładzie Sufragańca i Silnicy (Kielce).

Wydawnictwo Uniwersytetu Humanistyczno­Przyrodniczego Jana Kochanowskiego, Kielce.

8. Courant R, Friedrichs K, Lewy H (1967) On the partial difference equations of mathematical physics. IBM Journal of Research and Development 11: 215–234.

9. ESRI (2012) ArcGIS Desktop Help 10.1. Dostępny pod adresem:

http://resources.arcgis.com/en/help/main/10.1/, sprawdzony: wrzesień 2013.

10. European Environment Agency (2011) Mapping Guide for a European Urban Atlas. Dostępny pod adresem: http://www.eea.europa.eu/data­and­

maps/data/urban­atlas#tab­methodology, sprawdzony: sierpień 2013.

11. European Environment Agency (2013) Urban Atlas. Dostępny pod adresem:

http://www.eea.europa.eu/data­and­maps/data/urban­atlas, sprawdzony: sierpień 2013.

43 12. Geoscience Australia (2010) ANUGA User Manual, Relase 1.2.1. Dostępny pod

adresem: http://anuga.anu.edu.au/, sprawdzony: sierpień 2013.

13. Gmina Kielce (2010) Koncepcja zagospodarowania wód deszczowych dla Miasta Kielce. Dostępny pod adresem:

http://www.um.kielce.pl/gfx/kielce2/files/srodowisko/zagosp_wod_deszcz.pdf, sprawdzony: czerwiec 2013.

14. Gmina Kielce (2012) Wykonanie koncepcji ochrony przeciwpowodziowej w zlewni Bobrzy, Silnicy, Sufragańca i Lubrzanki na obszarze miasta Kielce.

Dostępne pod adresem: http://www.um.kielce.pl/srodowisko/art,5409,koncepcja­

ochrony­przeciwpowodziowej­w­zlewni­bobrzy­silnicy­sufraganca­i­

lubrzanki.html, sprawdzony: czerwiec 2013.

15. GUGIK (2013) Projekt ISOK. Dostępny pod adresem:

http://www.gugik.gov.pl/projekty/isok, sprawdzony: czerwiec 2013.

16. Hapuarachchi HAP., Wang QJ, Pagano TC (2011) A review of advances in flash flood forecasting. Hydrological Processes 25: 2771–2784.

17. Horritt MS, Bates PD (2001) Effects of spatial resolution on a raster based model of flood flow. Journal of Hydrology 253: 239­249.

18. Hunter NM, Bates PD, Neelz S, Pender G, Villanueva I, Wright NG, Liang D, Falconer RA, Lin B, Waller S, Crossley AJ, Mason DC (2008) Benchmarking 2D hydraulic models for urban flooding. Proceedings of the ICE ­ Water Management 161: 13­30.

19. ISOK (2013) Informatyczny System Osłony Kraju. Dostępny pod adresem:

http://isok.imgw.pl/, sprawdzony: październik 2013.

20. Klimaszewski M (1972) Geomorfologia Polski Tom 1 Polska Południowa Góry i Wyżyny. PWN, Warszawa.

21. Kondracki J (2000) Geografia regionalna Polski. Wydawnictwo Naukowe PWN, Warszawa.

22. Lan­Anh N, Boxall J, Saul A, Willems P (2012) Integrated river flow modelling:

a case study. Urban Water Journal 9: 1­18.

23. Leandro J, Chen A, Djordjević S, Savić D (2009) Comparison of 1D/1D and

44 1D/2D Coupled (Sewer/Surface) Hydraulic Models for Urban Flood Simulation, Journal of Hydraulic Engineering 135: 495–504.

24. Lhomme J, Bouvier C, Mignot E, Paquier A (2006) One­dimensional GIS­based model compared to two­dimensional model in urban floods simulations. Water Science and Technology 54: 83–91.

25. Lhomme J, Bouvier C, Perrin JL (2004) Applying a GIS­based

geomorphological routing model in urban catchments. Journal of Hydrology 299: 203­216.

26. Longley PA, Goodchild MF, Maguire DJ, Rhind DW (2008) GIS Teoria i praktyka. Wydawnictwo Naukowe PWN, Warszawa.

27. Majewski W (2008) Urban flash flood in Gdańsk, 2001. Annal of Warsaw University of Life Science – SGGW, Land Reclamation 39: 129­137.

28. Mark O, Weesakul S, Apirumanekul C, Aroonnet SB, Djordjević S (2004) Potential and limitations of 1D modelling of urban flooding. Journal of Hydrology 299: 284–299.

29. PCRaster Team (2011) PCRaster Documentation, Release 3.0.1. Dostępny pod adresem: http://pcraster.geo.uu.nl/downloads/, sprawdzony: wrzesień 2013.

30. Puzyrewski R, Sawicki J (1998) Podstawy mechaniki płynów i hydrauliki.

Wydawnictwo Naukowe PWN, Warszawa.

31. Roberts SG, Nielsen OM, Jakeman J (2008) Simulation of tsunami and flash floods. Modeling Simulation and Optimization of Complex Processes:

Proceedings of the Third International Conference on High Performance Scientific Computing, March 6–10, 2006, Hanoi, Vietnam: 489­499.

32. Schmitt TG, Thomas M, Ettrich N (2004) Analysis and modeling of flooding in urban drainage systems. Journal of Hydrology 299: 300–311.

33. Szymkiewicz R (2000) Modelowanie matematyczne przepływów w kanałach otwartych. Wydawnictwo Naukowe PWN, Warszawa.

34. Zevenbergen C, Cashman A, Evelpidou N, Pasche E, Garvin S, Ashley R (2011) Urban flood management. CRC Press.

45 Spis rycin

Ryc. 2.1 Obszar badań na tle mapy topograficznej ... 7

Ryc. 2.2 Ukształtowanie terenu badanej zlewni ... 9

Ryc. 2.3 Użytkowanie ziemi w zlewni Si­8, na podstawie European Urban Atlas (2013) ...11

Ryc. 2.4 Sieć kanalizacji deszczowej w zlewni SI­8 (Gmina Kielce 2010) ... 12

Ryc. 4.1 Wartości pikseli dla poszczególnych kierunków płynięcia wody ... 18

Ryc. 4.2 Schemat obszaru bezodpływowego oraz jego parametrów (źródło: PCRaster 2011) ... 19

Ryc. 4.3 Schemat obliczeń w uproszczonym modelu spływu powierzchniowego ... 21

Ryc. 4.4 Siatka trójkątów w metodzie objętości skończonych (źródło: Roberts i in. 2008) ... 22

Ryc. 4.5 Schemat działania funkcji Combine (ESRI 2012) ... 24

Ryc. 5.1 NMT z wprowadzoną informacją o wysokości budynków ... 25

Ryc. 5.2 Mapa wartości współczynnika szorstkości Manninga n dla obszaru objętego modelowaniem ... 26

Ryc. 5.3 Podtapiane i wylewające studzienki na obszarze zlewni Si­8 (Gmina Kielce 2010) ... 27

Ryc. 5.4 Hydrogramy wypływu wody ze studzienek... 27

Ryc. 5.5 Fragment mapy rastrowej przedstawiającej kierunki spływu ... 28

Ryc. 5.6 Wyniki symulacji modelu uproszczonego – maksymalna głębokość zalewu 29 Ryc. 5.7 Głębokość zalewu w miejscu największej koncentracji studzienek – model uproszczony ... 30

Ryc. 5.8 Wyniki symulacji modelu uproszczonego – maksymalna wartość przepływu ... 31

Ryc. 5.9 Przykładowe miejsce występowania bardzo wysokich wartości przepływu w modelu uproszczonym ... 32

Ryc. 5.10Wyniki symulacji modelu ANUGA – maksymalna głębokość zalewu ... 33

Ryc. 5.11Wyniki symulacji modelu ANUGA – maksymalna wartość przepływu ... 34

Ryc. 5.12Porównanie zasięgu zalewanych obszarów w obu modelach ... 35

Ryc. 5.13Różnica głębokości wody pomiędzy wynikami uzyskanymi z obu modeli ... 36

Ryc. 5.14Różnica maksymalnych wartości przepływów pomiędzy wynikami uzyskanymi z obu modeli ... 37

Spis tabel Tab. 3.1 Hierarchiczny podział na klasy użytkowania ziemi w bazie danych European Urban Atlas (European… 2011) ... 14

Tab. 4.1 Wartości współczynnika szorstkości Manninga dla poszczególnych klas użytkowania ziemi (Chow 1959) ... 15