MODELOVÁNÍ TRANSPORTNˇE-CHEMICKÝCH PROCES˚U

(1)

MODELOV ´ AN´I

TRANSPORTNˇ E-CHEMICK ´ YCH PROCES˚ U

Dizertaˇ cn´ı pr´ ace

Studijn´ı program: P3901 – Aplikované vˇedy v inˇzenýrstv´ı Studijn´ı obor: 3901V025 – Pˇr´ırodovˇedné inˇzenýrstv´ı Autor práce: Ing. Lukáˇs Zedek

Vedouc´ı pr´ace: doc. Ing. Jan ˇSembera, Ph.D.

(2)

MODELING REACTIVE TRANSPORT PROCESSES

Dissertation

Study programme: P3901 – Applied Sciences in Engineering Study branch: 3901V025 – Science Engineering

Author: Ing. Luk´aˇs Zedek

Supervisor: doc. Ing. Jan ˇSembera, Ph.D.

(3)

Prohl´ aˇ sen´ı

Byl jsem seznámen s t´ım, ˇze na mou dizertaˇcn´ı práci se plnˇe vzta- huje zákon ˇc. 121/2000 Sb., o právu autorském, zejména § 60 – ˇskoln´ı d´ılo.

Beru na vˇedom´ı, ˇze Technická univerzita v Liberci (TUL) neza- sahuje do mých autorských práv uˇzit´ım mé dizertaˇcn´ı práce pro vnitˇrn´ı potˇrebu TUL.

Uˇziji-li dizertaˇcn´ı pr´aci nebo poskytnu-li licenci k jej´ımu vyuˇzit´ı, jsem si vˇedom povinnosti informovat o t´eto skuteˇcnosti TUL;

v tomto pˇr´ıpadˇe má TUL právo ode mne poˇzadovat úhradu náklad˚u, které vynaloˇzila na vytvoˇren´ı d´ıla, aˇz do jejich skuteˇcné výˇse.

Dizertaˇcn´ı práci jsem vypracoval samostatnˇe s pouˇzit´ım uvedené literatury a na základˇe konzultac´ı s vedouc´ım mé dizertaˇcn´ı práce a konzultantem.

Datum:

Podpis:

(4)

Abstrakt

Dizertaˇcn´ı práce, kterou drˇz´ıte v rukou je vˇenována problematice poˇc´ıtaˇcové simulace vybraných pˇr´ırodn´ıch jev˚u a chemických reakc´ı, které ovliˇnuj´ı ˇs´ıˇren´ı zneˇciˇstˇen´ı podzemn´ı vodou. V práci popiso- vané poˇc´ıtaˇcové modely vyuˇz´ıvaj´ı inovativn´ı aplikace v jiných oblastech zavedených postup˚u a modelovac´ıch nástroj˚u. Velká pozornost je v dizertaci vˇenována matematické formulaci popisu problému reakˇcn´ıho transportu kontaminace. V rámci práce je navrˇzen a byl provˇeˇren origináln´ı popis problému reakˇcn´ıho transportu soustavou obyˇcejných diferenciáln´ıch rovnic. Navrˇzený popis diferenciáln´ımi rovnicemi umoˇzˇnuje zahrnout do modelu i rovnováˇzné chemické reakce, které jsou obvykle popisovány rovnicemi algebraickými.

D´ıky zamˇeˇren´ı práce na matematický popis simulovaných úloh je moˇzné navrˇzené postupy analogicky aplikovat na modely dalˇs´ıch pˇr´ırodn´ıch a jiných jev˚u jako jsou napˇr´ıklad ˇs´ıˇren´ı tepla, mechanismus ztenˇcován´ı ozonové vrstvy, zneˇciˇstˇen´ı vzduchu a mnoho dalˇs´ıch.

Abstract

The book you hold in your hand describes computer simulation of selected natural processes and reactive-transport problem which influence groundwater contaminant transport. Models described in this book describe inovative application of approaches and simulation tools which are widely spread and often used in different simulation fields. This dissertation is strongly focused on mathematical formulation of reactive-transport problem description. In one part of the book we proposed and described tests of an original ordinary differential formulation of reactive transport problem. The proposed formulation enables to involve equlibrium reactions which are obviously described by algebraic equations.

Thanks to a focus of dissertation on mathematical formulation of reactive transport problem, proposed approaches can be used ana- logicaly for simulation of similarly described natural and other processes such as heat transfer, ozone depletion chemical mechanism, air polution model and many others.

(5)

Podˇ ekov´ an´ı

Na tomto m´ıstˇe bych rád podˇekoval lidem, kteˇr´ı mˇe nˇekterým z n´ıˇze uvedených zp˚usob˚u podpoˇrili pˇri pˇr´ıpravˇe této dizertaˇcn´ı práce.

Dˇekuji svému ˇskoliteli Doc. Ing. Janu ˇSemberovi, Ph.D. za odborné rady a konzultace, které mi v pr˚ubˇehu studia ochotnˇe poskytoval.

Dˇekuji za finanˇcn´ı podporou prostˇrednictv´ım zamˇestn´an´ı v r´amci projekt˚u

”Vlivy poˇcas´ı a zmˇen klimatu na ˇzivotn´ı prostˇred´ı a zemn´ı stavby-KLIPRO“ ˇc. 100113281 (program Ziel 3—C´ıl 3, vedouc´ı projektu Doc. Ing. Jan ˇSembera, Ph.D.), TA ˇCR TA01021331 (vedouc´ı projektu Doc. Ing. Jiˇrina Královcová, Ph.D.), grantu GA ˇCR 102/08/H081 (vedouc´ı projektu Prof. Ing. Aleˇs Richter, CSc.) a d´ıky úˇcasti ve Studentské grantové soutˇeˇzi (2014) na MTI FM TUL.

Dále dˇekuji mému otci za trpˇelivost, d˚uvˇeru, podporu a v nepo- sledn´ı ˇradˇe také za pravidelnˇe pˇripravované nedˇeln´ı obˇedy.

(6)

Seznam obr´ azk˚ u

2.1 Rozdˇelen´ı l´atek z pohledu modelu popisu chemick´ych reakc´ı. . . 23

2.2 Schematické znázornˇen´ı principu výmˇeny látek mezi mobiln´ımi a imobiln´ımi póry . . . 38

3.1 Geometrie modelu kolonového experimentu ovlivnˇeného chemickými reakcemi. A = C₂Cl₄, B = M nO⁻₄, C = C₂Cl₃H, D = H⁺, E = OH⁻, F = Cl⁻. . . 41

3.2 Rozpadov´a ˇrada pro ilustraci konstrukce matice rozpad˚u. . . 43

3.3 Geometrie modelu pro simulaci rozpad˚u. . . 44

3.4 Vliv zapnut´ı radioaktivn´ıch rozpad˚u. . . 46

3.5 Schema pro v´ypoˇcet koeficient˚u. . . 50

3.6 V´yvoj koncentrac´ı izotop˚u. . . 52

3.7 Vliv volby délky ˇcasového kroku na výsledky simulace radioaktivn´ıho rozpadu. . . 53

3.8 Diskretizace geometrie popisovan´e oblasti. . . 54

3.9 Hypotetické úloˇziˇstˇe (ˇcervenˇe) a kontaminovaná oblast (modˇre). . . . 55

3.10 Srovn´an´ı pˇredpovˇezen´ych mrak˚u kontaminace. . . 56

3.11 Srovn´an´ı v´ysledk˚u simulace rozpadu Cesia na Baryum. . . 56

3.12 Délky trván´ı výpoˇct˚u v závilosti na poˇctu pouˇzitých procesor˚u. . . . 58

3.13 Zrychlen´ı v´ypoˇct˚u vlivem paralelizace. . . 58

3.14 Omezenou rozpustnost´ı ovlivnˇen´e, izotermami popsan´e sorpce. . . 60

3.15 Hraniˇcn´ı pˇr´ımka pro simulaci sorpce. . . 61

3.16 Body do interpolaˇcn´ı tabulky v otoˇcen´e souˇradn´e soustavˇe. . . 63

3.17 Porovn´an´ı v´ysledk˚u simulac´ı rozpouˇstˇen´ı kalcitu. . . 71

3.18 Srovnán´ı výsledk˚u z Octave (modrá) a z programu The Geochemist’s Workbench (ˇcervená). . . 72

3.19 Porovn´an´ı v´ysledk˚u simulac´ı s pouˇzit´ım ODE formulace modelu. . . . 75

3.20 S´ıt’ pro ˇreˇsen´ı modelu reakˇcn´ıho transportu. . . 77

3.21 Shoda výsledk˚u simulace transport˚u v r˚uzných smˇerech. N´ızká rychlost. 79 3.22 Porovnán´ı výsledk˚u simulac´ı reakˇcn´ıho transportu. . . 79

3.23 Nastaven´ı okrajov´ych podm´ınek s vyuˇzit´ım aktu´aln´ıch hodnot koncentrac´ı. . . 80

3.24 Vliv zjemnˇen´ı s´ıtˇe na porovnatelnost v´ysledk˚u simulac´ı. . . 81

3.25 Zmˇena sledovan´ych bod˚u pro model transportu v horizont´aln´ım smˇeru. 81 3.26 Vliv zjemnˇen´ı s´ıtˇe ve smˇeru osy y. . . 82

(7)

3.27 Vliv zjemnˇen´ı s´ıtˇe ve smˇeru osy x. . . 82 3.28 Vliv zjemnˇen´ı s´ıtˇe ve smˇeru proudˇen´ı na rozloˇzen´ı pH pod´el osy x. . 83

(8)

Seznam tabulek

3.1 Popis kinetick´ych reakc´ı. . . 47

3.2 V´ypoˇcetn´ı ˇcasy region´aln´ıho modelu, stoln´ı PC . . . 57

3.3 V´ypoˇcetn´ı ˇcasy region´aln´ıho modelu, cluster Hydra . . . 57

3.4 Pod´ıly sekvenciáln´ıch ˇcást´ı výpoˇct˚u regionáln´ıho modelu (fractioning). 58 3.5 Parametry modelu sorpce. . . 66

3.6 Casy simulac´ı transportu ovlivnˇˇ en´eho sorpc´ı. . . 66

3.7 Zmˇena znaˇcen´ı . . . 70

3.8 V´ypoˇcetn´ı ˇcasy model˚u reakˇcn´ıho transportu. . . 84

(9)

Obsah

Znaˇcen´ı a seznam zkratek . . . 11

Uvod´ 13 1 C´ıle dizertaˇcn´ı práce 15 2 Souˇcasný stav problematiky 16 2.1 Kategorizace chemických reakc´ı . . . 18

2.2 Známé postupy matematické formulace problému . . . 20

2.2.1 Terminologie a znaˇcen´ı . . . 20

2.2.2 Rovnice popisuj´ıc´ı probl´em . . . 24

2.2.3 Pˇr´ıstupy k numerick´emu ˇreˇsen´ı . . . 30

2.2.4 Reakˇcnˇe orientovan´e modelov´an´ı . . . 35

2.3 Existuj´ıc´ı software . . . 37

2.3.1 Flow123D a Semchem . . . 37

3 Pˇr´ıklady vlastn´ıch model˚u reakˇcn´ıho transportu 40 3.1 Vyuˇzit´ı metody ˇstˇepen´ı oper´atoru k propojen´ı program˚u Flow123D a Semchem . . . 41

3.2 Radioaktivn´ı rozpady . . . 43

3.2.1 Zp˚usob simulace radioaktivn´ıch rozpad˚u . . . 44

3.2.2 Vliv zapnut´ı radioaktivn´ıch rozpad˚u . . . 45

3.2.3 Kinetick´e reakce prvn´ıho ˇr´adu . . . 47

3.2.4 Rozˇsiˇritelnost postupu maticov´eho n´asoben´ı . . . 47

3.2.5 Zohlednˇen´ı rozpad˚u pˇredch´azej´ıc´ıch izotop˚u . . . 49

3.2.6 Casov´ˇ a n´aroˇcnost v´ypoˇct˚u, vliv paralelizace . . . 54

3.3 Rovnov´aˇzn´e sorpce . . . 59

3.3.1 Formulace modelu . . . 59

3.3.2 Reˇsen´ı modelu . . . 62ˇ 3.3.3 Implementace a testov´an´ı modelu . . . 65

3.4 Simulace paraleln´ıch chemick´ych reakc´ı . . . 68

3.4.1 Numerick´e ˇreˇsiˇce . . . 68

3.4.2 DAE model rozpouˇstˇen´ı kalcitu . . . 69

3.4.3 ODE model rozpouˇstˇen´ı kalcitu . . . 72

3.5 PDE Model reakˇcn´ıho transportu . . . 76

4 Pˇr´ınos pr´ace 85

(10)

A DAE soustava, reakce, DASPK 90

B DAE soustava, reakce, IDA 92

C DAE soustava, reakce, IDA-Jacobiho matice 94

D ODE soustava, reakce, CVODE 97

E PDE soustava, transport, CVODE 99

F Publikaˇcn´ı a jin´e v´ysledky studia 105

(11)

Znaˇ cen´ı a seznam zkratek

α_i modifikovaná rovnováˇzná konstanta i-té reakce

a^k_ij, b^k_ij, k ∈ {x, y, p} stechiometrick´e koeficienty komplexaˇcn´ıch, sorpˇcn´ıch a sr´aˇzec´ıch reakc´ı, [a^k_ij] = [b^k_ij] = 1

c_a koncentrace specie v kapaln´e f´azi, [c_a] = M · L⁻³ = kg · m⁻³

c_j koncentrace j-t´e specie rozpuˇstˇen´e ve vodˇe, [c_j] = M · L⁻³ = kg · m⁻³

C_j koncentrace j-t´e rozpuˇstˇen´e komponenty, [C_j] = M · L⁻³ = kg · m⁻³

−

→c^R vektor koncentrac´ı rozpuˇstˇených a sorbovaných speci´ı v ro- tované souˇradné soustavˇe

c_s koncentrace sorbovan´e specie, [c_s] = M · L⁻³ = kg · m⁻³ nebo [c_s] = M · M⁻¹ = kg · kg⁻¹

D tenzor hydrodynamick´e disperze, [D] = L²· T⁻¹ = m²· s⁻¹ DAE Differential-Algebraic Equation

f_izo(. . .) funkˇcn´ı pˇredpis adsorpce izotermou γ_i aktivitn´ı koeficient i-t´e specie, [γ_i] = [_c¹

i] = m³· kg⁻¹

GWB The Geochemist’s Workbench

k kinetick´a konstanta reakce

K_i, L, M rovnov´aˇzn´e konstanty reakc´ı

L(. . .) advekˇcnˇe difuzn´ı/disperzn´ı oper´ator L(c_i) = ∇ · (vc_i) − ∇ · (θD∇c_i)

m_i molalita i-t´e specie, [m_i] = N · L⁻³ = mol · kg⁻¹ M celkov´y poˇcet speci´ı M = M_a+ M_s+ M_p

M_a poˇcet vodou transportovan´ych speci´ı M_a= M_x+ N_a M_p poˇcet speci´ı ve formˇe sraˇzenin

M_s v reˇserˇsn´ı ˇc´asti, poˇcet netransportovan´ych, imobiln´ıch speci´ı M_s= M_y+ N_s

M_s v ˇcásti o vlastn´ıch modelech, molárn´ı hmotnost sorbované specie

M_x poˇcet speci´ı ve formˇe komplex˚u M_y poˇcet sorbovan´ych speci´ı

n por´ozita, [n] = 1

N celkový poˇcet komponentn´ıch speci´ı N = N_a+ N_s, plat´ı pro reˇserˇsn´ı ˇcást práce

N_a poˇcet rozpuˇstˇen´ych komponentn´ıch speci´ı N_s poˇcet komponentn´ıch speci´ı sorbentu ODE Ordinary Differential Equation

OS Operator Splitting

p tlaková výˇska, [p] = L = m p_j koncentrace j-té sraˇzené specie,

[p_j] = M · L⁻³ = kg · m⁻³

(12)

Znaˇ cen´ı a seznam zkratek

P_j koncentrace j-t´e komponenty ve sraˇzenin´ach, [Pj] = M · L⁻³ = kg · m⁻³

q darcyovsk´a rychlost, [q] = L · T⁻¹= m · s⁻¹

r_k reakˇcn´ı rychlost, [r_k] = M · L⁻³· T⁻¹ = kg · m⁻³· s⁻¹ ρw hustota vody (transportn´ıho m´edia), [ρw] = kg · m⁻³ s_j koncentrace j-t´e komponentn´ı specie sorbentu,

[s_j] = M · L⁻³ = kg · m⁻³

Sj koncentrace j-t´e sorbovan´e komponenty, [S_j] = M · L⁻³ = kg · m⁻³

S_w saturace, [S_w] = 1

SIA Sequential Iterative Approach SNIA Sequential Non-iterative Approach

τ_i poloˇcas rozpadu i-t´eho izotopu, [τ_i] = T = rok θ objemov´a vlhkost, θ = Sw· n, [θ] = 1

T-H-M-C Thermal-Hydrodynamical-Mechanical-Chemical x_j koncentrace j-t´eho rozpuˇstˇen´eho komplexu,

[xj] = M · L⁻³ = kg · m⁻³

y_j koncentrace j-t´e sorbovan´e specie, [y_j] = M · L⁻³= kg · m⁻³

Tj celkov´a koncentrace j-t´e komponenty Tj = Cj+ Sj+ Pj, [T_j] = M · L⁻³ = kg · m⁻³

v pr˚umˇern´a rychlost vody v p´orech, v = q/n, [v] = L·T⁻¹= m·s⁻¹ V, Velm objem, objem elementu, [V ] = [Velm] = L³ = m³

V_s objem pevn´e f´aze, [V_s] = L³ = m³

W_j celkov´a koncentrace j-t´e komponenty sorbentu, [Wj] = M · L⁻³= kg · m⁻³

X_i aktivita i-t´e specie, [X_i] = 1 ζ^(r), x, y, w, z rozsahy reakc´ı, [ζ^(r)] = N = mol

(13)

Uvod ´

Tato dizertaˇcn´ı práce pˇripravená v rámci absolvován´ı studijn´ıho oboru

”Pˇr´ıro- dovˇedné inˇzenýrstv´ı“ leˇz´ı na pomez´ı výzkumných smˇer˚u

”Modelov´an´ı pˇr´ırodn´ıch a technick´ych proces˚u“ a

”Aplikace numerických metod a vývoj nástroj˚u pro mode- lován´ı pˇr´ırodn´ıch proces˚u“. ˇSirˇs´ı tematický okruh práce nese oznaˇcen´ı

”Numerické modelován´ı sdruˇzených T-H-M-C proces˚u v horninovém prostˇred´ı“. Sepsaná práce je vˇenována návrhu a popisu efektivn´ıch postup˚u poˇc´ıtaˇcové simulace vybraných geochemických dˇej˚u a je rozˇclenˇena do ˇctyˇrech n´ıˇze popsaných kapitol.

Prvn´ı kapitola kr´atce, v nˇekolika bodech, pˇredstavuje c´ıle dizertaˇcn´ı pr´ace.

Druhá kapitola pˇrináˇs´ı souhrn základn´ıch znalost´ı a postup˚u vyuˇz´ıvaných pro simulaci reakˇcn´ıho transportu kontaminace podzemn´ı vodou. Kapitola obsahuje in- formace o obvyklém dˇelen´ı chemických reakc´ı v podzemn´ı vodˇe do nˇekolika skupin podle rychlosti reakc´ı a fáze reaktant˚u. Dále je v kapitole uvedeno nˇekolik mate- matických formulac´ı vyuˇz´ıvaných k poˇc´ıtaˇcovému modelován´ı geochemických reakc´ı r˚uzných druh˚u. Druhou kapitolu uzav´ırá krátké pˇredstaven´ı simulaˇcn´ıho nástroje Flow123D do nˇehoˇz byly v rámci praktické ˇcásti dizertaˇcn´ı práce pˇridány souˇcásti pro simulaci vybraných geochemických, popˇr. geofyzikáln´ıch interakc´ı.

Tˇeˇziˇstˇem práce je kapitola ˇc´ıslo 3, která popisuje vlastn´ı navrˇzené postupy a metody modelován´ı jev˚u, které ovlivˇnuj´ı, zpravidla zpomaluj´ı, transport zneˇciˇstˇen´ı podzemn´ı vodou.

V jedné z podkapitol naleznete origináln´ı matematickou formulaci popisu radioaktivn´ıch rozpad˚u a lineárn´ıch reakc´ı, která byla v nˇekolika kroc´ıch zpˇresˇnována a peˇclivˇe provˇeˇrována.

Dalˇs´ı podkapitola pˇredstavuje postup simulace rovnováˇzné sorpce, který umoˇzˇnuje pˇredej´ıt ˇcasovˇe nároˇcnému iteraˇcn´ımu ˇreˇsen´ı nelineárn´ıch rovnic d´ıky vyuˇzit´ı pˇredpoˇc´ı- tané interpolaˇcn´ı tabulky.

Následuj´ıc´ı podkapitola ukazuje aplikaci ˇreˇsiˇc˚u z bal´ıku SundialsTB na simulaci rozpouˇstˇen´ı kalcitu podzemn´ı vodou bez kontaktu s atmosférou. V této podkapitole je k matematickému popisu problému nejprve vyuˇzita soustava sestávaj´ıc´ı z algebraických spolu s obyˇcejnými diferenciáln´ımi rovnicemi (standardn´ı popis che- mických reakc´ı). Následnˇe je popis pˇreformulován do podoby soustavy obyˇcejných diferenciáln´ıch rovnic, které nesou informaci o kinetických i rovnováˇzných reakc´ıch (inovace). Výsledky model˚u rozpouˇstˇen´ı kalcitu byly srovnávány s výstupy z provˇe- ˇrených simulaˇcn´ıch program˚u.

Tˇret´ı kapitolu uzav´ırá ˇcást vˇenovaná rozˇs´ıˇren´ı modelu rozpouˇstˇen´ı kalcitu o ad- vekˇcn´ı transport látek podzemn´ı vodou.

Ctvrt´ˇ a kapitola shrnuje pˇr´ınosy odevzdan´e dizertaˇcn´ı pr´ace.

(14)

V pˇr´ılohách práce naleznete skripty pouˇzité v programu Octave k simulaci rozpouˇstˇen´ı kalcitu a k modelován´ı rozpouˇstˇen´ım kalcitu ovlivnˇeného transportu.

Kromˇe toho obsahuj´ı pˇr´ılohy seznam publikac´ı a v´yˇcet dalˇs´ıch aktivit v pr˚ubˇehu studia.

(15)

1. C´ıle dizertaˇ cn´ı pr´ ace

Vzhledem k velké ˇs´ıˇri problematiky poˇc´ıtaˇcové simulace reakˇcn´ıho transportu jsem se v dizertaˇcn´ı práci zamˇeˇril pˇredevˇs´ım na matematickou formulaci a na zp˚usoby ˇreˇsen´ı soustav popisuj´ıc´ıch daný problém. C´ıle pˇredkládané práce je moˇzné shrnout do nˇekolika následuj´ıc´ıch bod˚u:

• Zmapován´ı matematického pozad´ı poˇc´ıtaˇcového modelován´ı reakˇcn´ıho transportu kontaminant˚u podzemn´ı vodou.

• Návrh inovativn´ıho matematického popisu problému reakˇcn´ıho transportu.

• Pˇr´ıprava vlastn´ıch, efektivn´ıch algoritm˚u pro vybran´e pˇr´ıpady reakˇcn´ıho transportu kontaminace.

• Implementace navrˇzen´ych postup˚u a jejich pouˇzit´ı k simulaci zvolen´ych pˇr´ırod- n´ıch jev˚u.

• Porovnán´ı dosaˇzených výsledk˚u simulac´ı s výsledky z´ıskanými vyuˇzit´ım pro- vˇeˇrených, komerˇcn´ıch program˚u.

V této práci je uvaˇzován deterministický pˇr´ıstup/model (pro popis viz kapitola 2). Motivac´ı pro tuto volbu je zamýˇslené pouˇzit´ı výsledného modelovac´ıho nástroje pro stanoven´ı rizik spojených s ukládán´ım radioaktivn´ıho odpadu a pro predikci

´

uˇcinnosti sanaˇcn´ıch z´asah˚u na oblastech region´aln´ıho rozsahu.

Praktickou ˇcást´ı pro ovˇeˇren´ı teoretických závˇer˚u dizertaˇcn´ı práce je bud’to implementace a otestován´ı vybraných postup˚u do reakˇcn´ıho modulu simulaˇcn´ıho softwaru Flow123D nebo jejich zprovoznˇen´ı pomoc´ı jiného programového nástroje (Octave).

Uˇćelem rozˇs´ıˇren´ı programu Flow123D o simulaci vybraných geochemických interakc´ı bylo vytvoˇren´ı prostˇredku efektivn´ıho modelován´ı tˇechto jev˚u a nikoliv pˇr´ıprava plnohodnotného reakˇcnˇe-transportn´ıho modelovac´ıho nástroje jako je HYDROGE- OCHEM [1], The Geochemist’s Workbench (viz www.gwb.com) nebo PHT3D (viz http://www.pht3d.org/).

Vlastn´ı práce se soustˇredila na optimalizaci algoritm˚u a výpoˇcetn´ıch proces˚u pouˇzitých k simulaci reakˇcn´ıho transportu stejnˇe jako na inovaci matematického popisu problému.

(16)

2. Souˇ casn´ y stav problematiky

V souvislosti s aktuáln´ı potˇrebou likvidace ˇrady ekologických zátˇeˇz´ı vyvstává vˇzdy otázka vhodné, efektivn´ı volby sanaˇcn´ıho zásahu. K pˇredbˇeˇznému posouzen´ı

´

uˇcinnosti zásahu m˚uˇze pomoci dobˇre postavený poˇc´ıtaˇcový model ˇs´ıˇren´ı kontaminace. Poˇc´ıtaˇcové modely ˇs´ıˇren´ı kontaminace se vzájemnˇe liˇs´ı poˇctem do modelu zahrnutých pˇr´ırodn´ıch jev˚u, které mohou transport kontaminant˚u ovlivnit. Poˇcet a druh jev˚u zahrnovaných do poˇc´ıtaˇcových model˚u transportu kontaminace podzemn´ı vodou se s ˇcasem vyv´ıjel a vedl k modifikaci chápán´ı pojmu

”modelován´ı reakˇcn´ıho transportu“. Pˇrehled tohoto vývoje je moˇzné nalézt v [2]. V souˇcasnosti produkované modely reakˇcn´ıho transportu ˇcasto zohledˇnuj´ı tzv. T-H-M-C (Thermal, Hydrodynamic, Mechanical, Chemical) procesy doplnˇené nav´ıc o vliv ˇzivých orga- nism˚u na pr˚ubˇeh chemických reakc´ı (viz [3]). Do simulace proudˇen´ı podzemn´ı vody jde zahrnout vliv hustoty podzemn´ı vody (viz [4, 5]). Nˇekteré modelovac´ı nástroje (viz [6, 7]) zohledˇnuj´ı dokonce moˇzné, reakcemi vyvolané zmˇeny fyzikáln´ıch vlastnost´ı (porozita, propustnost) horninové matrice.

Tato dizertaˇcn´ı práce se soustˇredila na nalezen´ı a vyhodnocen´ı efektivnosti postup˚u pro simulaci transportu kontaminant˚u, ovlivnˇeného pouze chemickými reakcemi. V tomto smyslu je dále pouˇz´ıván i pojem

”reakˇcn´ı transport“.

Navzdory tomu, ˇze problematika reakˇcn´ıho transportu kontaminant˚u je v hle- dáˇcku modeláˇr˚u uˇz nˇekolik desetilet´ı, pˇredstavuje efektivn´ı algoritmizace problému pro modeláˇre neustávaj´ıc´ı výzvu. Kromˇe moˇznosti kontaminace podzemn´ı vody che- mikáliemi je tento zájem v dneˇsn´ı dobˇe motivován napˇr´ıklad snahou o nalezen´ı op- timáln´ıch poˇzadavk˚u na izolaci úloˇziˇst’ radioaktivn´ıho odpadu (viz [6, 8, 9]). Potˇreba zefektivˇnován´ı modelovac´ıch postup˚u trvá i pˇres neustávaj´ıc´ı nár˚ust výpoˇcetn´ıho výkonu poˇc´ıtaˇc˚u. Zvýˇsený výpoˇcetn´ı výkon umoˇzˇnuje pouˇzit´ı model˚u pro pˇredpo- vˇed’ ˇs´ıˇren´ı kontaminace v rámci vˇetˇs´ıch geografických celk˚u, které jsou v modelu reprezentovány s´ıtˇemi s vyˇsˇs´ım poˇctem element˚u.

Modely transportu chemikáli´ı kapalinou mohou být podle pˇr´ıstupu k modelo- ván´ı rozdˇeleny do dvou skupin. Jednou skupinou jsou modely s deterministickým základem a druhou skupinou jsou modely se základem v teorii pravdˇepodobnosti (stochastické). V deterministických modelech je proudˇen´ı podzemn´ı vody obvykle popsáno parciáln´ımi diferenciáln´ımi rovnicemi v prostoru a v ˇcase (Euler˚uv pˇr´ıstup).

Pˇr´ıtomnost kontaminantu v podzemn´ı vodˇe je interpretována prostˇrednictv´ım koncentrac´ı tˇechto látek nebo jejich analogi´ı vyjadˇruj´ıc´ıch mnoˇzstv´ı látky vztaˇzené na jednotkové mnoˇzstv´ı rozpouˇstˇedla. Neurˇcitost je do deterministického modelu zane- sena expertn´ımi odhady fyzikáln´ıch vlastnost´ı horniny a volbou parametr˚u mate- matického popisu pˇr´ırodn´ıch jev˚u. Deterministický pˇr´ıstup je vhodný tehdy, pokud

(17)

má simulované proudˇen´ı pˇreváˇznˇe charakter konvekce nebo také tehdy, kdyˇz je po- pisovaná oblast, na které k proudˇen´ı docház´ı, rozsáhlá.

Ve stochastických modelech je kapalina zpravidla vn´ımána jako soubor ˇcástic (Lagrange˚uv pˇr´ıstup) a trajektorie ˇcástic kontaminantu m˚uˇze být vyhodnocována na základˇe pravdˇepodobnosti výskytu (resp. prvn´ıho pr˚uchodu) ˇcástice v bodˇe. T´ımto zp˚usobem jsou konstruovány modely zaloˇzené na tzv.

”Theory of First Passage Processes“(viz [10]). Modely s pravdˇepodobnostn´ım základem se hod´ı pro simulaci transportn´ıch jev˚u s pˇrevládaj´ıc´ım difuzn´ım charakterem. Stochastické modely se naopak nehod´ı pro simulace transportn´ıch jev˚u na rozsáhlých oblastech, protoˇze popis kapaliny jako souboru ˇcástic by kladl neúnosné poˇzadavky na velikost pamˇeti poˇc´ıtaˇce.

Zefektivnˇen´ı postup˚u simulace reakˇcn´ıho transportu m˚uˇze být realizováno násle- duj´ıc´ımi zp˚usoby:

• Vyuˇzit´ım dostupných programových nástroj˚u pro ˇreˇsen´ı zkoumaného problému (ˇreˇsiˇce DAE, soustav nelineárn´ıch rovnic a symbolické ˇreˇsiˇce, viz kapitola 3.4).

• Zmˇenou matematiké formulace problému (viz kapitoly 3.2, 3.2.5, 3.3, 3.4.3 a 3.5). Tradiˇcn´ı pˇr´ıstupy jsou zachyceny v [11]. Automatizaci vybraných tra- diˇcn´ıch pˇr´ıstup˚u ukazuje [12].

• Pouˇzit´ım vhodnˇejˇs´ıch numerick´ych metod.

• Volbou efektivnˇejˇs´ı implementace v programu (PETSc, paralelizace).

Jednou z nejˇcastˇeji citovaných prac´ı o popisované problematice je [13] kde J. Ru- bin v roce 1983 sumarizoval základy modelován´ı reakˇcn´ıho transportu kontaminant˚u.

Také Rubinovou prac´ı, mimo jiné, se inspirovali modeláˇri na pensylvánské státn´ı uni- verzitˇe kde byl pod veden´ım G. T. Yeha a G. H. Iskry od konce osmdesátých let do roku 1995 vyvinut programový bal´ık KEMOD (Kinetic and Equilibrium MODel) vycházej´ıc´ı z freewarového produktu MINTEQ. Text v reˇserˇsn´ı ˇcásti tez´ı vycház´ı z velké ˇcásti z Rubinova ˇclánku [13] a dále pak ze zdroj˚u [11, 14, 12, 15] na jejichˇz vzniku se pod´ılel výˇse zm´ınˇený Gour-Tsyh Yeh a které obsahuj´ı popis rozvinut´ı a vyuˇzit´ı z dˇr´ıvˇejˇska známých technik pro ˇradu komplexn´ıch model˚u.

(18)

2.1 Kategorizace chemick´ ych reakc´ı

Z pohledu tvorby poˇc´ıtaˇcových model˚u reakˇcnˇe-transportn´ıch proces˚u je vhodné rozˇclenit chemické reakce do nˇekolika skupin (viz [13]), které se vyznaˇcuj´ı moˇznos- t´ı pouˇzit´ı stejného nebo alespoˇn podobného matematického popisu. Rozdˇelen´ı je následuj´ıc´ı:

1) Dostateˇcnˇe rychlé, vratné chemické reakce A) Homogenn´ı reakce

B) Heterogenn´ı reakce I) Povrchov´e reakce.

II) Klasick´e reakce.

2) Nedostateˇcnˇe rychlé a/nebo nevratné chemické reakce A) Homogenn´ı reakce

B) Heterogenn´ı reakce I) Povrchov´e reakce.

II) Klasick´e reakce.

Dostateˇcnost rychlosti chemických reakc´ı se vyhodnocuje na základˇe porovnán´ı rychlosti reakce s rychlost´ı transportu. Mezi tˇr´ıdami dostateˇcnˇe a nedostateˇcnˇe rychlých chemických reakc´ı nejsou pevnˇe vytyˇcené hranice a zaˇrazen´ı jedné, konkrétn´ı reakce do skupiny se m˚uˇze pˇr´ıpad od pˇr´ıpadu liˇsit. V rámci této dizertaˇcn´ı práce odpov´ıdá oznaˇcen´ı dostateˇcnˇe rychlé reakce vratným, rovnováˇzným reakc´ım. Ne- dostateˇcnˇe rychlými (tzv. kinetickými) reakcemi jsou nejˇcastˇeji m´ınˇeny reakce ne- vratné.

Za homogenn´ı je reakce oznaˇcována tehdy, pokud v jej´ım pr˚ubˇehu reaguj´ı jenom specie jedné fáze (tj. de facto jenom specie rozpuˇstˇené). Naopak jako heterogenn´ı je reakce oznaˇcována, kdyˇz spolu reaguj´ı specie r˚uzné fáze (tj. roztok s horninou ˇci se sraˇzeninou nebo kapalná, pˇr´ıpadnˇe pevná fáze s plynem). Praktickým pˇr´ıkladem heterogenn´ı reakce je rozpouˇstˇen´ı kalcitu popsané v kapitole 3.4.

Rozd´ıl mezi povrchovými a klasickými heterogenn´ımi reakcemi je nejlépe patrný z uveden´ı pˇr´ıkladu typických zástupc˚u tˇechto skupin. Reprezentativn´ımi zástupci povrchových reakc´ı jsou sorpce a desorpce (tj. reakce kapaliny s porézn´ım médiem, horninou, kterou kapalina proud´ı). Naproti tomu ke klasickým heterogenn´ım reakc´ım patˇr´ı sráˇzec´ı a rozpouˇstˇec´ı reakce. Rozpouˇstˇec´ı reakce m˚uˇze být schematicky popsána rovnic´ı

P1(s) + P2(aq) → P3(aq),

podle které reaguje specie P₁(s)-solid, pevná látka, se speci´ı P₂(aq)-aqueous, roz- puˇstˇená látka, za vzniku specie P3(aq)-rozpuˇstˇené látky.

Pr˚ubˇeh nedostateˇcnˇe rychlých reakc´ı je popsán obyˇcejnými diferenciáln´ımi rovnicemi (napˇr. (2.15)), které definuj´ı takzvanou reakˇcn´ı rychlost (r_i, [r_i] = mol · m⁻³·

(19)

s⁻¹). Reakˇcn´ı rychlost r_i popisuje rychlost zmˇeny koncentrace i-té specie. Stano- ven´ı reakˇcn´ı rychlosti má obvykle empirický charakter a vycház´ı z ˇcasových ˇrad namˇeˇrených dat o pozorovaných koncentrac´ıch.

Reakˇcn´ı sloˇzka model˚u, které zahrnuj´ı dostateˇcnˇe rychlé (rovnováˇzné) reakce, je popsána algebraickými rovnicemi. Tyto rovnice vyjadˇruj´ı pod´ıl mezi koncentracemi (resp. aktivitami) výchoz´ıch látek a produkt˚u. Zm´ınˇený pod´ıl z˚ustává po dosaˇzen´ı chemické rovnováhy systému zachován. Je jimi tedy definována tzv. rovnováˇzná konstanta reakce. Rovnováˇzná konstanta i-té reakce (popˇr. modifikovaná rovnováˇzná konstanta i-té reakce) se obvykle znaˇc´ı K_i (resp. α_i). Definice tˇechto konstant jsou uvedeny dále (viz (2.9) v kapitole 2.2.2).

(20)

2.2 Zn´ am´ e postupy matematick´ e formulace probl´ e- mu

2.2.1 Terminologie a znaˇ cen´ı

Pˇred vlastn´ım popisem matematických model˚u je tˇreba zavést dva pojmy. Je to pojem komponenta a pojem specie. Komponenty tvoˇr´ı mnoˇzinu takových chemických entit, jejichˇz unikátn´ımi kombinacemi je moˇzné popsat kaˇzdou specii, která figuruje v simulované reakci. Zároveˇn s t´ım nesm´ı být ˇzádná z komponent vyjádˇritelná kombinac´ı komponent ostatn´ıch.

Specie jsou produkty chemických reakc´ı mezi komponentami (resp. jsou vyjádˇritel- né jako kombinace/slouˇceniny komponent).

Následuj´ıc´ı znaˇcen´ı pro komponenty a specie je pˇrevzaté z [11]. N´ıˇze zavedené znaˇcen´ı bude pouˇzito v kapitolách 2.2.2 a 2.2.3 pro pˇredstaven´ı obvyklých zp˚usob˚u simulace základn´ıch typ˚u geochemických reakc´ı.

Koncentrace speci´ı:

• cjznaˇc´ı koncentraci j-té specie rozpuˇstˇené ve vodˇe. Jedná o takzvanou rozpuˇstˇe- nou komponentn´ı specii, která je tvoˇrena jen jednou komponentou a která nen´ı produktem ˇzádné z uvaˇzovaných chemických reakc´ı, [c_j] = kg · m⁻³.

• s_j je znaˇcka pro mnoˇzstv´ı j-té komponentn´ı specie sorbentu, tj. komponentn´ı specie v pevné horninové matrici. Komponentn´ı specie sorbentu m˚uˇze v chemic- kých reakc´ıch figurovat jen v roli výchoz´ı látky. Specie s_j je tvoˇrena pouze jednou komponentou, [s_j] = kg · m⁻³.

• x_j popisuje mnoˇzstv´ı j-tého mobiln´ıho komplexu (tj. specie) pˇr´ıtomného ve vodˇe. Komplex x_j je slouˇceninou nˇekolika rozpuˇstˇených komponentn´ıch speci´ı c_j, [x_j] = kg · m⁻³.

• y_jje mnoˇzstv´ı j-té sorbované specie vzniklé reakc´ı komponentn´ıch speci´ı c_jv ka- palné fázi s komponentn´ımi speciemi s_j v pevné fázi. K sorpˇcn´ım reakc´ım docház´ı na povrchu horniny v pórech, kterými protéká kapalina, [y_j] = kg·m⁻³.

• p_j je znaˇcka pro mnoˇzstv´ı j-té specie vzniklé sráˇzec´ı reakc´ı ve vodˇe pˇr´ıtomných rozpuˇstˇených komponentn´ıch speci´ı c_j nebo komplex˚u x_j, [p_j] = kg · m⁻³. Koncentrace komponent:

• C_j je celkové mnoˇzstv´ı j-té komponenty rozpuˇstˇené ve vodˇe bud’ ve formˇe rozpuˇstˇené komponentn´ı specie nebo jako souˇcást komplex˚u, [C_j] = kg · m⁻³.

• S_j znaˇc´ı mnoˇzstv´ı j-té komponenty, která která je sorbována z roztoku na povrch horninové matrice, [S_j] = kg · m⁻³.

• P_j je symbol pro mnoˇzstv´ı j-té komponenty, která je souˇcást´ı sraˇzenin, [Pj] = kg · m⁻³.

(21)

• T_j je celkové mnoˇzstv´ı j-té komponenty, která je v systému (resp. v elementu) pˇr´ıtomná v podobˇe rozpuˇstˇené komponentn´ı specie, jako souˇcást mobiln´ıho komplexu, sraˇzeniny nebo je sorbovaná na povrch horniny (viz (2.7)),

[T_j] = kg · m⁻³.

• Wj znaˇc´ı celkové analytické mnoˇzstv´ı j-té komponenty sorbentu v komplexech a v komponentn´ıch speci´ıch sorbent˚u (viz (2.8)), [W_j] = kg · m⁻³.

Fyzikáln´ı rozmˇery uvedených mnoˇzstv´ı ve vodˇe pˇr´ıtomných komponent a speci´ı jsou

[C_j] = [P_j] = [T_j] = [c_j] = [x_j] = [y_j] = [p_j] = M · L⁻³ = kg · m⁻³,

(tj. pomˇer hmotnosti j-té komponenty nebo specie k objemu kapaliny Vaq = V · θ, kde θ znaˇc´ı bezrozmˇernou objemovou vlhkost (moisture), která je k saturaci S_w ve vztahu θ = S_w· n, kde n znaˇc´ı porózitu horniny).

Rozmˇery veliˇcin popisuj´ıc´ıch koncentrace komponent a speci´ı sorbent˚u jsou dány pomˇerem hmotnosti sorbentu vztaˇzeného na objem pevné fáze (porézn´ı horniny, V_s= V · (1 − n), kde n je opˇet znaˇcka pro porózitu).

[W_j] = [S_j] = [s_j] = M · L⁻³ = kg · m⁻³. Znaˇcky pro poˇcty speci´ı r˚uzn´ych druh˚u:

Dále je zavedeno znaˇcen´ı pro poˇcty r˚uzných druh˚u speci´ı pˇr´ıtomných v sledo- vaném systému. Symboly N_i, i ∈ {a, s} jsou znaˇceny poˇcty komponentn´ıch speci´ı a symboly M_i, i ∈ {a, p, s, x, y} oznaˇcuj´ı poˇcty speci´ı vzniklých chemickými reakcemi.

• N_a znaˇc´ı poˇcet ve vodˇe pˇr´ıtomn´ych (a-aqueous), rozpuˇstˇen´ych komponentn´ıch speci´ı.

• N_sje symbolem pro poˇcet imobiln´ıch komponentn´ıch speci´ı sorbentu horninov´e matrice (s-solid).

• N oznaˇcuje celkový poˇcet komponen´ıch speci´ı, tedy N = N_a+ N_s. Chemickými reakcemi N_a ve vodˇe rozpuˇstˇených speci´ı vzniká

• M_x speci´ı v podobˇe komplex˚u a

• M_p speci´ı ve formˇe sraˇzenin.

Poˇcet vodou transportovaných, mobiln´ıch, speci´ı Ma je rovný souˇctu Mx mobiln´ıch komplex˚u (produkt˚u reakc´ı) a N_arozpuˇstˇených komponentn´ıch speci´ı (výcho- z´ıch látek reakc´ı)

• M_a = M_x+ N_a.

Chemické reakce mohou prob´ıhat také mezi speciemi v kapalné a speciemi v pevné fázi. Nˇekteré z Ma rozpuˇstˇených speci´ı a z Nssorbent˚u reaguj´ı do podoby My sorbo- vaných speci´ı. Pokud je do modelu zahrnuta sorpce, pak bude poˇcet speci´ı v pevné fázi rovný souˇctu N_s speci´ı sorbentu a M_y povrchovˇe komplexovaných nebo iontovou výmˇennou sorbovaných speci´ı. Vyjádˇreno vzorcem je

(22)

• M_s = N_s+ M_y.

Prob´ıhá-li sorpce mechanismem iontové výmˇeny, mohou se sorpˇcn´ıch proces˚u

´

uˇcastnit jak nˇekteré z Na ve vodˇe pˇr´ıtomných rozpuˇstˇených komponentn´ıch speci´ı, tak i vybrané z M_x komplex˚u. V takovém pˇr´ıpadˇe je poˇcet sorbovaných speci´ı

• M_y = M_x^k+ N_a^k,

kde horn´ı index u X_j^k, X_j ∈ {N_a, M_x} znaˇc´ı skuteˇcnost, ˇze se iontové výmˇeny vˇzdy nutnˇe neúˇcastn´ı vˇsechny rozpuˇstˇené ani vˇsechny komplexované specie, tedy ˇze M_x^k <= Mx, N_a^k <= Na.

Uvaˇzovaných M_s speci´ı v pevné fázi je imobiln´ıch.

Z výˇse uvedeného vyplývá, ˇze celkový poˇcet speci´ı M je souˇctem poˇct˚u speci´ı v kapalné fázi, poˇctu speci´ı v pevné fázi a M_p speci´ı ve formˇe sraˇzenin. Tedy,

• M = M_a+ M_s+ M_p.

Transportovány jsou rozpuˇstˇené komponentn´ı specie a komplexy, zat´ımco sraˇzené specie, sorbované komponentn´ı specie a komplexy transportovány nejsou.

Pˇr´ıklad pro ilustraci zaveden´eho znaˇcen´ı:

Pro vytvoˇren´ı pˇredstavy o rozloˇzen´ı speci´ı a komponent v podzemn´ı vodou za- plnˇené horninˇe poslouˇz´ı obrázek 2.1. Uvedený obrázek ilustruje ˇrez porézn´ı horninou.

Objem pór˚u, vyplnˇený vodou, je na obrázku vyznaˇcen b´ılou barvou. Pˇr´ıkladem rozpuˇstˇených komponentn´ıch speci´ı jsou specie c_j, j ∈ {1, . . . , N_a} sloˇzené vˇzdy pouze z jedné komponenty C_j.

Reakcemi speci´ı c_j, tj. kombinacemi vybran´ych komponent z C_i, i ∈ {1, . . . , N_a}, mohou vznikat rozpuˇstˇen´e, mobiln´ı komplexy x_j, j ∈ {1, . . . , M_x}.

C´ˇast obrázku s ˇsedým pozad´ım reprezentuje imobiln´ı produkty p_j, j ∈ {1, . . . , M_p} sráˇzec´ıch reakc´ı mezi rozpuˇstˇenými komponentn´ımi speciemi nebo/a mobiln´ımi komplexy.

Samotná hornina je vyznaˇcena zelenou barvou. Ve spodn´ım z rámeˇck˚u na ze- leném pozad´ı je zachycen pˇr´ıklad komponentn´ı specie sorbentu s_j, tvoˇrené jen komponentou S_j.

Cervenou barvou jsou vyznaˇˇ ceny oblasti se speciemi vzniklými sorpc´ı povrchovu komplexac´ı nebo iontovou výmˇenou. Pˇr´ıkladem takové specie je y_j, která se skládá z vybraných komponent C_j, j ∈ {1, . . . , N_a} a S_j, j ∈ {1, . . . , N_s}.

Kromˇe pˇr´ıklad˚u rozm´ısttˇen´ı r˚uzných druh˚u komponent a speci´ı v systému jsou na obrázku zachyceny také znaˇcky pouˇzité pro poˇcty uvaˇzovaných speci´ı (N_i, i ∈ {a, s}

a M_j, j ∈ {a, p, s, x, y}).

Doplˇnuj´ıc´ı znaˇcky pro popis transportu a reakc´ı:

Kromˇe symbol˚u pro koncentrace jsou v textu kapitoly 2.2 pouˇz´ıvány dalˇs´ı znaˇcky, které vystupuj´ı v chemických rovnic´ıch nebo jako parametry matematických model˚u.

• v je symbol pro pr˚umˇernou rychlost vody v p´orech, [v] = L · T⁻¹ = m · s⁻¹.

(23)

Obrázek 2.1: Rozdˇelen´ı látek z pohledu modelu popisu chemických reakc´ı.

• D je tenzor hydrodynamické disperze (materiálová konstanta), [D] = L² · T⁻¹ = m²· s⁻¹.

• L(. . .) znaˇc´ı advekˇcnˇe difuzn´ı/disperzn´ı oper´ator, L(ci) = ∇ · (vci) − ∇ · (θD∇ci).

• a^k_ij, b^k_ij k ∈ {x, y, p} (k znaˇc´ı typ reakce bud’ komplexotvorná, sorpˇcn´ı nebo sráˇzec´ı) je znaˇcka pro stechiometrický koeficient, který pˇr´ısluˇs´ı j-té specii v rovnici popisuj´ıc´ı i-tou chemickou reakci. Koeficienty a^k_ij jsou bezrozmˇerné, tedy

[a^k_ij] = [b^k_ij] = 1.

• K_i je oznaˇcen´ı rovnováˇzné konstanty i-té reakce.

• X_j znaˇc´ı aktivitu j-t´e specie, [X_j] = 1.

• γ_i oznaˇcuje aktivitn´ı koeficient pˇr´ısluˇsej´ıc´ı i-t´e specii, [γ_i] = [_c¹

i] = m³· kg⁻¹.

• m_i je znaˇcka pro molalitu, [m_i] = N · M⁻¹ = mol · kg⁻¹.

(24)

2.2.2 Rovnice popisuj´ıc´ı probl´ em

V literatuˇre je obvykle reakˇcn´ı transport popsán advekˇcnˇe difuzn´ımi parciáln´ımi diferenciáln´ımi rovnicemi (viz [16]) doplnˇenými o zdrojový ˇclen. Rovnice reakˇcn´ıho transportu má potom tvar

∂(θc_i)

∂t + ∇ · (vci) − ∇ · (θD∇ci) + S(~c, t) = 0. (2.1) Tato rovnice popisuje skuteˇcnost, ˇze koncentrace i-té specie se m˚uˇze v uvaˇzova- ném modelu zmˇenit vlivem advekce, difuze nebo vlivem jiného pˇr´ırodn´ıho jevu, který je zohlednˇen prostˇrednictv´ım zdrojového ˇclenu S(~c, t). S výjimkou ˇclenu ∂(θc_i)/∂t tak vˇsechny sˇc´ıtance na levé stranˇe pˇredstavuj´ı

”pˇr´ıspˇevky“ ke zmˇenˇe koncentrace i-té specie. V dalˇs´ım textu budou m´ıt rovnice transportu obvykle formálnˇe m´ırnˇe odliˇsný tvar neˇz (2.1). Mohou se liˇsit pouˇzitým znaˇcen´ım zavedeným v kapitole 2.2.1 a pˇr´ıtomnost´ı reakˇcn´ıho ˇclenu.

N´ıˇze uvedené rovnice jsou pˇrevzaty z [11] a v rámci nˇekolika následuj´ıc´ıch podkapitol popisuj´ı matematické modely transportu ovlivnˇené rovnováˇznými a kinetickými reakcemi.

V prvn´ıch dvou rovnic´ıch (2.2),(2.3) vystupuj´ı jako promˇenné koncentrace komponent. Uvedené rovnice popisuj´ı zmˇeny distribuce koncentrac´ı komponent v sys- tému v ˇcase a v prostoru vlivem transportu a v d˚usledku zmˇen saturace. Rovnice (2.2) je popisem vývoje koncentrac´ı komponent v kapalné fázi (ve vodˇe vyplˇnuj´ıc´ı póry), zat´ımco rovnice (2.3) je bilanˇcn´ı rovnice pro imobiln´ı komponenty v pevné fázi (v horninˇe s póry zaplnˇenými vodou):

θ∂T_j

∂t + ∂θ

∂t(Sj+ Pj) = −∇(v · Ci) + ∇[θD∇Ci], j = 1, . . . , Na (2.2) θ∂W_j

∂t + ∂θ

∂tW_j = 0, j = 1, . . . , N_s (2.3) Koncentrace komponent nejsou podle definice pojmu komponenta ovlivnˇeny pr˚u- bˇehem chemických reakc´ı a proto neobsahuje (2.2) zdrojový ˇclen zohledˇnuj´ıc´ı vliv chemických reakc´ı. Popis transportu parciáln´ımi diferenciáln´ımi rovnicemi s koncentracemi komponent na m´ıstˇe promˇenných je typický pro modely reakˇcn´ıho transportu s dostateˇcnˇe rychlými reakcemi.

Rovnice (2.4) aˇz (2.8) zachycuj´ı vztah mezi koncentracemi komponent (velk´ymi p´ısmeny) a koncentrace speci´ı (mal´ymi p´ısmeny):

C_j = c_j +

Mx

X

i=1

a^x_ijx_i, j = 1, . . . , N_a (2.4)

S_j =

My

X

i=1

a^y_ijy_i, j = 1, . . . , N_a (2.5)

P_j =

Mp

X

i=1

a^p_ijp_i, j = 1, . . . , N_a (2.6)

(25)

Tj = cj +

Mx

X

i=1

a^x_ijxi+

My

X

i=1

a^y_ijyi+

Mp

X

i=1

a^p_ijpi, j = 1, . . . , Na (2.7) Zaps´ana jin´ym zp˚usobem dostane (2.7) tvar Tj = Cj + Sj + Pj.

Pro celkov´a mnoˇzstv´ı komponent sorbentu plat´ı rovnice W_j = s_j +

My

X

i=1

a^y_ijy_i, j = 1, . . . , N_s. (2.8) Symboly a^k_ij, k ∈ {x, y, p} v rovnic´ıch (2.4) aˇz (2.8) jsou znaˇcky pro stechiomet- rické koeficienty stoj´ıc´ı v chemických rovnic´ıch u pˇr´ısluˇsných speci´ı.

Rovnice popisuj´ıc´ı rovnov´ahu

Rovnováˇzné reakce jsou vratné chemické reakce, které vedou k témˇeˇr okamˇzitému dosaˇzen´ı a k udrˇzen´ı takzvané dynamické rovnováhy mezi speciemi na levé a na pravé stranˇe chemické rovnice, která popisuje danou reakci. Reakce mezi speciemi na levé stranˇe rovnice a reakce mezi speciemi na pravé stranˇe rovnice prob´ıhaj´ı po dosaˇzen´ı rovnováhy stejnou rychlost´ı a koncentrace speci´ı v pozorovaném m´ıstˇe se tak vlivem reakce nemˇen´ı.

Uvaˇzujme rovnováˇznou reakci mezi speciemi A a B na jedné stranˇe a speciemi C a D na stranˇe druhé. V systému jsou dále také v rakci nefiguruj´ıc´ı specie E, F . Prob´ıhaj´ıc´ı reakce necht’ je popsána rovnic´ı

A + B C + D.

Koncentrace speci´ı A, B, C, D z˚ustávaj´ı d´ıky dynamické rovnováze za- chovány. Koncentrace speci´ı E, F se také nemˇen´ı, protoˇze ani E, ani F se uvaˇzované rovnováˇzné reakce neúˇcastn´ı. K popisu takovéto rov- nováhy slouˇz´ı tzv. rovnováˇzná konstanta reakce (viz [11]).

Pˇr´ıklad 1.:

Kaˇzdá z reakc´ı má svou rovnováˇznou konstantu (K_i), která je funkc´ı teploty a kterou je moˇzné analyticky odvodit na základˇe termodynamických zákon˚u [17].

Odvozená rovnováˇzná konstanta je pod´ılem souˇcin˚u mocnin koncentrac´ı produkt˚u a reaktant˚u rovnováˇzné reakce. V literatuˇre (viz [18]) je moˇzné nalézt vyjádˇren´ı rovnováˇzné konstanty reakce jako pod´ılu ne koncentrac´ı, ale takzvaných aktivit.

Aktivity se zavád´ı z toho d˚uvodu, ˇze specie se pˇri reakc´ıch zpravidla chovaj´ı takovým zp˚usobem, jako by jich reagovalo jiné, obvykle menˇs´ı, mnoˇzstv´ı neˇz je jich v systému pˇr´ıtomno. Takové chován´ı si je moˇzné pro pˇredstavu vysvˇetlit vlivem fyzikáln´ıch interakc´ı mezi ˇcásticemi chemikáli´ı v roztoku. Z uvedených d˚uvod˚u se nˇekdy pro aktivitu pouˇz´ıvá pojem efektivn´ı koncentrace. Aktivita je ovlivnˇena pˇr´ıtomnost´ı dalˇs´ıch speci´ı v roztoku a poˇc´ıtá se jako souˇcin aktivitn´ıho koeficientu a koncentrace rozpuˇstˇené specie

Xi = γi· ci, [Xi] = 1.

(26)

Aktivitn´ı koeficienty pro specie, které se úˇcastn´ı simulované reakce, je moˇzné vypoˇc´ıtat napˇr´ıklad na základˇe Debye-Hückelovy teorie elektrolyt˚u publikované roku 1923 v [19]. Alternativn´ı zp˚usob výpoˇctu aktivitn´ıch koeficient˚u je moˇzné nalézt v [20]. Fyzikáln´ı rozmˇer aktivitn´ıho koeficientu je [γ_i] = [_c¹

i].

Pokud je systém v dynamické chemické rovnováze, potom mezi aktivitami speci´ı v systému plat´ı vztah

K_i =

Na

Y

j=1

X_j^a^ij = konst., (2.9)

kde Ki je rovnováˇzná konstanta reakce a X_jâîk jsou aktivity vˇsech speci´ı v systému (vˇcetnˇe E, F ) umocnˇené na stechiometrické koeficienty a_ik z rovnice uvaˇzované reakce. Stechiometrické koeficienty z levé strany chemické rovnice rovnováˇzné reakce jsou podle konvence záporná celá ˇc´ısla. Koeficienty na pravé stranˇe rovnice jsou kladná celá ˇc´ısla. Z mocnˇen´ı aktivit na stechiometrický koeficient z rovnice reakce je patrné, ˇze nen´ı potˇreba násobit aktivitami speci´ı, které v reakci nevystupuj´ı. Stechi- ometrické koeficienty speci´ı, které v reakci nefiguruj´ı, jsou nulové. Formáln´ı zahrnut´ı nereaguj´ıc´ıch speci´ı E, F do modelu umoˇzˇnuje výrazné zjednoduˇsen´ı matematického popisu problému, protoˇze indexy aktivit speci´ı ve vzorci (2.9) mohou identifikovat vˇsech Na speci´ı v systému. Ve výˇse uvedeném pˇr´ıkladu jsou stechiometrické koeficienty −1 pro X_j ∈ {A, B}, 1 pro X_j ∈ {C, D} a 0 pro E, F .

Protoˇze se pro popis pˇr´ıtomnosti soeci´ı v urˇcitém m´ıstˇe pouˇz´ıvaj´ı ˇcastˇeji koncentrace neˇz aktivity, je vhodné pˇrepsat vztah (2.9) právˇe pro koncentrace:

K_i =

Na

Y

j=1

X_j^a^ij =

Na

Y

j=1

(γ_jc_j)^a^ij =

Na

Y

j=1

γ_j^a^ij ·

Na

Y

j=1

c^a_j^ij.

Souˇcin aktivitn´ıch koeficient˚u lze pˇrev´est na levou stranu rovnice, ˇc´ımˇz z´ısk´ame tzv.

modifikovanou rovnováˇznou konstantu α_i pˇr´ısluˇsej´ıc´ı i-té chemické reakci:

α_i = K_i·

Na

Y

j=1

γ_j^−a^ij = K_i Y

j∈R

γ_j^−a^ij Y

k∈P

γ_k^a^ik ,

kde R je mnoˇzina index˚u reaktant˚u a P je mnoˇzina index˚u produkt˚u.

Zahrnut´ı rychl´ych reakc´ı do rovnice reakˇcn´ıho transportu

Z rovnic popisuj´ıc´ıch rovnováˇzné reakce je patrné, ˇze tyto nedefinuj´ı, na rozd´ıl od rovnic pro popis kinetických reakc´ı, pˇr´ıspˇevek ke zmˇenˇe koncentrace sledované specie.

Z tohoto d˚uvodu nen´ı moˇzné pˇr´ımo zahrnout výˇse uvedené rovnice rovnováhy do zdrojového ˇclenu rovnice (2.1). Pˇritom je zˇrejmé, ˇze neˇz dosáhnou specie v systému rovnováhy, mus´ı se jejich koncentrace nˇejakou rychlost´ı mˇenit. Rychlost r takové zmˇeny m˚uˇze být vyjádˇrena následuj´ıc´ım zp˚usobem:

r = k · lnQ

K, (2.10)

(27)

kde K = lim_c_i_→c^eq

i , i∈{p,r}

Y

p∈P

c^−a_p ^p Y

r∈R

câ_r^r je rovnováˇzná konstanta reakce, Q = Y

p∈P

c^−a_p ^p Y

r∈R

câ_r^r je takzvaný iontový souˇcin a k pˇredstavuje empiricky zjiˇstˇenou kinetickou konstantu.

Vyjádˇren´ı rychlosti pr˚ubˇehu rovnováˇzných reakc´ı má význam z hlediska pˇred- stavy moˇznosti matematického popisu reakˇcn´ıho transportu s rovnováˇznými a ki- netickými reakcemi jednou rovnic´ı. V modelech se popisu (2.10) ˇcasto nevyuˇz´ıvá a rychlost rovnováˇzných reakc´ı je pokládána za nekoneˇcnou (viz odvozen´ı z termo- dynamických zákon˚u v [12]). Popis (2.10) plat´ı pro stav ve kterém je systém bl´ızký rovnováze.

Jiné vyjádˇren´ı pˇr´ıspˇevku rovnováˇzné reakce ke zmˇenˇe koncentrac´ı speci´ı je moˇzné naj´ıt v kapitole 3.4.3 vˇenované simulaci rozpouˇstˇen´ı kalcitu a v kapitole 3.5 popisuj´ıc´ı model rozpouˇstˇen´ı kalcitu spolu se simulac´ı advekˇcn´ıho transportu.

Základn´ı druhy rovnováˇzných reakc´ı

Po struˇcném teoretickém úvodu do základ˚u popisu rovnováˇzných chemických reakc´ı pˇricház´ı v této podkapitole na ˇradu pˇredstaven´ı obvyklých matematických formulac´ı nejˇcastˇeji simulovaných druh˚u rovnováˇzných reakc´ı. Rovnováˇzné reakce tˇrech základn´ıch typ˚u (komplexace, sorpce, sráˇzen´ı) mohou být pro matematické modely popsány rovnicemi (2.11) aˇz (2.14). Pro jednoduchost je v této podkapitole pˇredpokládáno, ˇze j-tému (resp. i-tému) produktu rovnováˇzné reakce pˇr´ısluˇs´ı ste- chiometrický koeficient rovný 1.

Pro rovnov´aˇznˇe prob´ıhaj´ıc´ı komplexace ve vodˇe plat´ı

x_j = α^x_i

Na

Y

k=1

c^a_k^x^ik, j = 1 . . . M_x, (2.11)

kde x_j znaˇc´ı koncentraci vznikaj´ıc´ıho j-tého komplexu, α^x_i je modifikovaná rov- nováˇzná konstanta (viz [11]) i-té komplexaˇcn´ı reakce a câ

x ik

k jsou koncentrace komponentn´ıch speci´ı, které se komplexace úˇcastn´ı, umocnˇené na pˇr´ısluˇsné stechiometrické koeficienty.

Pro sorpce prob´ıhaj´ıc´ı povrchovou komplexac´ı plat´ı

y_i = α^y_i(

Na

Y

k=1

c^a

y ik

k )(

Ns

Y

k=1

s^b

y ik

k ), i = 1 . . . M_y, (2.12) kde y_i je koncentrace i-tého vznikaj´ıc´ıho komplexu, α^y_i je pˇr´ısluˇsná modifikovaná rov- nováˇzná konstanta (viz [11]) a násobeny jsou mocniny koncentrac´ı komponentn´ıch speci´ı v kapalné (câ

y ik

k ) a v pevn´e f´azi (s^b

y ik

k ).

Sorpce iontovu v´ymˇenou je obvykle popisov´ana funkc´ı

y_i = f_i(c₁, . . . , c_N_a), i = 1, . . . , M_s (2.13)

(28)

a pro sr´aˇzec´ı/rozpouˇstˇec´ı reakce plat´ı pˇredpis

1 = α^p_i

Na

Y

k=1

c^a

p ik

k , i = 1, . . . , M_p (2.14)

Rovnice (2.14) neobsahuje koncentrace pi, protoˇze je pˇredpokládána aktivita speci´ı v pevné fázi rovna 1. T´ım se popis sráˇzec´ıch reakc´ı liˇs´ı od ostatn´ıch heterogenn´ıch reakc´ı jako jsou napˇr´ıklad sorpce nebo iontová výmˇena na povrchu horniny.

Rovnice popisuj´ıc´ı kinetiku

Rada reakc´ı (kinetick´ˇ ych) prob´ıh´a v´yraznˇe pomaleji neˇz transportn´ı procesy.

Z d˚uvodu relativn´ı pomalosti kinetických reakc´ı, je tˇreba popsat jejich pr˚ubˇeh. K popisu pr˚ubˇehu kinetické reakce se vyuˇz´ıvá veliˇcina nazývaná rychlost reakce (anglicky

”reaction rate“).

Rychlost i-té reakce se znaˇc´ı r_i a popisuje reakc´ı vyvolanou ˇcasovou zmˇenu mnoˇzstv´ı specie (koncentrace, molality, molarity) vztaˇzeného na jednotkový stechi- ometrický koeficient. Pro jednoduchost je, pro produkty popisovaných reakc´ı, dále uvaˇzován stechiometrický koeficient 1. Matematicky je rychlost reakce popsána, na pˇr´ıkladu komplexace (vratné reakce) ve vodném roztoku, vzorcem

r_i^x = ∂xi

∂t = −k^b_i^xx_i+ k_i^f^x

Na

Y

k=1

c^a_k^ik^x, i = 1, . . . , K_x. (2.15)

Vratné reakce mohou prob´ıhat obˇema smˇery. Proto vystupuj´ı v rovnici (2.15) dvˇe kinetické konstanty k_iâ, a ∈ {b_x, f_x}. k_i^b^x je znaˇcka pro kinetickou konstantu reakce prob´ıhaj´ıc´ı ve smˇeru, ve kterém se komplex x_i rozpadá (pro nevratné reakce k_i^b^x = 0, v rovnici zprava doleva) a k^f_i^x je kinetická konstanta reakce prob´ıhaj´ıc´ı ve smˇeru ve kterém komplex x_i vzniká (v rovnici zleva doprava). Kinetickou konstantu reakce pro i-tou specii je moˇzné urˇcit z ˇcasové ˇrady koncentrac´ı dané specie empirickými postupy. Fyzikáln´ı rozmˇer kinetických konstant je v uvaˇzovaném pˇr´ıpadˇe

[k^j_i] = T⁻¹ = s⁻¹.

Pokud je mechanismus vzniku nˇekter´e ze speci´ı kinetick´y, tak je pro odpov´ıdaj´ıc´ı popis tohoto mechanismu v modelu tˇreba upravit dvˇe vˇeci.

• Za prvé je tˇreba nahradit rovnice popisu chemické rovnováhy takové, ve kterých daná specie vystupuje. Náhrada se provede bud’ rovnic´ı popisuj´ıc´ı transport (pro komplexaci ve vodˇe, viz (2.16)) nebo za rovnici vyjadˇruj´ıc´ı zákon za- chován´ı hmoty (v pˇr´ıpadˇe sorpce nebo sráˇzen´ı, viz (2.17), (2.20)).

• Za druhé je tˇreba vyjádˇrit rychlost reakce pomoc´ı výrazu ve kterém vystupuj´ı pouze koncentrace speci´ı a komponentn´ıch speci´ı (tj. výchoz´ıch látek a produkt˚u chemických reakc´ı). K vyjádˇren´ı rychlosti reakc´ı nen´ı moˇzné vyuˇz´ıt koncentrace komponent, protoˇze mnoˇzstv´ı komponent se reakc´ı nemˇen´ı.

(29)

V dalˇs´ım textu je pouˇzit advekˇcnˇe disperzn´ı diferenci´aln´ı oper´ator

L(c) = −v · ∇c + ∇ · (θD · ∇c) = −∇(v · c + θD · ∇c), kde v je pórová rychlost prou- dˇen´ı a D je tenzor hydrodynamické disperze (tj. tenzor materiálových konstant hor- ninové matrice).

N´ıˇze uvedené rovnice (2.16-2.20) jsou pˇrevzaty z [11]. Pro simulaci kinetického pr˚ubˇehu komplexaˇcn´ıch reakc´ı ve vodném prostˇred´ı je tˇreba nahradit popis rov- nováhy (2.11) rovnic´ı transportu doplnˇenou o zdrojový ˇclen

θ∂x_i

∂t = L(x_i) + θr^x_i, i = 1, . . . , M_x, (2.16) kde M_x je poˇcet komplexovan´ych speci´ı a r^x_i je rychlost komplexaˇcn´ı reakce (viz (2.15)). [r^x_i] = M · L⁻³· T⁻¹ = kg · m⁻³· s⁻¹.

Pro simulaci kinetického pr˚ubˇehu sorpˇcn´ı reakce je tˇreba nahradit popis pˇr´ısluˇsné rovnováhy bilanˇcn´ı rovnic´ı

θ∂yi

∂t +∂θ

∂ty_i = θr_i^y, i = 1, . . . , M_y, (2.17) kde r_i^y je rychlost kinetické sorpce i-té specie prob´ıhaj´ıc´ı mechanismem iontové výmˇeny a M_y je, analogicky k pˇredchoz´ımu pˇr´ıpadu, poˇcet speci´ı sorbovaných na povrch horniny modelovaným mechanismem.

Pro rychlosti sorpc´ı prob´ıhaj´ıc´ıch mechanismem iontov´e v´ymˇeny plat´ı

r^y_i = −k^b_i^yy_i+ k^f_i^y

Na

Y

k=1

c^a_k^ik^y

Ns

Y

k=1

s^b_k^ik^y, i = 1, . . . , M_y (2.18)

a pro rychlosti sorpc´ı prob´ıhaj´ıc´ıch prostˇrednictv´ım povrchov´e komplexace uvaˇzuje- me

r^y_i = −k_i^b^yyi+ k_i^f^yfi(c1, c2, . . . , cNs), i = 1, . . . , My. (2.19) Pˇri simulaci sráˇzec´ı reakce je, pro modelován´ı kinetického pr˚ubˇehu reakce, rovnice popisu rovnováˇzné reakce (2.14) nahrazena zákonem zachován´ı

θ∂p_i

∂t +∂θ

∂tp_i = θr_i^p, i = 1, . . . , M_p. (2.20) Rychlost r^p_i kinetických sráˇzec´ıch reakc´ı i-té specie je dána vztahem

r^p_i = −k_i^b^pp_i+ k_i^f^p

Na

Y

k=1

c^a_k^ik^p, i = 1, . . . , M_p. (2.21)

(30)

2.2.3 Pˇ r´ıstupy k numerick´ emu ˇ reˇ sen´ı

Výˇse zapsané vztahy (2.2)-(2.8) a (2.11)-(2.14) tvoˇr´ı soustavu popisuj´ıc´ıch hyd- rogeologický transport spolu s dostateˇcnˇe rychlými, rovnováˇznými komplexacemi, redoxn´ımi reakcemi, neutralizacemi, sorpcemi (podle druhu sorpce popsané popsané bud’ (2.12) nebo (2.13)) a sráˇzec´ımi/rozpouˇstˇec´ımi reakcemi. S vyuˇzit´ım informac´ı z kapitoly 2.2.2 m˚uˇze být simulován kinetický pr˚ubˇeh reakc´ı.

Soustavu lze nˇekolika zp˚usoby zjednoduˇsit. Podrobný popis technik zmenˇsen´ı di- menze uvedené soustavy je poskytnut v [13, 12]. Kl´ıˇcovou otázkou pro volbu mezi n´ıˇze uvedenými postupy ˇreˇsen´ı soustavy popisuj´ıc´ı reakˇcn´ı transport je výbˇer tak- zvaných primárn´ıch a sekundárn´ıch promˇenných. Pˇredstavované postupy totiˇz spoˇc´ı- vaj´ı, velmi zjednoduˇsenˇe ˇreˇceno, v r˚uzných rozdˇelen´ıch soustavy na dva subsystémy se dvˇema sadami promˇenných. Oba subsytémy mohou být ˇreˇseny, podle výbˇeru metody, v r˚uzné m´ıˇre samostatnˇe. Nejprve jsou v kaˇzdém ˇcasovém kroku vypoˇcteny primárn´ı promˇenné a sekundárn´ı promˇenné jsou urˇceny dodateˇcnˇe.

Popisované postupy se vzájemnˇe liˇs´ı ˇcasovou nároˇcnost´ı za jejich pomoci provádˇe- ných simulac´ı a druhy do model˚u zahrnutelných typ˚u chemických reakc´ı. Jednot- livým pˇr´ıstup˚um k numerickému ˇreˇsen´ı jsou dále ponechány p˚uvodn´ı anglické názvy pˇrevzaté z [11]. Pro ˇreˇsen´ı soustav nelineárn´ıch algebraických rovnic, které v obecném pˇr´ıpadˇe zanáˇs´ı do modelu jeho reakˇcn´ı sloˇzka, se pouˇz´ıvá Newton-Raphsonova nebo Picardova metoda (viz [21]).

Differential and Algebraic Equation Approach-DAE

Prvn´ı druh formulace matematického popisu reakˇcnˇe transportn´ıch proces˚u nese oznaˇcen´ı Diferential and Algebraic Equation Approach (DAE). Jak uˇz název nazna- ˇcuje, soustava popisuj´ıc´ı problém je v pˇr´ıpadˇe tohoto pˇr´ıstupu sloˇzena dohromady jak z diferenciáln´ıch tak algebraických, obvykle nelineárn´ıch rovnic. Pro ˇreˇsen´ı takové soustavy rovnic je nutné vyuˇz´ıt odpov´ıdaj´ıc´ı numerické metody nebo ˇreˇsiˇce (napˇr.

Newton-Raphsonova metoda, viz kapitola 2.3.1 nebo ˇreˇsiˇce jako DASPK, IDA, atp.

viz kapitola 3.4). Jakoˇzto primárn´ı promˇenné mohou být pˇri pouˇzit´ı této tˇr´ıdy metod vybrány bud’ koncentrace vˇsech speci´ı nebo koncentrace komponentn´ıch speci´ı (c_j, s_j) spolu s koncentracemi sraˇzených speci´ı. Soustava (resp. matice soustavy) se- stavená postupem DAE má, v porovnán´ı s pouˇzit´ım dále pˇredstavených metod, vel- kou dimenzi a jej´ıˇreˇsen´ı m˚uˇze být ˇcasovˇe nároˇcné. Pˇres pˇr´ıpadnou ˇcasovou nároˇcnost se metoda znovu zaˇc´ıná dostávat do zájmu modeláˇr˚u v souvislosti s reakˇcnˇe orien- tovanými (tzv.

”reaction-based“) postupy (viz kapitolu 2.2.4 a zdroje [22],[12]) a ke slovu se hlás´ı také ˇreˇsiˇce soustav diferenciáln´ıch a algebraických rovnic (DAE-solvery [23] viz kapitolu 3.4). Pro ilustraci tradiˇcn´ıho mechanismu pouˇzit´ı DAE (viz [11]) následuje pˇr´ıklad.