-
Notifications
You must be signed in to change notification settings - Fork 21
/
Copy pathcarbon_chemistry.jl
213 lines (171 loc) · 7.91 KB
/
carbon_chemistry.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
using Roots
using OceanBioME.Models: teos10_polynomial_approximation
struct CarbonChemistry{P0, PC, PB, PS, PF, PP, PSi, PW, IS, PKS, PRho}
ionic_strength :: IS
solubility :: P0
carbonic_acid :: PC
boric_acid :: PB
water :: PW
sulfate :: PS
fluoride :: PF
phosphoric_acid :: PP
silicic_acid :: PSi
calcite_solubility :: PKS
density_function :: PRho
end
"""
CarbonChemistry(FT = Float64;
ionic_strength = IonicStrength(),
solubility = K0(),
carbonic_acid = (K1 = K1(), K2 = K2()),
boric_acid = KB(),
water = KW(),
sulfate = KS(; ionic_strength),
fluoride = KF(; ionic_strength),
phosphoric_acid = (KP1 = KP1(), KP2 = KP2(), KP3 = KP3()),
silicic_acid = KSi(; ionic_strength),
calcite_solubility = KSP_calcite(),
density_function = teos10_polynomial_approximation)
Carbon chemistry model capable of solving for sea water pCO₂ from DIC and
total alkalinity or DIC and pH.
Default form from Dickson, A.G., Sabine, C.L. and Christian, J.R. (2007),
Guide to Best Practices for Ocean CO 2 Measurements. PICES Special Publication 3, 191 pp.
See each parameters documentation for origional sources.
Example
=======
```jldoctest
julia> using OceanBioME
julia> carbon_chemistry = CarbonChemistry()
`CarbonChemistry` model which solves for pCO₂ and pH
julia> pCO₂ = carbon_chemistry(; DIC = 2000, Alk = 2000, T = 10, S = 35)
1308.0843992121615
julia> pH = carbon_chemistry(; DIC = 2000, Alk = 2000, T = 10, S = 35, return_pH = true)
7.502534641304366
julia> pCO₂_higher_pH = carbon_chemistry(; DIC = 2000, T = 10, S = 35, pH = 7.5)
1315.6558976217746
```
"""
function CarbonChemistry(FT = Float64;
ionic_strength = IonicStrength{FT}(),
solubility = K0{FT}(),
carbonic_acid = (K1 = K1(FT), K2 = K2(FT)),
boric_acid = KB(FT),
water = KW(FT),
sulfate = KS(FT; ionic_strength),
fluoride = KF(FT; ionic_strength),
phosphoric_acid = (KP1 = KP1(FT), KP2 = KP2(FT), KP3 = KP3(FT)),
silicic_acid = KSi(FT; ionic_strength),
calcite_solubility = KSP_calcite(FT),
density_function = teos10_polynomial_approximation) # the denisity function *is* going to cause type instability but I can't see a way to fix it
return CarbonChemistry(ionic_strength, solubility, carbonic_acid, boric_acid, water,
sulfate, fluoride, phosphoric_acid, silicic_acid, calcite_solubility, density_function)
end
"""
alkalinity_residual(H, p)
Returns the difference between total alkalinity computed from `H`` (hydrogen ion
concentration), `DIC`, `borate`, `sulfate`, `phosphate`, `silicate`, and `fluoride`
concentration and chemical equilibrium constants specified in `p`, and the specified
total `Alk`alinity.
TAlk = [HCO₃⁻] + 2[CO₃²⁻] + [B(OH)₄⁻] + [OH⁻] + [HPO₄²⁻] + 2[PO₄³⁻] + [SiO(OH)₃⁻]
+ [NH₃] + [HS⁻] - [H⁺] - [HSO₄⁻] - [HF] - [H₃PO₄] + minor acids and bases
Concentrations diagnosed as specified in Dickson et. al best practice descried in
`CarbonChemistry` docstring.
Note ammonia (NH₃) is not currently included.
"""
@inline function alkalinity_residual(H, p)
carbonate_denom = H^2 + p.K1 * H + p.K1 * p.K2
phosphorus_denom = H^3 + p.KP1 * H^2 + p.KP1 * p.KP2 * H + p.KP1 * p.KP2 * p.KP3
sulfate_denom = 1 + p.sulfate / p.KS
bicarbonate = p.K1 * H * p.DIC / carbonate_denom
carbonate = 2 * p.DIC * p.K1 * p.K2 / carbonate_denom
borate = p.boron / (1 + H / p.KB)
hydroxide = p.KW / H
hydrogen_phosphate = p.phosphate * p.KP1 * p.KP2 * H / phosphorus_denom
phosphate = 2 * p.phosphate * p.KP1 * p.KP2 * p.KP3 / phosphorus_denom
silicate = p.silicate / (1 + H / p.KSi)
free_hydrogen = - H / sulfate_denom
hydrogen_suplfate = - p.sulfate / (1 + p.KS / H / sulfate_denom)
hydrogen_fluoride = -p.fluoride / (1 + p.KF / H)
phosphoric_acid = -p.phosphate * H^3 / phosphorus_denom
return (bicarbonate
+ carbonate
+ borate
+ hydroxide
+ hydrogen_phosphate
+ phosphate
+ silicate
+ free_hydrogen
+ hydrogen_suplfate
+ hydrogen_fluoride
+ phosphoric_acid
- p.Alk)
end
"""
(p::CarbonChemistry)(; DIC, T, S, Alk = 0, pH = nothing,
return_pH = false,
boron = 0.000232 / 10.811 * S / 1.80655,
sulfate = 0.14 / 96.06 * S / 1.80655,
fluoride = 0.000067 / 18.9984 * S / 1.80655,
silicate = 0,
phosphate = 0,
upper_pH_bound = 14,
lower_pH_bound = 0)
Calculates `fCO₂` in sea water with `DIC`, `Alk`alinity, `T`emperature, and `S`alinity
unless `pH` is specified, in which case intermediate computation of `pH` is skipped and
`pCO₂` is calculated from the `DIC`, `T`, `S` and `pH`.
Alternativly, `pH` is returned if `return_pH` is `true`.
"""
@inline function (p::CarbonChemistry)(; DIC, T, S, Alk = 0, pH = nothing,
P = nothing,
lon = 0,
lat = 0,
return_pH = false,
boron = 0.000232 / 10.811 * S / 1.80655,
sulfate = 0.14 / 96.06 * S / 1.80655,
fluoride = 0.000067 / 18.9984 * S / 1.80655,
silicate = 0,
phosphate = 0,
upper_pH_bound = 14,
lower_pH_bound = 0)
ρₒ = p.density_function(T, S, ifelse(isnothing(P), 0, P), lon, lat)
# Centigrade to kelvin
T += 273.15
# mili-equivalents / m³ to equivalents / kg
Alk *= 1e-3 / ρₒ
# mmol / m³ to mol / kg
DIC *= 1e-3 / ρₒ
phosphate *= 1e-3 / ρₒ
silicate *= 1e-3 / ρₒ
# ionic strength
Is = p.ionic_strength(S)
# compute equilibrium constants
K1 = p.carbonic_acid.K1(T, S; P)
K2 = p.carbonic_acid.K2(T, S; P)
KB = p.boric_acid(T, S; P)
KW = p.water(T, S; P)
KS = p.sulfate(T, S, Is; P)
KF = p.fluoride(T, S, Is, KS; P)
KP1 = p.phosphoric_acid.KP1(T, S; P)
KP2 = p.phosphoric_acid.KP2(T, S; P)
KP3 = p.phosphoric_acid.KP3(T, S; P)
KSi = p.silicic_acid(T, S, Is)
params = (; DIC, Alk, boron, sulfate, fluoride, silicate, phosphate,
K1, K2, KB, KW, KS, KF, KP1, KP2, KP3, KSi)
# solve equilibrium for hydrogen ion concentration
H = solve_for_H(pH, params, upper_pH_bound, lower_pH_bound)
# compute solubility equilibrium constant
K0 = p.solubility(T, S)
# compute pCO₂
CO₂ = DIC * H ^ 2 / (H ^ 2 + K1 * H + K1 * K2)
fCO₂ = (CO₂ / K0) * 10 ^ 6 # μatm
# compute pH
pH = -log10(H)
return ifelse(return_pH, pH, fCO₂)
end
# solves `alkalinity_residual` for pH
solve_for_H(pH, args...) = 10.0 ^ - pH
solve_for_H(::Nothing, params, upper_pH_bound, lower_pH_bound) =
find_zero(alkalinity_residual, (10.0 ^ - upper_pH_bound, 10.0 ^ - lower_pH_bound), Bisection(); atol = 1e-10, p = params)
# display
summary(::IO, ::CarbonChemistry) = string("`CarbonChemistry` model")
show(io::IO, ::CarbonChemistry) = print(io, "`CarbonChemistry` model which solves for pCO₂ and pH")