AnAlizA wrAżliwości w modelowAniu hydrogeochemicznym nA przykłAdzie systemów wodA–skAłA–gAz
sensitivity AnAlysis in hydrogeochemicAl modelling on the exAmple of wAter–rock–gAs systems
Krzysztof Labus1, Katarzyna suchodoLsKa1
Abstrakt. Przedstawiono przykład analizy wrażliwości rezultatów otrzymywanych w modelowaniu systemów woda–skała–gaz do celów sekwestracji gazów kwaśnych (CO2 i H2S) na zmiany parametrów kinetyki reakcji. Z rezultatami stanowiącymi poziom odniesienia porównywano wyniki uzyskane po uwzględnieniu stałych tempa reakcji dla skalenia potasowego i klinochloru, podwyższonych lub po- mniejszonych o jeden rząd wielkości. Różnice w zidentyfikowanych efektach petroteksturalnych i fizykochemicznych są pozornie niewiel- kie, gdyż wynikają z modyfikacji pojedynczego tylko parametru. W rzeczywistości, przy występowaniu wielu parametrów opisujących modelowany system, efekty z nimi związane mogą ze sobą interferować.
Słowa kluczowe: analiza wrażliwości, modelowanie hydrogeochemiczne, interakcje woda–skała–gaz, sekwestracja geologiczna.
Abstract. The paper presents an example of sensitivity analysis of results obtained in modelling of water–rock–gas systems, for the purpose of acid gas sequestration (CO2 and H2S), for changes in reaction kinetics parameters. The reference level was compared to the results, including the elevated or reduced reaction rate constants of K-feldspar and clinochlore. Differences in the identified petrotextural and physicochemical effects are seemingly insignificant, as they result from the modification of a single parameter only. In fact, when there are many parameters describing the modelled system, the associated effects may interfere with each other.
Key words: sensitivity analysis, hydrogeochemical modelling, water–rock–gas interactions, geological sequestration.
1 Politechnika Śląska, Wydział Górnictwa, Inżynierii Bezpieczeństwa i Automatyki Przemysłowej, ul. Akademicka 2; 44-100 Gliwice; e-mail:
krzysztof.labus@polsl.pl, katarzyna.suchodolska@polsl.pl.
WStęp
W przypadku stosowania modelowania komputerowego do symulacji złożonych procesów, jakimi są interakcje w sys- temach hydrogeochemicznych, analiza wrażliwości okazuje się odpowiednim narzędziem dla oceny wpływu zróżnicowa- nia wielkości parametrów na uzyskiwane rezultaty. Jest ona zatem istotna z punktu widzenia weryfikacji wyników modelu i pozwala ponadto na określenie przydatności oprogramowa- nia zastosowanego do danego celu obliczeniowego.
Pod pojęciem analiza wrażliwości zwykle rozumie się badanie zmian rezultatu pod wpływem zmiany danych wej- ściowych. Działania takie mogą przynieść kilka następują- cych korzyści:
– możliwa jest weryfikacja poprawności wykonanych wcześniej prognoz eksperckich;
– istnieje możliwość ujawnienia mało istotnych zmien- nych lub niepotrzebnych komplikacji modelu;
– dane wejściowe mogą zostać uszeregowane pod wzglę- dem znaczenia ich wpływu na wynik.
W ostatnich dwóch przypadkach liczba zmiennych może zostać znacznie zmniejszona, dzięki czemu, zwłaszcza dla złożonych kodów, skrócony zostaje czas obliczeń a ponadto ewentualna analiza niepewności może być przeprowadzona prościej i szybciej (Ekberg, 1999). Z drugiej strony istnieją sy- tuacje, w których ten rodzaj przesiewu parametrów może spo- wodować, że niektóre ważne procesy, np. efekty sprzężone, mogą pozostać niewykryte podczas prowadzonych symulacji.
W modelowaniu hydrogeochemicznym dróg reakcji w układach woda–skała–gaz najczęściej uwzględniane są parametry opisujące każdy z elementów takiego systemu (fig. 1).
Objętość wód porowych w modelowanym systemie może być określona na podstawie pomiaru porowatości efektywnej skały, np. metodą porozymetryczną, przy zało- żeniu pełnego lub częściowego nasycenia porów; parametry fizykochemiczne i temperatura przyjmowane są na podsta- wie bezpośrednich oznaczeń i analiz (często są to dane ar- chiwalne).
W przypadku mieszanin gazów najczęściej zakłada się ich skład, wyrażony w postaci ciśnień parcjalnych.
Skała opisywana jest przez skład fazowy szkieletu ziar- nowego, określany metodami mikroskopowymi lub za po- mocą dyfraktometrii rentgenowskiej. Powierzchnia właści- wa skały może zostać zmierzona metodami adsorpcyjnymi, zaś reaktywna powierzchnia właściwa minerałów może być oszacowana na podstawie określenia wymiarów, kształtu i morfologii ziaren przy pomocy technik mikroskopowych.
Tempo (stała szybkości reakcji) rozpuszczania i krystali- zacji dla poszczególnych faz mineralnych jest na ogół przyj- mowana na podstawie danych literaturowych, które naj- częściej różnią się od siebie, w zależności od zastosowanej metody badawczej i nie zawsze uwzględniają zróżnicowanie warunków pH środowiska, w którym reakcje te przebiegają.
Szybkość rozpuszczania i krystalizacji minerałów jest iloczynem ich reaktywnej powierzchni właściwej i stałej ki- netyki reakcji. Pośród wymienionych wcześniej parametrów
zwłaszcza te właśnie wielkości obarczone są wysoką nie- pewnością i obejmują szeroki zakres wartości. W związku z tym, na potrzeby niniejszej pracy, przeprowadzono szereg symulacji (analiz) wrażliwości, skupiając się na wpływie zróżnicowania stałych kinetyki reakcji dla analizowanych faz mineralnych.
DAne WejścioWe i metoDyKA moDeloWAniA
W niniejszej pracy zaprezentowano przykład analizy wrażliwości na zmiany parametrów kinetyki reakcji, w od- niesieniu do rezultatów otrzymywanych w modelowaniu hydrogeochemicznym systemów woda–skała–gaz w celu sekwestracji gazów kwaśnych (mieszaniny dwutlenku wę- gla i siarkowodoru) w poziomach wodonośnych. Badania prowadzone nad zagadnieniami interakcji gazów kwaśnych ze skałami wynikają z przekonania o konieczności ograni- czenia efektu cieplarnianego (np. Holloway, 2005; Xu i in., 2007), oraz lokowania gazów odpadowych po odsiarczaniu gazu ziemnego (Lubaś, Szott, 2010). Wynikają one tak- że z konieczności intensyfikacji wydobycia ropy naftowej i gazu lub zastosowaniem szczelinujących cieczy energety- zujących (Sinal, Lancaster, 1987; Labus i in., 2016).
Obliczenia przedstawiono na podstawie próbki skały o znacznej przepuszczalności, reprezentującej formację dę- bowiecką (miocen) z rejonu Górnośląskiego Zagłębia Wę- glowego (GZW). Charakteryzuje się ona strukturą psefitową i słabym wysortowaniem. W składzie obecne są różnie obto- czone okruchy drobnokrystalicznych skał krzemionkowych, węglanowych i krzemionkowo-ilastych. Spoiwo węglano- wo-ilaste stanowi około 15% objętości skały. Skład mine- ralogiczny próbki skał określono metodą dyfraktometryczną (tab. 1); do identyfikacji faz mineralnych wykorzystano tak- że mikroskop elektronowy z analizatorem EDX. Parametry przestrzeni porowej badanej próbki scharakteryzowano me- todą porozymetrii rtęciowej.
Modele interakcji w systemie woda–skała–gaz zostały wykonane w dwóch etapach, zgodnie z metodyką opisaną w poprzednich pracach (m.in. Labus i in., 2015; Labus, Su- chodolska, 2015). Etap pierwszy pozwalał na analizę zmian zachodzących w poziomach wodonośnych pod wpływem iniekcji CO2 i/lub H2S. Etap drugi umożliwiał prześledzenie długotrwałych efektów sekwestracji, w tym zmian porowa- tości skał oraz sekwestracji mineralnej C i/lub S. W niniej- szym artykule omówiono wyłącznie przykład dotyczący tego właśnie etapu.
Zastosowany w symulacjach skład wód porowych usta- lono na podstawie modelowania równowagowego solanki złożowej o składzie podanym w tabeli 2, z minerałami bu- dującymi szkielet ziarnowy modelowanej skały. Ciśnienie zatłaczania mieszaniny CO2 i H2S zadano wg reżimu hy- drostatycznego i a jego wartości wprowadzono jako aktyw- ność ciśnieniową, zgodnie z wymaganiami oprogramowania zastosowanego w modelu – Geochemist’s Workbench 9.0 (GWB) (Bethke, 2008; Bethke, Yeakel, 2013). Temperatu- Fig. 1. najczęściej stosowane parametry opisujące systemy
woda–skała–gaz w modelach hydrogeochemicznych The most commonly used parameters describing water–rock–gas
systems in hydrogeochemical models objętość volume Woda porowa
Pore water parametry fizykochemiczne
physicochemical parameters temperatura
temperature
GazGas ciśnienia parcjalne składników
partial pressure skład mineralny mineral composition Skała
Rock powierzchnia właściwa
specific surface tempo reakcji reaction rate
rę oszacowano na podstawie bezpośrednich pomiarów i da- nych archiwalnych (tab. 2). W obliczeniach wykorzystano uproszczone równanie kinetyki rozpuszczania/krystaliza- cji (Lasaga, 1984), zgodnie z którym minerały krystalizują w momencie, gdy roztwór jest nimi przesycony, a ulegają rozpuszczaniu, gdy jest niedosycony:
gdzie:
rk – szybkość reakcji [mol·s–1], rk > 0 – rozpuszczanie, rk < 0 – krystalizacja,
AS – reaktywna powierzchnia minerału [cm2], kT – stała reakcji [mol·cm–2·s–1],
Q – iloczyn aktywności,
K – stała równowagi reakcji rozpuszczania.
Powierzchnię właściwą ziaren minerałów obliczono przy założeniu kulistości ziaren o określonej średnicy na podsta- wie modelu sferycznego, dla znanych mas cząsteczkowych oraz objętości molowych każdego z minerałów:
gdzie:
A – powierzchnia sfery o danej średnicy, v – objętość molowa minerału,
V – objętość sfery o danej średnicy, MW – masa molowa minerału.
Wykorzystane w modelowaniu stałe szybkości reakcji minerałów (kT) przyjęto na podstawie danych literaturowych (Palandri, Kharaka, 2004) – tabela 1.
W celu przeprowadzenia analizy wrażliwości rezultatów modelowania systemu woda–skała–gaz na zmiany parame- trów kinetyki reakcji, w przykładowym poziomie wodono- śnym, podczas symulowanej sekwestracji gazów, wykonano odpowiednie obliczenia. Uwzględniły one wpływ zróżnico- wania stałej szybkości reakcji dla skalenia potasowego oraz klinochloru na przebieg reakcji rozkładu i krystalizacji faz mineralnych w trakcie modelowanego okresu składowania gazów kwaśnych oraz efekty petroteksturalne (zmiany poro- watości) i wielkość mineralnej sekwestracji pierwiastka węgla w formacji skalnej. Przyjęto założenie, że rezultaty stanowią- ce poziom odniesienia (scenariusz bazowy) wynikają z wyko- rzystania w modelowaniu wielkości stałych szybkości reakcji zaczerpniętych z literatury. Z tymi rezultatami porównywano wyniki uzyskane dla pozostałych dwóch scenariuszy, odpo- wiadających uwzględnieniu stałych szybkości reakcji pod- wyższonych lub pomniejszonych o jeden rząd wielkości.
W przeprowadzonej analizie wrażliwości uwzględniono klinochlor (o udziale w matrycy skały na poziomie 10%), któ- rego rozkład powoduje uwolnienie znacznej ilości jonów Mg2+, umożliwiając mineralną sekwestrację CO2 w krystalizującym, wtórnym dolomicie. Z drugiej strony, aby wykazać, że istotne zmiany hydrochemiczne nie zależą wyłącznie od udziału obję-
tabela 1 parametry modelu, skład mineralogiczny matrycy skalnej
oraz powierzchnia właściwa minerałów i logarytm stałej szybkości reakcji uwzględnione w modelowaniu Model parameters, mineral composition of the rock matrix:
specific surface and logarithm of the reaction rate constant of minerals, applied in modelling
Porowatość [–] 0,224
Przepuszczal-
ność [mD] 205
Objętość modelu
[cm3] 10000
Minerały Skład Udział
[% vol.]
Pow.
właś.
[cm2/g] log kT
Kwarc SiO2 34,0 22,66 –13,99
Muskowit KAl3Si3O10(OH)2 22,0 105,98 –11,85
Syderyt FeCO3 0,2 14,83 –3,76
Albit NaAlSi3O8 6,8 22,94 –9,87
Klinochlor Mg5Al2Si3O10(OH)8 10,0 1117,91 –9,79
K-Skaleń KAlSi3O8 2,0 23,47 –10,06
Kalcyt CaCO3 9,0 22,14 –0,3
Ankeryt
Fe 0,55 Ca Fe(CO3)2 5,0 20,94 –3,76
flogopit KAlMg3Si3O10(OH)2 11,0 1076,02 –12,40 kT – stała szybkości reakcji / reaction rate constant [mol/m2/s]
tabela 2 pierwotny skład wód porowych oraz skład mieszaniny gazów
Primary composition of pore water and composition of the gas mixture
Parametr Wartość
Temperatura [°C] 30
pH 7,8
Cl– [mol/kg] 2,831
SO4 2– [mol/kg] 7,8·10–4
HCO3 – [mol/kg] 0,002
Ca2+ [mol/kg] 0,939
Mg2+ [mol/kg] 0,053
Na+ [mol/kg] 0,842
K+ [mol/kg] 0,005
Ciśnienie [bar] 75
Aktywność ciśnieniowa [bar] CO2 H2S
41,8 14,5
AS
rk= kT (1− )Q
K [1]
[2]
SSM V= AMWv
tościowego danego składnika w matrycy skalnej analizowano także skaleń potasowy, mimo że jest go w analizowanym syste- mie względnie niewiele – 2%. W przypadku wyższych propor- cji H2S do CO2 w mieszaninie gazów, wśród minerałów wtór- nych powstających w analizowanym systemie można byłoby spodziewać się powstawania minerałów siarczkowych.
WyniKi moDeloWAniA
Na figurze 2A przedstawiono zróżnicowanie zmian po- rowatości próbki na etapie składowania CO2, w zależności od stałej szybkości reakcji przyjętej do modelowania dla skalenia potasowego. O ile w okresie do około 1000 lat tendencja zmian jest dla wszystkich scenariuszy podob- na – gwałtowny wzrost porowatości, po którym nastę- puje jej intensywny spadek a następnie szybki wzrost, to w miarę upływu modelowanego czasu zaznaczają się wyraźne różnice. Porowatość obliczona dla przypadku bazowego przejawia stabilny, powolny trend wzrosto- wy do wartości niemal 0,247. Dla przypadku powięk- szonego tempa reakcji widoczny jest niewielki spadek (z minimum około 1500 lat na skali czasu), a następnie wzrost i stabilizacja od około 4000 do 5500 lat, po której na- stępuje spadek porowatości do wartości 0,247. Przy obniżo- nym tempie reakcji wzrost porowatości trwa nadal, do czasu około 4200 lat, a następnie zaznacza się liniowy spadek do wartości poniżej 0,246. Końcowe porowatości obliczone dla modelowanego etapu składowania CO2 dla wszystkich sce- nariuszy są do siebie zbliżone, różnice pomiędzy nimi nie przekraczają 0,1 punktu procentowego (fig. 2A).
W przypadku klinochloru, w okresie do około 1800 lat modelowanego etapu, tendencje zmian porowatości dla wszystkich scenariuszy są do siebie zbliżone – fig. 2B. Na po- czątku modelowanego etapu zachodzi gwałtowny wzrost po- rowatości (najintensywniejszy dla scenariusza powiększonej stałej reakcji), po nim następuje spadek, a następnie wzrost porowatości. Dalej, w miarę upływu modelowanego czasu, zaznaczają się pewne różnice. Porowatość obliczona dla przy- padku bazowego przejawia stabilny, powolny trend wzrosto- wy do wartości około 0,2465. Dla przypadku powiększonej stałej szybkości reakcji zauważalny jest spadek (z minimum około 6000 lat na skali czasu), a następnie niewielki wzrost porowatości do wartości nieco powyżej 0,25. Przy obniżo- nej stałej szybkości reakcji porowatość wzrasta nieznacznie i w konsekwencji nie przekracza wartości 0,245. Różnice pomiędzy skrajnymi, końcowymi wartościami porowatości obliczonymi dla kolejnych scenariuszy modelowanego etapu składowania CO2 sięgają 0,5 punktu procentowego (fig. 2B).
Dla wyjaśnienia przyczyn różnic pomiędzy efektami opisanych powyżej zjawisk w analizowanych systemach, porównano przebieg reakcji dla scenariusza bazowego oraz pozostałych dwóch scenariuszy, w których uwzględniono stałe szybkości reakcji skalenia potasowego i klinochloru podwyższone lub pomniejszonych o jeden rząd wielkości.
Przykładowe wyniki symulacji przedstawiono na figurze 3.
Dla scenariusza bazowego (fig. 3A, 3B), w początko- wym okresie, trwającym do około 1 roku, zmiany porowato- ści próbki, w warunkach symulowanej sekwestracji CO2, są związane z rozkładem klinochloru i kalcytu oraz powstawa- niem dolomitu i kaolinitu, przy czym jony glinu pobierane są z roztworu wód porowych, który został w nie wzbogaco- ny na etapie iniekcji gazów:
Mg5Al2Si3O10(OH)8 + Al3++ 5CaCO3 + 7HCO3– + 2H2S = 2,5H2O + +1,5Al2Si2O5(OH)4 + 5CaMg(CO3)2 + 2CH4 + 2SO42–
Następnie, do około 2 roku, zauważalny jest rozkład no- wopowstałego kaolinitu, co sprzyja krystalizacji gibbsytu oraz saponitu-Ca, a krystalizacja dolomitu zostaje wyhamo- wana:
Mg5Al2Si3O10(OH)8 + 1,2H2O + 1,4Al2Si2O5(OH)4 + 0,52CaCO3 =
= 1,6Ca0,16Mg3Al0,33Si3,67O10(OH)2 + 4,3Al(OH)3 + 0,26CaMg(CO3)2
Kolejny etap, trwający do czasu około 100 lat, obejmu- je rozkład klinochloru, albitu oraz degradację saponitu-Ca i niewielkich ilości dolomitu, co pozwala na kontynuację krystalizacji gibbsytu, wytrącanie się saponitu-Na oraz re- krystalizację kalcytu:
Mg5Al2Si3O10(OH)8 + 1,6H2O + 1,7Ca0,16Mg3Al0,33Si3,67O10(OH)2 + + 0,3CaMg(CO3)2 + 1,2NaAlSi3O8 = 3,5Na0,33Mg3Al0,33Si3,67O10(OH)2 + 2,6Al(OH)3 + 0,6CaCO3
Etap końcowy obejmuje kilka głównych reakcji przebie- gających z różną intensywnością, w efekcie których noto- wany jest rozkład gibbsytu, flogopitu, skalenia potasowego i albitu oraz krystalizacja muskowitu i saponitu-Na:
2H2O + 2NaAlSi3O8 + H2S + KAlSi3O8 +3CO2 = KAl3Si3O10(OH)2 +
+ CH4 +2Na+ + 6SiO2+ SO42–
2Al(OH)3 + KAlMg3Si3O10(OH)2 + 3H2S + 3CO2 = KAl3Si3O10(OH)2
+ 3CH4 +3Mg2+ + 2SO42–
1,1Al(OH)3 + 1,2NaAlSi3O8 + 0,4H2S + KAlMg3Si3O10(OH)2 + + 0,4CO2 = Na0,33Mg3Al0,33Si3,67O10(OH)2 + 0,2H2O + 0,9Na+ + + KAl3Si3O10(OH)2 + 0,4SO42– + 0,4CH4
[3]
klinochlor kaolinit kalcyt
saponit-Ca gibbsyt dolomit [4]
klinochlor kalcyt
kaolinit dolomit
klinochlor saponit-Ca
dolomit albit saponit-Na
gibbsyt kalcyt [5]
[6]
albit skaleń-K muskowit
gibbsyt flogopit muskowit
[7]
gibbsyt albit flogopit
saponit-Na
muskowit [8]
Dla przyspieszonego tempa reakcji skalenia potasowego zmiany objętości matrycy skały wynikają z bilansu objętości pomiędzy ulegającymi rozkładowi minerałami pierwotnymi (fig. 3C), głównie klinochloru, albitu, skalenia potasowego i kalcytu, oraz minerałami krystalizującymi w systemie, ta- kimi jak: saponit-Na i saponit-Ca , muskowit, kaolinit, gibb- syt, paragonit i dolomit.
Całkowity rozkład skalenia potasowego następuje już po około 200 latach modelowanego etapu, co widoczne jest na wykresie jako poziomy przebieg linii zmian objętości lub ilości minerału (fig. 3D). Zauważalne są także i inne różni- ce w stosunku do scenariusza bazowego, polegające między innymi na:
– szybszym zakończeniu rozkładu klinochloru (czas oko- ło 80 lat);
– kontynuowanym rozkładzie flogopitu i paragonitu oraz trwającej krystalizacji saponitu-Na i muskowitu.
Z punktu widzenia wielkości mineralnej sekwestracji pierwiastka węgla podczas modelowanego składowania CO2 (fig. 2C) znaczenie mają zjawiska, które nie były obserwo- wane w scenariuszu bazowym. Są to rekrystalizacja kalcytu, powiązana z degradacją dolomitu, widoczne na figurze 3C i 3D, jako wygięcia odpowiednich krzywych w kierunku osi X, w czasie około 1500 lat oraz w końcowym okresie mode- lowanego etapu (około 6000 lat).
W przypadku scenariusza uwzględniającego przyspie- szone tempo reakcji klinochloru zmiany objętości matrycy skały są efektem bilansu pomiędzy ulegającymi rozkładowi minerałami pierwotnymi, głównie flogopitu, albitu, klino- chloru, kalcytu i skalenia potasowego, oraz minerałami kry- Fig. 2. Zróżnicowanie zmian porowatości próbki i wielkości mineralnej sekwestracji pierwiastka węgla na etapie składowania
co2, w zależności od stałej szybkości reakcji [mol/cm2/s] dla skalenia potasowego i klinochloru
Differences in sample porosity and mineral carbon sequestration at the CO2 storage stage, depending on the reaction rate constants [mol/cm2/s] for K-feldspar and clinochlore
K -skaleń
0 +1000 +2000 +3000 +4000 +5000 +6000 +7000 +8000 +9000+1e4 .238
.24 .242 .244 .246 .248 .25 .252
Porowatość
10-13
bazowa
10-15
Porosity base
B
0 +1000 +2000 +3000 +4000 +5000 +6000 +7000 +8000 +9000 +1e4 0
1 2 3 4 5 6 7 8 9
Pierwiastek C w skalezmiana ilości – [mol]
bazowa
10-14 10-16
base
C
pierwiastek C w skalezmiana ilości – [mol] ]lom [ atled kco nbarinorc–
0 +1000 +2000 +3000 +4000 +5000 +6000 +7000 +8000 +9000+1e4 0
1 2 3 4 5 6 7 8 9
bazowa 10-13
10-15
base
D klinochlor /clinochlore
Czas [lata]
Time [yr]
Pierwiastek C w skalezmiana ilości – [mol]
/ K-feldspar
Porowatość Porosity
Czas [lata]
Czas [lata]
Czas [lata]
Time [yr] Time [yr]
Time [yr]
Carbon in rock – delta [mol] Carbon in rock – delta [mol]
Fig. 3. Zmiany objętości i ilości faz mineralnych podczas realizacji scenariusza bazowego oraz scenariuszy, w których stałe szybkości reakcji skalenia potasowego i klinochloru podwyższono o jeden rząd wielkości
Changes in the volume and quantity of mineral phases during the baseline scenario, and the scenarios in which the K-feldspar and clinochlore reaction rate constants were increased by one order of magnitude bazowa / baseK-Skaleń / K-feldsparklinochlor /clinochlore
+.1 +1 +10 +100 +1000 +1e4
–1e3 –500 0 500 1000 1500
Czas [lata]
dolomite muscovite
k-feldspar
albit/albite klinochlor 14A/clinochlore 14A - kalcyt/ calcite
flogopit saponit-Ca/saponite-Ca gibbsite
saponite-Na
paragonit/ paragonite
Minerały icśoli anaimz– [cm ]3 Some minerals – delta [cm ] 3
saponit-Na
muskowit
gibbsyt
dolomit kaolinit/ kaolinite
skaleń potasowy phlogopite
Time [yr]
A
+.1 +1 +10 +100 +1000 +1e4
–8 –6 –4 –2 0 2 4 6 8 10 12
dolomit/dolomite muscovite
- phlogopite kaolinite
gibbsite
saponite-Na saponit-Na gibbsyt
muskowit
paragonit/ paragonite saponit-Ca/saponite-Ca
kalcyt/ calcite klinochlor 14A clinochlore 14A
skaleń potasowy k-feldspar
flogopit
albit/albite kaolinit
B
+.1 +1 +10 +100 +1000 +1e4
–1e3 –500 0 500 1000 1500
dolomit/dolomite
muscovite
- flogopit/phlogopite gibbsyt/gibbsite
saponite-Na
paragonit/paragonite
klinochlor 14A/clinochlore 14A albit/albite kalcyt/ calcite
kaolinit/ kaolinite
saponit-Na
muskowit
skaleń potasowy/ k-feldspar saponit-Ca/saponite-Ca
C
+.1 +1 +10 +100 +1000 +1e4
–8 –6 –4 –2 0 2 4 6 8 10
muscovite
- saponite-Na
dolomit/dolomite saponit-Ca/saponite-Ca
gibbsyt/gibbsite
kaolinit kaolinite
klinochlor 14A clinochlore 14A
kalcyt/ calcite
flogopit/phlogopite skaleń potasowy/ k-feldspar
saponit-Na
muskowit
albit/albite paragonit/paragonite
D
+.1 +1 +10 +100 +1000 +1e4
–1500 –1e3 –500 0 500 1000 1500
muscovite
- fengit/phengite
saponite-Na
albit/albite klinochlor 14A
clinochlore 14A kalcyt/ calcite
flogopit/phlogopite skaleń potasowy k-feldspar gibbsyt/gibbsite
Kaolinit kaolinite
dolomit/dolomite
saponit-Na
muskowit saponit-Ca/saponite-Ca
E
+.1 +1 +10 +100 +1000 +1e4
–8 –6 –4 –2 0 2 4 6 8 10
muscovite
-
saponite-Na
dolomit/dolomite
saponit-Na muskowit
fengit/phengite Kaolinit/kaolinite
gibbsyt/gibbsite saponit-Ca/saponite-Ca
klinochlor 14A clinochlore 14A
flogopit/phlogopite albit/albite kalcyt/ calcite
skaleń potasowy k-feldspar
F
Minerały icśoli anaimz– [mol ] Some minerals – delta [mol ] some minerals – delta [cm]
Minerały icśoli anaimz– [cm ]3 Some minerals – delta [cm ] some minerals – delta [cm]
3Minerały icśoli anaimz– [cm ]3 Some minerals – delta [cm ]
some minerals – delta [cm]
3 Minerały icśoli anaimz– [mol ] Some minerals – delta [mol ]
some minerals – delta [cm]
Minerały icśoli anaimz– [mol ] Some minerals – delta [mol ]
Czas [lata]
Time [yr]
Czas [lata]
Time [yr]
Czas [lata]
Time [yr]
Czas [lata]
Time [yr]
Czas [lata]
Time [yr]
stalizującymi w systemie, takimi jak: saponit-Na, muskowit, fengit, dolomit, saponit-Ca i kaolinit (fig. 3E).
Szybki rozkład klinochloru (fig. 3F) zostaje wyhamowa- ny już po około jednym roku, mimo że jego zasób w matry- cy skały został wyczerpany tylko w połowie. Dalej zmiany objętości lub ilości tego minerału są już nieznaczne. Zauwa- żalne są także i inne różnice w stosunku do scenariusza ba- zowego, polegające między innymi na:
– rozkładzie znacznych ilości flogopitu i krystalizacji muskowitu;
– krystalizacji znacznych ilości fengitu:
KAlMg3Si3O10(OH)2 + 4H2O + NaAlSi3O8 + 2H2S + 2CO2 +
+ 2KAlSi3O8 + 3KAlMgSi4O10(OH)2 + Al3+ + 2SO4 2– + 2CH4 +Na+ – rozkładzie kaolinitu, który mimo początkowego wzrostu
objętości, ulega degradacji i w rezultacie, od czasu około 1000 lat od początku etapu, jego objętość w matrycy skalnej jest niższa od początkowej (linia na wykresie wę- druje poniżej wartości zerowej);
– znikomej roli gibbsytu, który w niewielkich ilościach krystalizuje a następnie rozkłada się w początkowym okresie modelowanego etapu.
W mineralnej sekwestracji pierwiastka węgla (fig. 2D) znaczenie mają zjawiska rekrystalizacji kalcytu, powiązanej z degradacją dolomitu (fig. 3E i 3F). Są one jednak rozło- żone w czasie nieco inaczej te, które opisano wcześniej dla przypadku skalenia potasowego.
poDSumoWAnie
Przedstawiono przykład analizy wrażliwości rezultatów modelowania systemu woda–skała–gaz dla celów sekwestra- cji gazów kwaśnych na zmiany stałej kinetyki reakcji wybra- nych minerałów, w przykładowym poziomie wodonośnym.
Pośród najczęściej uwzględnianych parametrów uwzględ- nianych w modelowaniu hydrogeochemicznym dróg reakcji w układach woda–skała gaz te właśnie wielkości obarczone są wysoką niepewnością i obejmują szeroki zakres wartości.
Z rezultatami stanowiącymi poziom odniesienia (scena- riusz bazowy) porównywano wyniki uzyskane dla pozostałych dwóch scenariuszy, odpowiadających uwzględnieniu podwyż- szonych lub pomniejszonych o jeden rząd wielkości stałych szybkości reakcji dla skalenia potasowego i klinochloru.
O ile w przypadku obydwu minerałów w początkowym okresie modelowania (do około 1000 dla skalenia potaso- wego i 1800 lat dla klinochloru) tendencja zmian jest dla wszystkich scenariuszy podobna – gwałtowny wzrost poro- watości, po którym następuje jej intensywny spadek, a na- stępnie szybki wzrost, to w miarę upływu modelowanego czasu zaznaczają się wyraźne różnice. Końcowe porowato- ści obliczone dla wszystkich scenariuszy, dla skalenia pota- sowego, są do siebie zbliżone – różnice pomiędzy nimi nie
przekraczają 0,1 punktu procentowego (p.p.), podczas gdy dla klinochloru sięgają 0,5 p.p.
Zaobserwowane różnice wynikają ze zróżnicowania przebiegów reakcji, wpływających na ilości wody znajdu- jącej się w systemie oraz regulujących bilanse objętości mi- nerałów krystalizujących lub rozkładających się w trakcie każdego z modelowanych scenariuszy.
Z punktu widzenia wielkości mineralnej sekwestracji pierwiastka węgla podczas modelowanego składowania CO2 znaczenie mają zjawiska, które nie były obserwowane w scenariuszu bazowym: rekrystalizacja kalcytu powiązana z degradacją dolomitu.
W celu uzyskania wglądu w naturę różnic pomiędzy efektami uzyskiwanymi w kolejnych scenariuszach zidenty- fikowano główne reakcje w analizowanych systemach.
Dla przyspieszonego tempa reakcji skalenia potasowego zmiany objętości matrycy skały wynikają z bilansu objętości pomiędzy ulegającymi rozkładowi minerałami pierwotnymi, głównie klinochloru, albitu, skalenia potasowego i kalcytu, oraz minerałami krystalizującymi w systemie, takimi jak: sa- ponit-Na i saponit-Ca, muskowit, kaolinit, gibbsyt, parago- nit i dolomit. W stosunku do scenariusza bazowego, główne różnice polegają między innymi na:
– szybszym zakończeniu rozkładu klinochloru;
– kontynuowanym rozkładzie flogopitu i paragonitu oraz trwającej krystalizacji saponitu-Na i muskowitu.
W scenariuszu przyspieszonego tempa reakcji dla klino- chloru zaobserwowano zmiany wynikające głównie z roz- kładu flogopitu, albitu, klinochloru, kalcytu i skalenia pota- sowego, oraz krystalizacji: saponitu-Na, muskowitu, fengitu, dolomitu, saponitu-Ca i kaolinitu. Szybki rozkład klinochloru, był wyhamowany mimo, że jego zasób w matrycy nie zostaje całkowicie wyczerpany. Dalsze różnice w stosunku do scena- riusza bazowego, polegają między innymi na:
– rozkładzie znacznych ilości flogopitu i fengitu i krysta- lizacji muskowitu;
– początkowym wzroście objętości kaolinitu, po którym następuje jego degradacja;
– małej roli gibbsytu, który w niewielkich ilościach kry- stalizuje a następnie się rozkłada.
Różnice w opisywanych efektach petroteksturalnych i fizykochemicznych wydają się z pozoru niewielkie. Nale- ży jednak wziąć pod uwagę, że wynikają one z modyfikacji pojedynczego tylko parametru – stałej szybkości reakcji, dla wybranych faz mineralnych. W rzeczywistości przy wystę- powaniu wielu parametrów opisujących modelowany sys- tem, efekty z nimi związane mogą ze sobą interferować.
Przeprowadzona analiza wrażliwości pozwoliła na okreś- lenie wpływu zróżnicowania wielkości parametrów kinety- ki reakcji na uzyskiwane rezultaty modelowania interakcji w analizowanych systemach hydrogeochemicznych.
Przedstawione rezultaty uzyskano w ramach projektu sfinansowanego ze środków Narodowego Centrum Nauki przyznanych na podstawie decyzji numer DEC-2012/05/B/
ST10/00416.
flogopit albit
skaleń-K fengit [9]
literAturA
BETHKE C.M., 2008 – Geochemical and biogeochemical reaction modeling. Cambridge Univ. Press., Cambridge.
BETHKE C.M., YEAKEL S., 2013 – The Geochemist’s Work- bench Release 9.0. Reactive Transport Modeling Guide. Aque- ous Solutions, Illinois.
EKBERG C., 1999 – Sensitivity analysis and simulation uncer- tainties in predictive geochemical modelling. Freiberg Online Geoscience, 2: 1–109.
HOLLOWAY S., 2005 – Underground sequestration of carbon di- oxide e a viable greenhouse gas mitigation option. Energy, 30:
2318–2333.
LASAGA A.C., 1984 – Chemical kinetics of water–rock interac- tions. J. Geoph. Res., 89: 4009–4025.
LABUS K., SUCHODOLSKA K., 2015 – Rezultaty badań ekspe- rymentalnych i hydrogeochemicznego modelowania efektów interakcji między gazami kwaśnymi (CO2 i H2S) a środowi- skiem geologicznym wybranych skał GZW. Prz. Geol., 63, 10/2: 887–832.
LABUS K., TARKOWSKI R., WDOWIN M., 2015 – Modeling gas–rock–water interactions in carbon dioxide storage capac- ity assessment: a case study of Jurassic sandstones in Poland.
Int. J. Env. Sci. Tech., 12, 8: 2493–2502.
LABUS K., KASZA P., TUREK M., DYDO P., JAKóBIK-KOLON A., WILK K., LEŚNIAK G., 2016 – Laboratory experiments on forming the chemical composition of flowback from hydraulic fracturing with energized fluid. Procedia Earth Planet. Sci., 17: 582–585.
LUBAŚ J., SZOTT W., 2010 – 15-year experience of acid gas stor- age in the natural gas structure of Borzęcin – Poland. Nafta- Gaz, 66, 5: 333–338.
PALANDRI J.L., KHARAKA Y.K., 2004 – A compilation of rate parameters of water-mineral interaction kinetics for applica- tion to geochemical modeling. W: US Geological Survey Open File Report 2004 –1068:1–64.
SINAL M.L., LANCASTER G., 1987 – Liquid CO2 Fracturing: Ad- vantages And Limitations. J. Can. Petrol. Techn., 26, 5: 26–30.
XU T., APPS J.A., PRUESS K., YAMAMOTO H., 2007 – Numeri- cal modeling of injection and mineral trapping of CO2 with H2S and SO2 in a sandstone formation. Chem. Geol., 242: 319–346.
SummAry In the case of using computer modelling to simulate complex processes, such as interactions in hydrogeochemi- cal systems, the sensitivity analysis proves to be a suitable tool for assessing the impact of varying parameter sizes on the obtained results.
This paper presents an example of sensitivity analysis for changes in reaction kinetics parameters in relation to the results obtained in hydrogeochemical modelling of water–
rock–gas systems, for the purpose of sequestering acidic ga- ses (mixtures of carbon dioxide and hydrogen sulphide) in water–rock–gas systems.
The calculations were made to take into account the in- fluence of differentiation of reaction rate for K-feldspar and clinochlore on the reactions of decomposition and crystal- lization of mineral phases during the modelled storage of acidic gases, and petrotextural effects (changes in porosity),
as well as the amount of mineral sequestration of carbon in the rock formation. These baseline scenario results were compared to the results obtained for the other two scenarios, corresponding to the constant reaction rates increased or de- creased by one order of magnitude.
The differences result from the diversity of reactions, af- fecting the amount of water in the system and the volume balances of the minerals crystallizing or decomposing during each of the modelled scenarios.
Differences in the described petrotextural and physico- chemical effects seem to be small. However, it should be taken into account that they result from the modification of a single parameter only – the rate of reaction for selected mineral phases. In fact, when there are many parameters to describe the modelled system, the associated effects may in- terfere with each other.
