<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 new 'errors' add-on for CO2SYS-Matlab to propagate uncertainties. It uses CO2SYS-Matlab and the add-on routine errors.m (which itself calls another add-on routine derivnum.m) in octave, GNU's clone of Matlab. One can simply inspect the HTML version of this file or take 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.
If you are visualizing this after clicking on the link to this file on github, you are visualizing the HTML version of a jupyter notebook. Alternatively, you may run cells interactively and modify them if you have jupyter notebook installed on your machine. To install that software, just download and the anaconda open software installer for your computing platform (Windows, OS X, or Linux) from https://www.continuum.io/downloads and then follow the easy install instructions at
https://docs.continuum.io/anaconda/install#
Then just download this jupyter notebook file as well as the 3 routines in the src directory (CO2SYS.m, errors.m, and derivnum.m). Afterwards, you'll only need to install octavemagic, available in the oct2py package with the 1-line command in the following section.
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 the following command:
jupyter notebook
A new window or tab should appear in your browser, showing the files in the directory from where you launched the above command. Then just click on one of the .ipynb files, such as this one.
Once the notebook file appears in your browser, move to any of its cells with your mouse. Run a cell by clicking on it and hitting Ctrl-Return. Alternatively, type Shift-Return to run a cell and then move to the next one. More information on all the interactive commankds is available in the Jupyter Notebook Quick Start Guide: http://jupyter-notebook-beginner-guide.readthedocs.io/en/latest/execute.html
At the top of the notebook, you'll see 8 Tabs (File, Edit, View, Insert, ...). Those tabls provide commands that will allow you to do whatever you like. Under the Help tab you'll find keyboard shortcuts for commands. Alternatively, a cheat sheet for short cuts to commands within jupyter notebook is available at https://zenodo.org/record/44973/files/ipynb-cheat-sheet.pdf . Or use the command palette after typing Ctrl-Shift-P.
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 ("~/Software/MATLAB/CO2SYS-MATLAB/src")
In [3]:
%%octave
help errors
In [4]:
%%octave
# Standard input for CO2SYS:
# --------------------------
# 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 = 18; % 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 for CO2SYS:
# --------------------------
# 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
#Default uncertainties (pK units): epK0 , epK1, epK2, epKb, epKw, epKspA, epKspC
epK = [0.002, 0.0075, 0.015, 0.01, 0.01, 0.02, 0.02];
#Default uncertainty for Total boron (0.01, i.e., a 1% relative uncertainty)
eBt = 0.02;
In [5]:
%%octave
% With no errors from Ks
epK = 0.0;
eBt = 0.0;
[e, ehead, enice] = errors (PAR1, PAR2, PAR1TYPE, PAR2TYPE, SAL, TEMPIN, TEMPOUT, PRESIN, PRESOUT, SI, PO4,...
ePAR1, ePAR2, eSAL, eTEMP, eSI, ePO4, epK, eBt, r, ...,
pHSCALEIN, K1K2CONSTANTS, KSO4CONSTANTS);
# Print results
# e
# ehead
# enice
# Print (nicely formatted):
printf("%s %s %s %s %s %s %s %s %s \n", ehead{1:9});
printf("%f %f %f %f %f %f %f %f %f \n", e(1:9));
In [6]:
%%octave
% With std errors from constants except for Bt
epK = [0.002, 0.0075, 0.015, 0.01, 0.01, 0.02, 0.02];
eBt = 0.0;
[ek, ekhead, eknice] = errors (PAR1, PAR2, PAR1TYPE, PAR2TYPE, SAL, TEMPIN, TEMPOUT, PRESIN, PRESOUT, SI, PO4, ...
ePAR1, ePAR2, eSAL, eTEMP, eSI, ePO4, epK, eBt, r, ...
pHSCALEIN, K1K2CONSTANTS, KSO4CONSTANTS);
printf("%s %s %s %s %s %s %s %s %s \n", ekhead{1:9});
printf("%f %f %f %f %f %f %f %f %f \n", ek(1:9));
In [7]:
%%octave
% With std errors from constants & Bt
epK = [0.002, 0.0075, 0.015, 0.01, 0.01, 0.02, 0.02];
eBt = 0.02;
[ekb, ekbhead] = errors (PAR1, PAR2, PAR1TYPE, PAR2TYPE, SAL, TEMPIN, TEMPOUT, PRESIN, PRESOUT, SI, PO4, ...
ePAR1, ePAR2, eSAL, eTEMP, eSI, ePO4, epK, eBt, r, ...
pHSCALEIN, K1K2CONSTANTS, KSO4CONSTANTS);
printf("%s %s %s %s %s %s %s %s %s \n", ekbhead{1:9});
printf("%f %f %f %f %f %f %f %f %f \n", ekb(1:9));
Same calculation, but with a simpler way to specify the defaults for epK and eBt
In [8]:
%%octave
% With std errors from constants & Bt
epK = '';
eBt = '';
[ekb, ekbhead] = errors (PAR1, PAR2, PAR1TYPE, PAR2TYPE, SAL, TEMPIN, TEMPOUT, PRESIN, PRESOUT, SI, PO4, ...
ePAR1, ePAR2, eSAL, eTEMP, eSI, ePO4, epK, eBt, r, ...
pHSCALEIN, K1K2CONSTANTS, KSO4CONSTANTS);
printf("%s %s %s %s %s %s %s %s %s \n", ekbhead{1:9});
printf("%f %f %f %f %f %f %f %f %f \n", ekb(1:9));
In [9]:
%%octave
printf("%s %s %s %s %s %s %s %s %s \n", ekbhead{1:9});
%printf("%f %f %f %f %f %f %f %f %f \n", e(1:9));
%printf("%f %f %f %f %f %f %f %f %f \n", ek(1:9));
printf("%f %f %f %f %f %f %f %f %f \n", ek(1:9) ./ e(1:9));
Conclusion: Accounting for uncertainties from the equilibrium constants increases propagated uncertainties by 1.7 to 5.7 times for the At-Ct pair with default uncertainties.
In [10]:
%%octave
printf("%s %s %s %s %s %s %s %s %s \n", ekbhead{1:9});
%printf("%f %f %f %f %f %f %f %f %f \n", ek(1:9));
%printf("%f %f %f %f %f %f %f %f %f \n", ekb(1:9));
printf("%f %f %f %f %f %f %f %f %f \n", ekb(1:9) ./ ek(1:9));
Conclusion: Accounting for uncertainty in total boron increases propagated absolute uncertainties by less than 1% when using the default uncertainty (eBt = 0.01) but by up to 3.5% when using the eBt=0.02, i.e., with the At-Ct pair.
In [11]:
%%octave
% With std errors from Ks
[d, dhead, dnice] = CO2SYS (PAR1, PAR2, PAR1TYPE, PAR2TYPE, SAL, TEMPIN, TEMPOUT, PRESIN, PRESOUT, SI, PO4, ...
pHSCALEIN, K1K2CONSTANTS, KSO4CONSTANTS);
In [12]:
%%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)
OmegaCa = d(:,15); %15 - OmegaCa input ()
OmegaAr = d(:,16); %16 - OmegaAr input ()
xCO2 = d(:,17); %17 - xCO2 input (ppm)
H = 10**(-pH) * 1e9; # Convert pH to H+ and then from mol/kg to nmol/kg
dar = [H, pCO2, fCO2, HCO3, CO3, CO2, OmegaCa, OmegaAr, xCO2];
%dar
In [13]:
%%octave
perr = 100* ekb(1:9) ./ dar;
# Show errors in percent of base value computed with CO2SYS
printf("%s %s %s %s %s %s %s %s %s \n", ehead{1:9});
printf("%4.2f %6.2f %7.2f %8.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n", perr(1:9));
An absolute change in pH is equivalent to a relative change in H+. That is, it can be shown that for small changes in H+ (dH):
\begin{equation} d\text{pH} = -\frac{1}{\text{ln}(10)} \frac{d\text{H}}{\text{H}} \end{equation}To get to this result,
1) start with the basic definition: $\text{pH} = - \text{log}_{10} \, [\text{H}^+]$,
2) convert the logarithm base 10 to a natural log with $\text{log}_{10}(x) = \text{ln}(x) / \text{ln}(10)$,
3) take the derivative of each side, and
4) plug in the definition of the derivative of a natural log, i.e., $d \text{ln}(x) = dx/x$.
In [14]:
%%octave
# We drop the minus sign below because errors are positive by definition (1 sigma)
# Note that "(0.01 * perr)" is the fractional error in hydrogen ion concentration (i.e., dH/H)
dpH = (0.01 * perr(1))/log(10);
dpH
In [15]:
%%octave
# Standard input for CO2SYS:
# --------------------------
# Input Variables:
PAR1 = 2300; % ALK
PAR2 = 8.1; % pH
PAR1TYPE = 1; % 1=ALK, 2=DIC, 3=pH, 4=pCO2, 5=fCO2
PAR2TYPE = 3; % Same 5 choices as PAR1TYPE
SAL = 35; % Salinity
TEMPIN = 18; % Temperature (input)
TEMPOUT = 18; % 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 = 1; %% pH scale (1=total, 2=seawater, 3=NBS, 4=Free)
K1K2CONSTANTS = 10; %% 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 for CO2SYS:
# --------------------------
# 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 = 0.01; % 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
#Default uncertainties: epK0 , epK1, epK2, epKb, epKw, epKspA, epKspC
epK = [0.002, 0.0075, 0.015, 0.01, 0.01, 0.02, 0.02];
#Default uncertainty for Total Boron: 0.01 (i.e., a 1% relative error)
eBt = 0.02;
In [16]:
%%octave
% With std errors from constants & Bt
epK = '';
eBt = '';
[e, ehead, eunits] = errors (PAR1, PAR2, PAR1TYPE, PAR2TYPE, SAL, TEMPIN, TEMPOUT, PRESIN, PRESOUT, SI, PO4, ...
ePAR1, ePAR2, eSAL, eTEMP, eSI, ePO4, epK, eBt, r, ...
pHSCALEIN, K1K2CONSTANTS, KSO4CONSTANTS);
printf("%s %s %s %s %s %s %s %s %s \n", ehead{1:9});
printf("%s %s %s %s %s %s %s %s %s \n", eunits{1:9});
printf("%f %f %f %f %f %f %f %f %f \n", e(1:9));
In [17]:
%%octave
PAR1 = 7.91; % pH_T
PAR2 = 500; % pCO2
PAR1TYPE = 3;
PAR2TYPE = 4;
SAL = 35;
TEMPIN = 18;
TEMPOUT = 18;
PRESIN = 0;
PRESOUT = PRESIN;
SI = 0;
PO4 = 0;
pHSCALEIN = 1; % total scale
K1K2CONSTANTS = 10; % Lueker 2000
KSO4CONSTANTS = 1; % KSO4 of Dickson & TB of Uppstrom 1979
ePAR1 = 0.02; % pH
ePAR2 = 5; % pCO2
eSAL = 0;
eTEMP = 0;
eSI = 0;
ePO4 = 0;
epK=0;
eBt=0;
r=0;
[e, ehead, eunits] = errors (PAR1,PAR2,PAR1TYPE,PAR2TYPE,SAL,TEMPIN,TEMPOUT,PRESIN,PRESOUT,SI,PO4,...
ePAR1,ePAR2,eSAL,eTEMP,eSI,ePO4,epK,eBt,r,pHSCALEIN,K1K2CONSTANTS,KSO4CONSTANTS);
printf("%s %s %s %s %s %s %s %s %s \n", ehead{1:9});
printf("%s %s %s %s %s %s %s %s %s \n", eunits{1:9});
printf("%f %f %f %f %f %f %f %f %f \n", e(1:9));
In [18]:
%%octave
ehead, eunits