forked from ooici/ion-functions
-
Notifications
You must be signed in to change notification settings - Fork 12
/
co2_functions.py
481 lines (364 loc) · 16.5 KB
/
co2_functions.py
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
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
#!/usr/bin/env python
"""
@package ion_functions.data.pco2_functions
@file ion_functions/data/pco2_functions.py
@author Christopher Wingard
@brief Module containing CO2 instrument family related functions
"""
import numpy as np
import numexpr as ne
from ion_functions.utils import fill_value
# wrapper functions to extract parameters from SAMI-II CO2 instruments (PCO2W)
# and process these extracted parameters to calculate pCO2
def pco2_abs434_ratio(light):
"""
Description:
Extract the absorbance ratio at 434 nm from the pCO2 instrument light
measurements. This will extract the CO2ABS1_L0 data product from the
instrument light measurement arrays.
Implemented by:
2013-04-20: Christopher Wingard. Initial code.
2014-02-19: Christopher Wingard. Updated comments.
Usage:
a434ratio = pco2_abs434_ratio(light)
where
a434ratio = optical absorbance Ratio at 434 nm (CO2ABS1_L0) [unitless]
light = array of light measurements
References:
OOI (2012). Data Product Specification for Partial Pressure of CO2 in
Seawater. Document Control Number 1341-00510.
https://alfresco.oceanobservatories.org/ (See: Company Home >>
OOI >> Controlled >> 1000 System Level >>
1341-00490_Data_Product_SPEC_PCO2WAT_OOI.pdf)
"""
light = np.atleast_2d(light)
a434ratio = light[:, 6]
return a434ratio
def pco2_abs620_ratio(light):
"""
Description:
Extract the absorbance ratio at 620 nm from the pCO2 instrument light
measurements. This will extract the CO2ABS2_L0 data product from the
instrument light measurement arrays.
Implemented by:
2013-04-20: Christopher Wingard. Initial code.
2014-02-19: Christopher Wingard. Updated comments.
Usage:
a620ratio = pco2_abs620_ratio(light)
where
a620ratio = optical absorbance Ratio at 620 nm (CO2ABS2_L0) [unitless]
light = array of light measurements
References:
OOI (2012). Data Product Specification for Partial Pressure of CO2 in
Seawater. Document Control Number 1341-00510.
https://alfresco.oceanobservatories.org/ (See: Company Home >>
OOI >> Controlled >> 1000 System Level >>
1341-00490_Data_Product_SPEC_PCO2WAT_OOI.pdf)
"""
light = np.atleast_2d(light)
a620ratio = light[:, 7]
return a620ratio
def pco2_blank(raw_blank):
"""
Description:
Calculates the absorbance blank at 434 or 620 nm from the SAMI2-pCO2
instrument.
Implemented by:
2013-04-20: Christopher Wingard. Initial code.
2014-02-19: Christopher Wingard. Updated comments.
2014-02-28: Christopher Wingard. Updated to except raw blank values
from a sparse array.
2018-03-04: Christopher Wingard. Updated to correctly calculate the
blank based on new code from the vendor.
Usage:
blank = pco2_blank(raw_blank)
where
blank = optical absorbance blank at 434 or 620 nm [unitless]
raw_blank = raw optical absorbance blank at 434 or 620 nm [counts]
References:
OOI (2012). Data Product Specification for Partial Pressure of CO2 in
Seawater. Document Control Number 1341-00510.
https://alfresco.oceanobservatories.org/ (See: Company Home >>
OOI >> Controlled >> 1000 System Level >>
1341-00490_Data_Product_SPEC_PCO2WAT_OOI.pdf)
"""
blank = raw_blank / 16384.
return blank
def pco2_thermistor(traw, sami_bits):
"""
Description:
Convert the thermistor data from counts to degrees Centigrade from the
pCO2 instrument.
Implemented by:
2013-04-20: Christopher Wingard. Initial code.
2014-02-19: Christopher Wingard. Updated comments.
2023-01-12: Mark Steiner. Add sami_bits arg to handle hardware upgrades
Usage:
therm = pco2_thermistor(traw, sami_bits)
where
therm = converted thermistor temperature [degC]
traw = raw thermistor temperature (CO2THRM_L0) [counts]
sami_bits = number of bits on the SAMI hardware
References:
OOI (2012). Data Product Specification for Partial Pressure of CO2 in
Seawater. Document Control Number 1341-00510.
https://alfresco.oceanobservatories.org/ (See: Company Home >>
OOI >> Controlled >> 1000 System Level >>
1341-00490_Data_Product_SPEC_PCO2WAT_OOI.pdf)
"""
# reset inputs to arrays
traw = np.atleast_1d(traw)
sami_bits = np.atleast_1d(sami_bits)
# convert raw thermistor readings from counts to degrees Centigrade
# conversion depends on whether the SAMI is older 12 bit or newer 14 bit hardware
if sami_bits[0] == 14:
Rt = ne.evaluate('log((traw / (16384. - traw)) * 17400.)')
else:
Rt = ne.evaluate('log((traw / (4096. - traw)) * 17400.)')
InvT = ne.evaluate('0.0010183 + 0.000241 * Rt + 0.00000015 * Rt**3')
therm = ne.evaluate('(1 / InvT) - 273.15')
return therm
def pco2_battery(braw, sami_bits):
"""
Description:
Convert the battery voltage data from counts to volts from the
pCO2 instrument.
Implemented by:
2023-02-23: Mark Steiner. Initial code.
Usage:
volts = pco2_battery_voltage(vraw, sami_bits)
where
volts = converted battery voltage [degC]
braw = raw battery voltage [counts]
sami_bits = number of bits on the SAMI hardware
"""
# reset inputs to arrays
braw = np.atleast_1d(braw)
sami_bits = np.atleast_1d(sami_bits)
# convert raw battery readings from counts to Volts
if sami_bits[0] == 14:
volts = ne.evaluate('braw * 3. / 4000.')
else:
volts = ne.evaluate('braw * 15. / 4096.')
return volts
def pco2_pco2wat(mtype, light, therm, ea434, eb434, ea620, eb620,
calt, cala, calb, calc, a434blank, a620blank):
"""
Description:
Function to calculate the L1 PCO2WAT core data from the pCO2 instrument
if the measurement type is 4 (pCO2 measurement), otherwise it is a
blank and return a fill value.
Implemented by:
2013-04-20: Christopher Wingard. Initial code.
2014-02-19: Christopher Wingard. Updated comments.
2014-03-19: Christopher Wingard. Optimized using feedback provided by
Chris Fortin.
2017-04-04: Pete Cable. Updated algorithm to use thermistor/blank counts
as indicated in the DPS and the usage below.
Usage:
pco2 = pco2_pco2wat(mtype, light, therm, ea434, eb434, ea620, eb620,
calt, cala, calb, calc, a434blank, a620blank)
where
pco2 = measured pco2 in seawater (PCO2WAT_L1) [uatm]
mtype = measurement type, where 4 == actual measurement and 5 == a
blank measurement [unitless]
light = array of light measurements
therm = PCO2W thermistor temperature (CO2THRM_L0) [counts]
ea434 = Reagent specific calibration coefficient
eb434 = Reagent specific calibration coefficient
ea620 = Reagent specific calibration coefficient
eb620 = Reagent specific calibration coefficient
calt = Instrument specific calibration coefficient for temperature
cala = Instrument specific calibration coefficient for the pCO2 measurements
calb = Instrument specific calibration coefficient for the pCO2 measurements
calc = Instrument specific calibration coefficient for the pCO2 measurements
a434blank = Blank measurements at 434 nm (CO2ABS1_L0) [counts]
a620blank = Blank measurements to 620 nm (CO2ABS2_L0) [counts]
References:
OOI (2012). Data Product Specification for Partial Pressure of CO2 in
Seawater. Document Control Number 1341-00510.
https://alfresco.oceanobservatories.org/ (See: Company Home >>
OOI >> Controlled >> 1000 System Level >>
1341-00490_Data_Product_SPEC_PCO2WAT_OOI.pdf)
"""
# reset inputs to arrays
# measurements
mtype = np.atleast_1d(mtype)
light = np.atleast_2d(light)
therm = np.atleast_1d(therm)
# calibration coefficients
ea434 = np.atleast_1d(ea434)
eb434 = np.atleast_1d(eb434)
ea620 = np.atleast_1d(ea620)
eb620 = np.atleast_1d(eb620)
calt = np.atleast_1d(calt)
cala = np.atleast_1d(cala)
calb = np.atleast_1d(calb)
calc = np.atleast_1d(calc)
# blank measurements
a434blank = np.atleast_1d(a434blank)
a620blank = np.atleast_1d(a620blank)
# calculate the pco2 value
pco2 = pco2_calc_pco2(light, therm, ea434, eb434, ea620, eb620,
calt, cala, calb, calc, a434blank, a620blank)
# reset dark measurements to the fill value
m = np.where(mtype == 5)[0]
pco2[m] = fill_value
return pco2
# L1a PCO2WAT calculation
def pco2_calc_pco2(light, therm, ea434, eb434, ea620, eb620,
calt, cala, calb, calc, a434blank, a620blank):
"""
Description:
OOI Level 1 Partial Pressure of CO2 (pCO2) in seawater core data
product, which is calculated from the Sunburst SAMI-II CO2 instrument
(PCO2W).
Implemented by:
20??-??-??: J. Newton (Sunburst Sensors, LLC). Original Matlab code.
2013-04-20: Christopher Wingard. Initial python code.
2014-02-19: Christopher Wingard. Updated comments.
2014-03-19: Christopher Wingard. Optimized.
2018-03-04: Christopher Wingard. Updated to correctly calculate pCO2 using
newly formulated code provided by the vendor. Original vendor code
incorrectly calculated the blank correction. Applies additional
corrections to calculations to avoid errors thrown when running a
blank measurement.
2023-01-12: Mark Steiner. Arg therm in degrees C instead of counts
Usage:
pco2 = pco2_pco2wat(light, therm, ea434, eb434, ea620, eb620,
calt, cala, calb, calc, a434blank, a620blank)
where
pco2 = measured pco2 in seawater (PCO2WAT_L1) [uatm]
light = array of light measurements
therm = PCO2W thermistor temperature (CO2THRM_L1) [degrees C]
ea434 = Reagent specific calibration coefficient
eb434 = Reagent specific calibration coefficient
ea620 = Reagent specific calibration coefficient
eb620 = Reagent specific calibration coefficient
calt = Instrument specific calibration coefficient for temperature
cala = Instrument specific calibration coefficient for the pCO2 measurements
calb = Instrument specific calibration coefficient for the pCO2 measurements
calc = Instrument specific calibration coefficient for the pCO2 measurements
a434blank = Blank measurements at 434 nm (CO2ABS1_L0) [counts]
a620blank = Blank measurements to 620 nm (CO2ABS2_L0) [counts]
References:
OOI (2012). Data Product Specification for Partial Pressure of CO2 in
Seawater. Document Control Number 1341-00510.
https://alfresco.oceanobservatories.org/ (See: Company Home >>
OOI >> Controlled >> 1000 System Level >>
1341-00490_Data_Product_SPEC_PCO2WAT_OOI.pdf)
"""
# set constants -- original vendor formulation, reset below
# ea434 = ea434 - 29.3 * calt
# eb620 = eb620 - 70.6 * calt
# e1 = ea620 / ea434
# e2 = eb620 / ea434
# e3 = eb434 / ea434
# set the e constants, values provided by Sunburst
e1 = 0.0043
e2 = 2.136
e3 = 0.2105
# Extract variables from light array
Ratio434 = light[:, 6] # 434nm Ratio
Ratio620 = light[:, 7] # 620nm Ratio
# correct the absorbance ratios using the blanks
AR434 = (Ratio434 / a434blank)
AR620 = (Ratio620 / a620blank)
# map out blank measurements and spoof the ratios to avoid throwing an error
m = np.where(AR434 == AR620)[0]
AR434[m] = 0.99999
AR620[m] = 0.99999
# Calculate the final absorbance ratio
A434 = -1 * np.log10(AR434) # 434 absorbance
A620 = -1 * np.log10(AR620) # 620 absorbance
Ratio = A620 / A434 # Absorbance ratio
# calculate pCO2
V1 = Ratio - e1
V2 = e2 - e3 * Ratio
RCO21 = -1 * np.log10(V1 / V2)
RCO22 = (therm - calt) * 0.008 + RCO21
Tcoeff = 0.0075778 - 0.0012389 * RCO22 - 0.00048757 * RCO22**2
Tcor_RCO2 = RCO21 + Tcoeff * (therm - calt)
pco2 = 10.**((-1. * calb + (calb**2 - (4. * cala * (calc - Tcor_RCO2)))**0.5) / (2. * cala))
pco2[m] = fill_value # reset the blanks captured earlier to a fill value
return np.real(pco2)
def pco2_ppressure(xco2, p, std=1013.25):
"""
Description:
OOI Level 1 Partial Pressure of CO2 in Air (PCO2ATM) or Surface
Seawater (PCO2SSW) core date product is computed by using an
equation that incorporates the Gas Stream Pressure (PRESAIR) and the
CO2 Mole Fraction in Air or Surface Seawater (XCO2ATM or XCO2SSW,
respectively). It is computed using data from the pCO2 air-sea (PCO2A)
family of instruments.
Implemented by:
2014-10-27: Christopher Wingard. Initial python code.
Usage:
ppres = pco2_ppressure(xco2, p, std)
where
ppres = partial pressure of CO2 in air or surface seawater [uatm]
(PCO2ATM_L1 or PCO2SSW_L1)
xco2 = CO2 mole fraction in air or surface seawater [ppm] (XCO2ATM_LO
or XCO2SSW_L0)
p = gas stream pressure [mbar] (PRESAIR_L0)
std = standard atmospheric pressure set to default of 1013.25 [mbar/atm]
References:
OOI (2012). Data Product Specification for Partial Pressure of CO2 in
Air and Surface Seawater. Document Control Number 1341-00260.
https://alfresco.oceanobservatories.org/ (See: Company Home >>
OOI >> Controlled >> 1000 System Level >>
1341-00260_Data_Product_SPEC_PCO2ATM_PCO2SSW_OOI.pdf)
"""
ppres = ne.evaluate('xco2 * p / std')
return ppres
def pco2_co2flux(pco2w, pco2a, u10, t, s):
"""
Description:
OOI Level 2 core date product CO2FLUX is an estimate of the CO2 flux
from the ocean to the atmosphere. It is computed using data from the
pCO2 air-sea (PCO2A) and bulk meteorology (METBK) families of
instruments.
Implemented by:
2012-03-28: Mathias Lankhorst. Original Matlab code.
2013-04-20: Christopher Wingard. Initial python code.
Usage:
flux = pco2_co2flux(pco2w, pco2a, u10, t, s)
where
flux = estimated flux of CO2 from the ocean to atmosphere [mol m-2 s-1]
(CO2FLUX_L2)
pco2w = partial pressure of CO2 in sea water [uatm] (PCO2SSW_L1)
pco2a = partial pressure of CO2 in air [uatm] (PCO2ATM_L1)
u10 = normalized wind speed at 10 m height from METBK [m s-1] (WIND10M_L2)
t = sea surface temperature from METBK [deg_C] (TEMPSRF_L1)
s = sea surface salinity from METBK [psu] (SALSURF_L2)
References:
OOI (2012). Data Product Specification for Flux of CO2 into the
Atmosphere. Document Control Number 1341-00270.
https://alfresco.oceanobservatories.org/ (See: Company Home >>
OOI >> Controlled >> 1000 System Level >>
1341-00270_Data_Product_SPEC_CO2FLUX_OOI.pdf)
"""
# convert micro-atm to atm
pco2a = pco2a / 1.0e6
pco2w = pco2w / 1.0e6
# Compute Schmidt number (after Wanninkhof, 1992, Table A1)
Sc = 2073.1 - (125.62 * t) + (3.6276 * t**2) - (0.043219 * t**3)
# Compute gas transfer velocity (after Sweeney et al. 2007, Fig. 3 and Table 1)
k = 0.27 * u10**2 * np.sqrt(660.0 / Sc)
# convert cm h-1 to m s-1
k = k / (100.0 * 3600.0)
# Compute the absolute temperature
T = t + 273.15
# Compute solubility (after Weiss 1974, Eqn. 12 and Table I).
# Note that there are two versions, one for units per volume and
# one per mass. Here, the volume version is used.
# mol atm-1 m-3
T100 = T / 100
K0 = 1000 * np.exp(-58.0931 + (90.5069 * (100/T)) + (22.2940 * np.log(T100)) +
s * (0.027766 - (0.025888 * T100) + (0.0050578 * T100**2)))
# mol atm-1 kg-1
#K0 = np.exp(-60.2409 + (93.4517 * (100/T)) + (23.3585 * np.log(T100)) +
# s * (0.023517 - (0.023656 * T100) + (0.0047036 * T100**2)))
# Compute flux (after Wanninkhof, 1992, eqn. A2)
flux = k * K0 * (pco2w - pco2a)
return flux