<img align="left" width="60%" src="http://www.lsce.ipsl.fr/Css/img/banniere_LSCE_75.png" >
LSCE/IPSL, CEA-CNRS-UVSQ, Gif-sur-Yvette, France
Abstract: This notebook shows how to use the standard CO2SYS routine of CO2SYS-Matlab to calculate derived variables of marine CO2 system. Basic and more advanced use is illustrated. For example, for the latter it is demonstrated how to export CO2SYS results into python and then show them off in more convenient form and manipulate them more easily with python's pandas library. This CO2SYS-Matlab routine is used in octave, a GNU's clone of Matlab. You can simply inspect the HTML version of this file or try to go further by taking advantage of being able to execute its commands interactively in your browser. The latter requires jupyter notebook and its oct2py add-on that includes the octavemagic interface for octave (see below).
To install the octavemagic funtionality, we must install oct2py, with the following command at the Unix prompt:
conda install -c conda-forge oct2py=3.6.5
Then launch the notebook as usual with jupyter notebook
Details on using octavemagic are here: https://ipython.org/ipython-doc/2/config/extensions/octavemagic.html
In [1]:
%load_ext oct2py.ipython
In [2]:
%%octave
#addpath()
#addpath ("~/analysis/Projects/epitalon/ErrorPropagation")
#addpath ("~/IPython_Notebooks/epitalon/Matlab/CO2sys")
addpath ("~/Software/MATLAB/CO2SYS-MATLAB/src")
In [3]:
%octave X = [1 2 3 4] ; mean(X)
Out[3]:
In [4]:
%%octave
X = [1 2; 3 4]
mean(X)
In [5]:
%%octave
# Input Variables:
PAR1 = 2300; % ALK
PAR2 = 2000; % DIC
PAR1TYPE = 1; % 1=ALK, 2=DIC, 3=pH, 4=pCO2, 5=fCO2
PAR2TYPE = 2; % Same 5 choices as PAR1TYPE
SAL = 35; % Salinity
TEMPIN = 25; % Temperature (input)
TEMPOUT = 25; % Temperature (output)
PRESIN = 0; % Pressure (input)
PRESOUT = PRESIN; % Pressure (output)
SI = 60; % Total dissolved inorganic silicon (Sit)
PO4 = 2; % Total dissoloved inorganic Phosphorus (Pt)
# Input Parameters:
pHSCALEIN = 2; % pH scale (1=total, 2=seawater, 3=NBS, 4=Free)
K1K2CONSTANTS = 15; % set for K1 & K2: (a) 10=Lueker et al. (2000); (b) 14=Millero (2010)
KSO4CONSTANTS = 1; % KSO4 of Dickson (1990a) & Total dissolved boron (Bt) from Uppstrom (1974)
# Input variables for error propagation:
r = 0.0; % Correlation between **uncertainties** in PAR1 and PAR2 (-1 < r < 1)
ePAR1 = 2; % uncertainty in PAR1 (same units as PAR1)
ePAR2 = 2; % uncertainty in PAR2 (same units as PAR2)
eSAL = 0; % uncertainty in Salinity
eTEMP = 0; % uncertainty in Temperature
eSI = 4; % uncertainty in Sit
ePO4 = 0.1; % uncertainty in Pt
In [6]:
%%octave
[d, dhead, dnice] = CO2SYS (PAR1, PAR2, PAR1TYPE, PAR2TYPE, SAL, TEMPIN, TEMPOUT, PRESIN, PRESOUT, SI, PO4,...
pHSCALEIN, K1K2CONSTANTS, KSO4CONSTANTS);
Display the Header information to see index numbers for each computed variable
In [7]:
%%octave
dhead
Display the Nice Header information to also see units
In [8]:
%%octave
dnice
In [9]:
%%octave
K1 = d(54) ; pK1 = -log10(K1) #Check value from Waters et al. (2014) is pK1 = 5.8404 (seawater scale)
K2 = d(55) ; pK2 = -log10(K2) #Check value from Waters et al. (2014) is pK2 = 8.9662 (seawater scale)
These values should check exactly.
IMPORTANT: The constants from Waters et al. (2014) are an update of those from Millero (2010), for which inconsistencies were identified (Orr and Epitalon, 2015; Orr et al. 2015). I have tested constants from Waters et al. (2014), and they no longer show these inconsistencies. Therefore, they should be used in place of the constants from Millero (2010), which should no longer be used.
In [10]:
%%octave
pHSCALEIN = 1; % pH scale (1=total, 2=seawater, 3=NBS, 4=Free)
K1K2CONSTANTS = 10; % K1 & K2 (two best choices): 10=Lueker et al. (2000); 15=Waters et al. (2014)
KSO4CONSTANTS = 1; % KSO4 of Dickson (1990a) & Total dissolved boron (Bt) from Uppstrom (1974)
[d, dhdr, dnice] = CO2SYS (PAR1, PAR2, PAR1TYPE, PAR2TYPE, SAL, TEMPIN, TEMPOUT, PRESIN, PRESOUT, SI, PO4,...
pHSCALEIN, K1K2CONSTANTS, KSO4CONSTANTS);
Show it all
In [11]:
%%octave
d
In [12]:
%%octave
help CO2SYS
In [13]:
%%octave
help errors
Select some output
In [14]:
%%octave
dhead(1:6)
d(1:6)
A bit nicer to read
In [15]:
%%octave
printf ("%s %s %s %s %s %s\n", dhead{1:6})
printf ("%f %f %f %f %f %f\n", d(1:6))
Get all the usual output, but only at the input T and P (don't show results at output T and P)
In [16]:
%%octave
#printf ("%s %s %s %s %s %s\n", Headers(1:6))
printf ("%s %s %s %s %s %s %s %s %s %s\n", dhead{1:8}, dhead{14}, dhead{16})
printf ("%f %f %f %f %f %f %f %f %f %f\n", d(1:8), d(14), d(16))
Notes:
In [17]:
# import python's numpy library
import numpy as np
# Export d to python, renaming it D
%octave -o D D=d(1:6);
# Export dhead to python, renaming id Dhead
%octave -o Dhead Dhead=dhead;
D = np.array(D)
Dhead = np.array(Dhead)
Dhead = Dhead[0:6]
print(Dhead)
print(D)
In [18]:
import pandas as pd
#import numpy as np
pad = pd.DataFrame(D, columns=Dhead)
pad
Out[18]:
Make same calculation but at 4 temperatures
In [19]:
%%octave
TEMPIN = [0:10:30];
TEMPIN
# If any input array has more than 1 member, all others will be duplicated to reach that number
# (input arrays must have either n=1 or the same as the maximum n)
[d, dhdr, dnice] = CO2SYS (PAR1, PAR2, PAR1TYPE, PAR2TYPE, SAL, TEMPIN, TEMPOUT, PRESIN, PRESOUT, SI, PO4,...
pHSCALEIN, K1K2CONSTANTS, KSO4CONSTANTS);
# Show calculated pCO2
d(:,4)
Define output arrays knowing the CO2SYS index numbers (get these from dhdr or dnice)
In [20]:
%%octave
ALK = d(:,1); %01 - TAlk (umol/kgSW)
DIC = d(:,2); %02 - TCO2 (umol/kgSW)
pH = d(:,3); %03 - pHin ()
pCO2 = d(:,4); %04 - pCO2 input (uatm)
fCO2 = d(:,5); %05 - fCO2 input (uatm)
HCO3 = d(:,6); %06 - HCO3 input (umol/kgSW)
CO3 = d(:,7); %07 - CO3 input (umol/kgSW)
CO2 = d(:,8); %08 - CO2 input (umol/kgSW)
Hfree = d(:,13); %13 - Hfree input (umol/kgSW)
RF = d(:,14); %14 - RevelleFactor input ()
OmegaCa = d(:,15); %15 - OmegaCa input ()
OmegaAr = d(:,16); %16 - OmegaAr input ()
xCO2 = d(:,17); %17 - xCO2 input (ppm)
Display the results in a concise manner (formatted)
In [21]:
%%octave
printf('%s %s %s %s %s %s %s %s %s %s %s %s \n',
'T', 'ALK', 'DIC', 'pH', 'pCO2', 'fCO2', 'xCO2', 'HCO3', 'CO3', 'CO2', 'RF', 'OmegaAr')
for i=1 : length(pCO2)
printf('%7.3f %7.1f %7.1f %7.4f %6.1f %6.1f %6.1f, %7.1f %7.1f %6.2f %f %6.1f\n',
TEMPIN(i), ALK(i), DIC(i),
pH(i), pCO2(i), fCO2(i), xCO2(i), HCO3(i), CO3(i), CO2(i), RF(i), OmegaAr(i));
end
Write out those same results to a output file (CSV file)
In [22]:
%%octave
fileout = 'test.csv';
fid = fopen(fileout, 'w');
fprintf(fid,'%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s \n',
'T', 'ALK', 'DIC', 'pH', 'pCO2', 'fCO2', 'xCO2', 'HCO3', 'CO3', 'CO2', 'RF', 'OmegaAr');
for i=1 : length(pCO2)
fprintf(fid,'%f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f\n',
TEMPIN(i), ALK(i), DIC(i),
pH(i), pCO2(i), fCO2(i), xCO2(i), HCO3(i), CO3(i), CO2(i), RF(i), OmegaAr(i));
end
fclose('all');
In [23]:
#!echo $PWD
!cat test.csv
Following cell is in python (default), because it does not start with %%octave
In [24]:
df=pd.read_csv('test.csv', sep=',')
df
Out[24]:
In [25]:
!cat Alk-pH-P-vary-for-CO2sys-Matlab.csv
In [26]:
%%octave
filein = 'Alk-pH-P-vary-for-CO2sys-Matlab.csv';
datin = csvread(filein, 1, 0);
# Transpose array so that columns become rows (and vice versa)
datin = transpose(datin);
par1 = datin(1,:);
par2 = datin(2,:);
par1type = datin(3,:);
par2type = datin(4,:);
sal = datin(5,:);
tempin = datin(6,:);
presin = datin(7,:);
# Specify remaining input data
tempout = 0; % not used
presout = 0; % not used
sil = 0; % Concentration of Sit in the sample (in umol/kg)
po4 = 0; % Concentration of Pt in the sample (in umol/kg)
pHscale = 1; % pH scale at which the input pH is reported ("1" means "Total Scale")
k1k2c = 10; % K1 & K2 ("10" means "Lueker et al., 2000")
kso4c = 1; % KSO4 & Bt ("1" means "KSO4 from Dickson 1990a & Bt from Uppstrom 1974")
b = CO2SYS(par1, par2, par1type, par2type, sal, tempin, tempout, presin,presout,sil,po4,pHscale,k1k2c,kso4c);
In [27]:
%%octave
ALK = b(:,1); %01 - TAlk (umol/kgSW)
DIC = b(:,2); %02 - TCO2 (umol/kgSW)
pH = b(:,3); %03 - pHin ()
pCO2 = b(:,4); %04 - pCO2 input (uatm)
fCO2 = b(:,5); %05 - fCO2 input (uatm)
HCO3 = b(:,6); %06 - HCO3 input (umol/kgSW)
CO3 = b(:,7); %07 - CO3 input (umol/kgSW)
CO2 = b(:,8); %08 - CO2 input (umol/kgSW)
Hfree = b(:,13); %13 - Hfree input (umol/kgSW)
RF = b(:,14); %14 - RevelleFactor input ()
OmegaCa = b(:,15); %15 - OmegaCa input ()
OmegaAr = b(:,16); %16 - OmegaAr input ()
xCO2 = b(:,17); %17 - xCO2 input (ppm)
In [28]:
%%octave
printf('%s %s %s %s %s %s %s %s %s %s %s %s \n',
'T', 'ALK', 'DIC', 'pH', 'pCO2', 'fCO2', 'xCO2', 'HCO3', 'CO3', 'CO2', 'RF', 'OmegaAr')
for i=1 : length(pCO2)
printf('%7.3f, %7.1f %7.1f %7.4f %6.1f %6.1f %6.1f, %7.1f %7.1f %6.2f %7.3f %6.3f\n',
tempin(i), ALK(i), DIC(i),
pH(i), pCO2(i), fCO2(i), xCO2(i), HCO3(i), CO3(i), CO2(i), RF(i), OmegaAr(i));
end
Pull arrays, then combine into a pandas DataFrame
In [29]:
%octave_pull tempin ALK DIC pH pCO2 fCO2 xCO2 HCO3 CO3 CO2 RF OmegaAr
#print(tempin)
pieces = [ALK, DIC, pH, pCO2, fCO2, xCO2, HCO3, CO3, CO2, RF, OmegaAr]
pieces = np.array(pieces)
pieces = np.squeeze(pieces)
pieces = np.transpose(pieces)
cols = list(['ALK', 'DIC', 'pH', 'pCO2', 'fCO2', 'xCO2', 'HCO3', 'CO3', 'CO2', 'RF', 'OmegaAr'])
#pad = pd.DataFrame(pieces, index=list(tempin), columns=cols)
df = pd.DataFrame(pieces, columns=cols)
df
Out[29]:
Quick statistical summary
In [30]:
df.describe()
Out[30]:
Select results where $\Omega_{Ar} < 1$
In [31]:
df[df.OmegaAr < 1]
Out[31]:
For more information on pandas (in Python) see for example http://pandas.pydata.org/pandas-docs/stable/10min.html
Plot sensitivity of pH and pCO2 to changes in DIC, while keeping everything else constant.
(ALK=2300, Si=50, PO4=2, dissociation constats: Mehrbach Refit)')
In [32]:
%%octave -s 1000,400
par1type = 1; % The first parameter supplied is of type "1", which is "alkalinity"
par1 = 2300; % value of the first parameter
par2type = 2; % The first parameter supplied is of type "1", which is "DIC"
par2 = [2000:5:2200]; % value of the second parameter, which is a long vector of different DIC's!
sal = 35; % Salinity of the sample
tempin = 10; % Temperature at input conditions
presin = 0; % Pressure at input conditions
tempout = 0; % Temperature at output conditions - doesn't matter in this example
presout = 0; % Pressure at output conditions - doesn't matter in this example
sil = 50; % Concentration of silicate in the sample (in umol/kg)
po4 = 2; % Concentration of phosphate in the sample (in umol/kg)
pHscale = 1; % pH scale at which the input pH is reported ("1" means "Total Scale") - doesn't matter in this example
k1k2c = 4; % Choice of dissociation constants K1 and K2 ("4" means "Mehrbach refit")
kso4c = 1; % Choice of KSO4 ("1" means "Dickson")
% Do the calculation. See CO2SYS's help for syntax and output format
A = CO2SYS(par1,par2,par1type,par2type,sal,tempin,tempout,presin,presout,sil,po4,pHscale,k1k2c,kso4c);
figure; clf
subplot(1,2,1)
plot(par2,A(:,4),'r.-') # The calculated pCO2's are in the 4th column of the output A of CO2SYS
xlabel('DIC'); ylabel('pCO2 [uatm]')
subplot(1,2,2)
plot(par2,A(:,3),'b.-') # The calculated pH's are in the 3rd column of the output A of CO2SYS
xlabel('DIC'); ylabel('pH')
Millero, F. J.: Carbonate constants for estuarine waters, Mar. Freshwater Res., 61, 139–142, doi:10.1071/MF09254, 2010.
Orr, J. C. and Epitalon, J.-M.: Improved routines to model the ocean carbonate system: mocsy 2.0, Geosci. Model Dev., 8, 485–499, doi:10.5194/gmd-8-485-2015, 2015.
Orr, J. C., Gattuso, J.-P., and Epitalon, J.-M.: Comparison of ten packages that compute ocean carbonate chemistry, Biogeosciences, 12, 1483–1510, doi:10.5194/bg-12-1483-2015, 2015.
Waters, J., Millero, F., and Woosley, R.: Corrigendum to “The free proton concentration scale for seawater pH”,[MARCHE: 149 (2013) 8-22], Mar. Chem., 165, 66–67, 2014.