-
Notifications
You must be signed in to change notification settings - Fork 11
/
CO2SYS.m
1914 lines (1810 loc) · 85.9 KB
/
CO2SYS.m
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
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
function [DATA,HEADERS,NICEHEADERS]=CO2SYS(PAR1,PAR2,PAR1TYPE,PAR2TYPE,SAL,TEMPIN,TEMPOUT,PRESIN,PRESOUT,SI,PO4,pHSCALEIN,K1K2CONSTANTS,KSO4CONSTANTS);
%**************************************************************************
%
% First CO2SYS.m version: 1.1 (Sep 2011)
% Current CO2SYS.m version: 2.0 (20 Dec 2016)
%
% CO2SYS is a MATLAB-version of the original CO2SYS for DOS.
% CO2SYS calculates and returns the state of the carbonate system of
% oceanographic water samples, if supplied with enough input.
%
% Please note that this software is intended to be exactly identical to the
% DOS and Excel versions that have been released previously, meaning that
% results obtained should be very nearly indentical for identical input.
% Additionally, several of the dissociation constants K1 and K2 that have
% been published since the original DOS version was written are implemented.
% For a complete list of changes since version 1.0, see below.
%
% For much more info please have a look at:
% Lewis, E., and D. W. R. Wallace. 1998. Program Developed for
% CO2 System Calculations. ORNL/CDIAC-105. Carbon Dioxide Information
% Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy,
% Oak Ridge, Tennessee.
% http://cdiac.ornl.gov/oceans/co2rprt.html
%
%**************************************************************************
%
% **** SYNTAX:
% [RESULT,HEADERS,NICEHEADERS]=CO2SYS(PAR1,PAR2,PAR1TYPE,PAR2TYPE,...
% ...SAL,TEMPIN,TEMPOUT,PRESIN,PRESOUT,SI,PO4,pHSCALEIN,...
% ...K1K2CONSTANTS,KSO4CONSTANTS)
%
% **** SYNTAX EXAMPLES:
% [Result] = CO2SYS(2400,2200,1,2,35,0,25,4200,0,15,1,1,4,1)
% [Result,Headers] = CO2SYS(2400, 8,1,3,35,0,25,4200,0,15,1,1,4,1)
% [Result,Headers,Niceheaders] = CO2SYS( 500, 8,5,3,35,0,25,4200,0,15,1,1,4,1)
% [A] = CO2SYS(2400,2000:10:2400,1,2,35,0,25,4200,0,15,1,1,4,1)
% [A] = CO2SYS(2400,2200,1,2,0:1:35,0,25,4200,0,15,1,1,4,1)
% [A] = CO2SYS(2400,2200,1,2,35,0,25,0:100:4200,0,15,1,1,4,1)
%
% **** APPLICATION EXAMPLE (copy and paste this into command window):
% tmps=0:40; sals=0:40; [X,Y]=meshgrid(tmps,sals);
% A = CO2SYS(2300,2100,1,2,Y(:),X(:),nan,0,nan,1,1,1,9,1);
% Z=nan(size(X)); Z(:)=A(:,4); figure; contourf(X,Y,Z,20); caxis([0 1200]); colorbar;
% ylabel('Salinity [psu]'); xlabel('Temperature [degC]'); title('Dependence of pCO2 [uatm] on T and S')
%
%**************************************************************************
%
% INPUT:
%
% PAR1 (some unit) : scalar or vector of size n
% PAR2 (some unit) : scalar or vector of size n
% PAR1TYPE () : scalar or vector of size n (*)
% PAR2TYPE () : scalar or vector of size n (*)
% SAL () : scalar or vector of size n
% TEMPIN (degr. C) : scalar or vector of size n
% TEMPOUT (degr. C) : scalar or vector of size n
% PRESIN (dbar) : scalar or vector of size n
% PRESOUT (dbar) : scalar or vector of size n
% SI (umol/kgSW) : scalar or vector of size n
% PO4 (umol/kgSW) : scalar or vector of size n
% pHSCALEIN : scalar or vector of size n (**)
% K1K2CONSTANTS : scalar or vector of size n (***)
% KSO4CONSTANTS : scalar or vector of size n (****)
%
% (*) Each element must be an integer,
% indicating that PAR1 (or PAR2) is of type:
% 1 = Total Alkalinity
% 2 = DIC
% 3 = pH
% 4 = pCO2
% 5 = fCO2
%
% (**) Each element must be an integer,
% indicating that the pH-input (PAR1 or PAR2, if any) is at:
% 1 = Total scale
% 2 = Seawater scale
% 3 = Free scale
% 4 = NBS scale
%
% (***) Each element must be an integer,
% indicating the K1 K2 dissociation constants that are to be used:
% 1 = Roy, 1993 T: 0-45 S: 5-45. Total scale. Artificial seawater.
% 2 = Goyet & Poisson T: -1-40 S: 10-50. Seaw. scale. Artificial seawater.
% 3 = HANSSON refit BY DICKSON AND MILLERO T: 2-35 S: 20-40. Seaw. scale. Artificial seawater.
% 4 = MEHRBACH refit BY DICKSON AND MILLERO T: 2-35 S: 20-40. Seaw. scale. Artificial seawater.
% 5 = HANSSON and MEHRBACH refit BY DICKSON AND MILLERO T: 2-35 S: 20-40. Seaw. scale. Artificial seawater.
% 6 = GEOSECS (i.e., original Mehrbach) T: 2-35 S: 19-43. NBS scale. Real seawater.
% 7 = Peng (i.e., originam Mehrbach but without XXX) T: 2-35 S: 19-43. NBS scale. Real seawater.
% 8 = Millero, 1979, FOR PURE WATER ONLY (i.e., Sal=0) T: 0-50 S: 0.
% 9 = Cai and Wang, 1998 T: 2-35 S: 0-49. NBS scale. Real and artificial seawater.
% 10 = Lueker et al, 2000 T: 2-35 S: 19-43. Total scale. Real seawater.
% 11 = Mojica Prieto and Millero, 2002. T: 0-45 S: 5-42. Seaw. scale. Real seawater
% 12 = Millero et al, 2002 T: -1.6-35 S: 34-37. Seaw. scale. Field measurements.
% 13 = Millero et al, 2006 T: 0-50 S: 1-50. Seaw. scale. Real seawater.
% 14 = Millero 2010 T: 0-50 S: 1-50. Seaw. scale. Real seawater.
% 15 = Waters, Millero, & Woosley 2014 T: 0-50 S: 1-50. Seaw. scale. Real seawater.
%
% (****) Each element must be an integer that
% indicates the KSO4 dissociation constants that are to be used,
% in combination with the formulation of the borate-to-salinity ratio to be used.
% Having both these choices in a single argument is somewhat awkward,
% but it maintains syntax compatibility with the previous version.
% 1 = KSO4 of Dickson 1990a & TB of Uppstrom 1974 (PREFERRED)
% 2 = KSO4 of Khoo et al 1977 & TB of Uppstrom 1974
% 3 = KSO4 of Dickson 1990a & TB of Lee 2010
% 4 = KSO4 of Khoo et al 1977 & TB of Lee 2010
%
%**************************************************************************%
%
% OUTPUT: * an array containing the following parameter values (one row per sample):
% * a cell-array containing crudely formatted headers
% * a cell-array containing nicely formatted headers
%
% POS PARAMETER UNIT
%
% 01 - TAlk (umol/kgSW)
% 02 - TCO2 (umol/kgSW)
% 03 - pHin ()
% 04 - pCO2 input (uatm)
% 05 - fCO2 input (uatm)
% 06 - HCO3 input (umol/kgSW)
% 07 - CO3 input (umol/kgSW)
% 08 - CO2 input (umol/kgSW)
% 09 - BAlk input (umol/kgSW)
% 10 - OH input (umol/kgSW)
% 11 - PAlk input (umol/kgSW)
% 12 - SiAlk input (umol/kgSW)
% 13 - Hfree input (umol/kgSW)
% 14 - RevelleFactor input ()
% 15 - OmegaCa input ()
% 16 - OmegaAr input ()
% 17 - xCO2 input (ppm)
% 18 - pH output ()
% 19 - pCO2 output (uatm)
% 20 - fCO2 output (uatm)
% 21 - HCO3 output (umol/kgSW)
% 22 - CO3 output (umol/kgSW)
% 23 - CO2 output (umol/kgSW)
% 24 - BAlk output (umol/kgSW)
% 25 - OH output (umol/kgSW)
% 26 - PAlk output (umol/kgSW)
% 27 - SiAlk output (umol/kgSW)
% 28 - Hfree output (umol/kgSW)
% 29 - RevelleFactor output ()
% 30 - OmegaCa output ()
% 31 - OmegaAr output ()
% 32 - xCO2 output (ppm)
% 33 - pH input (Total) ()
% 34 - pH input (SWS) ()
% 35 - pH input (Free) ()
% 36 - pH input (NBS) ()
% 37 - pH output (Total) ()
% 38 - pH output (SWS) ()
% 39 - pH output (Free) ()
% 40 - pH output (NBS) ()
% 41 - TEMP input (deg C) ***
% 42 - TEMPOUT (deg C) ***
% 43 - PRES input (dbar or m) ***
% 44 - PRESOUT (dbar or m) ***
% 45 - PAR1TYPE (integer) ***
% 46 - PAR2TYPE (integer) ***
% 47 - K1K2CONSTANTS (integer) ***
% 48 - KSO4CONSTANTS (integer) ***
% 49 - pHSCALE of input (integer) ***
% 50 - SAL (psu) ***
% 51 - PO4 (umol/kgSW) ***
% 52 - SI (umol/kgSW) ***
% 53 - K0 input ()
% 54 - K1 input ()
% 55 - K2 input ()
% 56 - pK1 input ()
% 57 - pK2 input ()
% 58 - KW input ()
% 59 - KB input ()
% 60 - KF input ()
% 61 - KS input ()
% 62 - KP1 input ()
% 63 - KP2 input ()
% 64 - KP3 input ()
% 65 - KSi input ()
% 66 - K0 output ()
% 67 - K1 output ()
% 68 - K2 output ()
% 69 - pK1 output ()
% 70 - pK2 output ()
% 71 - KW output ()
% 72 - KB output ()
% 73 - KF output ()
% 74 - KS output ()
% 75 - KP1 output ()
% 76 - KP2 output ()
% 77 - KP3 output ()
% 78 - KSi output ()
% 79 - TB (umol/kgSW)
% 80 - TF (umol/kgSW)
% 81 - TS (umol/kgSW)
% 82 - TP (umol/kgSW)
% 83 - TSi (umol/kgSW)
%
% *** SIMPLY RESTATES THE INPUT BY USER
%
% In all the above, the terms "input" and "output" may be understood
% to refer to the 2 scenarios for which CO2SYS performs calculations,
% each defined by its own combination of temperature and pressure.
% For instance, one may use CO2SYS to calculate, from measured DIC and TAlk,
% the pH that that sample will have in the lab (e.g., T=25 degC, P=0 dbar),
% and what the in situ pH would have been (e.g., at T=1 degC, P=4500).
% A = CO2SYS(2400,2200,1,2,35,25,1,0,4200,1,1,1,4,1)
% pH_lab = A(3); % 7.84
% pH_sea = A(18); % 8.05
%
%**************************************************************************
%
% This is version 2.0 (uploaded to CDIAC at SEP XXth, 2011):
%
% **** Changes since 2.0
% - slight changes to allow error propagation
% - new option to choose K1 & K2 from Waters et al. (2014): fixes inconsistencies with Millero (2010) identified by Orr et al. (2015)
%
% **** Changes since 1.01 (uploaded to CDIAC at June 11th, 2009):
% - Function cleans up its global variables when done (if you loose variables, this may be the cause -- see around line 570)
% - Added the outputting of K values
% - Implementation of constants of Cai and Wang, 1998
% - Implementation of constants of Lueker et al., 2000
% - Implementation of constants of Mojica-Prieto and Millero, 2002
% - Implementation of constants of Millero et al., 2002 (only their eqs. 19, 20, no TCO2 dependency)
% - Implementation of constants of Millero et al., 2006
% - Implementation of constants of Millero et al., 2010
% - Properly listed Sal and Temp limits for the available constants
% - added switch for using the new Lee et al., (2010) formulation of Total Borate (see KSO4CONSTANTS above)
% - Minor corrections to the GEOSECS constants (gave NaN for some output in earlier version)
% - Fixed decimal point error on [H+] (did not get converted to umol/kgSW from mol/kgSW).
% - Changed 'Hfreein' to 'Hfreeout' in the 'NICEHEADERS'-output (typo)
%
% **** Changes since 1.00 (uploaded to CDIAC at May 29th, 2009):
% - added a note explaining that all known bugs were removed before release of 1.00
%
%**************************************************************************
%
% CO2SYS originally by Lewis and Wallace 1998
% Converted to MATLAB by Denis Pierrot at
% CIMAS, University of Miami, Miami, Florida
% Vectorization, internal refinements and speed improvements by
% Steven van Heuven, University of Groningen, The Netherlands.
% Questions, bug reports et cetera: svheuven@gmail.com
%
%**************************************************************************
%**************************************************************************
% NOTHING BELOW THIS SHOULD REQUIRE EDITING BY USER!
%**************************************************************************
% Declare global variables
global pHScale WhichKs WhoseKSO4 Pbar
global Sal sqrSal TempK logTempK TempCi TempCo Pdbari Pdbaro;
global FugFac VPFac PengCorrection ntps RGasConstant;
global fH RT;
global K0 K1 K2 KW KB KF KS KP1 KP2 KP3 KSi;
global TB TF TS TP TSi F;
% Added by JM Epitalon
% For computing derivative with respect to Ks, one has to call CO2sys with a perturbed K
% Requested perturbation is passed through the following global variables
global PertK % Id of perturbed K
global Perturb % perturbation
% Input conditioning
% Determine lengths of input vectors
veclengths=[length(PAR1) length(PAR2) length(PAR1TYPE)...
length(PAR2TYPE) length(SAL) length(TEMPIN)...
length(TEMPOUT) length(PRESIN) length(PRESOUT)...
length(SI) length(PO4) length(pHSCALEIN)...
length(K1K2CONSTANTS) length(KSO4CONSTANTS)];
if length(unique(veclengths))>2
disp(' '); disp('*** INPUT ERROR: Input vectors must all be of same length, or of length 1. ***'); disp(' '); return
end
% Make column vectors of all input vectors
PAR1 =PAR1 (:);
PAR2 =PAR2 (:);
PAR1TYPE =PAR1TYPE (:);
PAR2TYPE =PAR2TYPE (:);
SAL =SAL (:);
TEMPIN =TEMPIN (:);
TEMPOUT =TEMPOUT (:);
PRESIN =PRESIN (:);
PRESOUT =PRESOUT (:);
SI =SI (:);
PO4 =PO4 (:);
pHSCALEIN =pHSCALEIN (:);
K1K2CONSTANTS=K1K2CONSTANTS(:);
KSO4CONSTANTS=KSO4CONSTANTS(:);
% Find the longest column vector:
ntps = max(veclengths);
% Populate column vectors
PAR1(1:ntps,1) = PAR1(:) ;
PAR2(1:ntps,1) = PAR2(:) ;
PAR1TYPE(1:ntps,1) = PAR1TYPE(:) ;
PAR2TYPE(1:ntps,1) = PAR2TYPE(:) ;
SAL(1:ntps,1) = SAL(:) ;
TEMPIN(1:ntps,1) = TEMPIN(:) ;
TEMPOUT(1:ntps,1) = TEMPOUT(:) ;
PRESIN(1:ntps,1) = PRESIN(:) ;
PRESOUT(1:ntps,1) = PRESOUT(:) ;
SI(1:ntps,1) = SI(:) ;
PO4(1:ntps,1) = PO4(:) ;
pHSCALEIN(1:ntps,1) = pHSCALEIN(:) ;
K1K2CONSTANTS(1:ntps,1) = K1K2CONSTANTS(:) ;
KSO4CONSTANTS(1:ntps,1) = KSO4CONSTANTS(:) ;
% Assign input to the 'historical' variable names.
pHScale = pHSCALEIN;
WhichKs = K1K2CONSTANTS;
WhoseKSO4 = KSO4CONSTANTS;
p1 = PAR1TYPE;
p2 = PAR2TYPE;
TempCi = TEMPIN;
TempCo = TEMPOUT;
Pdbari = PRESIN;
Pdbaro = PRESOUT;
Sal = SAL;
sqrSal = sqrt(SAL);
TP = PO4;
TSi = SI;
RGasConstant = 83.1451; % ml bar-1 K-1 mol-1, DOEv2
%RGasConstant = 83.14472; % ml bar-1 K-1 mol-1, DOEv3
% Generate empty vectors for...
TA = nan(ntps,1); % Talk
TC = nan(ntps,1); % DIC
PH = nan(ntps,1); % pH
PC = nan(ntps,1); % pCO2
FC = nan(ntps,1); % fCO2
% Assign values to empty vectors.
F=(p1==1); TA(F)=PAR1(F)/1e6; % Convert from micromol/kg to mol/kg
F=(p1==2); TC(F)=PAR1(F)/1e6; % Convert from micromol/kg to mol/kg
F=(p1==3); PH(F)=PAR1(F);
F=(p1==4); PC(F)=PAR1(F)/1e6; % Convert from microatm. to atm.
F=(p1==5); FC(F)=PAR1(F)/1e6; % Convert from microatm. to atm.
F=(p2==1); TA(F)=PAR2(F)/1e6; % Convert from micromol/kg to mol/kg
F=(p2==2); TC(F)=PAR2(F)/1e6; % Convert from micromol/kg to mol/kg
F=(p2==3); PH(F)=PAR2(F);
F=(p2==4); PC(F)=PAR2(F)/1e6; % Convert from microatm. to atm.
F=(p2==5); FC(F)=PAR2(F)/1e6; % Convert from microatm. to atm.
% Generate the columns holding Si, Phos and Sal.
% Pure Water case:
F=(WhichKs==8);
Sal(F) = 0;
% GEOSECS and Pure Water:
F=(WhichKs==8 | WhichKs==6);
TP(F) = 0;
TSi(F) = 0;
% All other cases
F=~F;
TP(F) = TP(F)./1e6;
TSi(F) = TSi(F)./1e6;
% The vector 'PengCorrection' is used to modify the value of TA, for those
% cases where WhichKs==7, since PAlk(Peng) = PAlk(Dickson) + TP.
% Thus, PengCorrection is 0 for all cases where WhichKs is not 7
PengCorrection=zeros(ntps,1); F=WhichKs==7; PengCorrection(F)=TP(F);
% Calculate the constants for all samples at input conditions
% The constants calculated for each sample will be on the appropriate pH scale!
Constants(TempCi,Pdbari);
% Added by JM Epitalon
% For computing derivative with respect to Ks, one has to perturb the value of one K
% Requested perturbation is passed through global variables PertK and Perturb
if (~ isempty(PertK))
switch PertK
case {'K0'}
K0 = K0 + Perturb;
case {'K1'}
K1 = K1 + Perturb;
case {'K2'}
K2 = K2 + Perturb;
case {'KB'}
KB = KB + Perturb;
case {'KW'}
KW = KW + Perturb;
case {'BOR'}
TB = TB + Perturb;
end
end
% Make sure fCO2 is available for each sample that has pCO2.
F=find(p1==4 | p2==4); FC(F) = PC(F).*FugFac(F);
% Generate vector for results, and copy the raw input values into them. This
% copies ~60% NaNs, which will be replaced for calculated values later on.
TAc = TA;
TCc = TC;
PHic = PH;
PCic = PC;
FCic = FC;
% Generate vector describing the combination of input parameters
% So, the valid ones are: 12,13,15,23,25,35
Icase = 10*min(p1,p2) + max(p1,p2);
% Calculate missing values for AT,CT,PH,FC:
% pCO2 will be calculated later on, routines work with fCO2.
F=Icase==12; % input TA, TC
if any(F)
[PHic(F) FCic(F)] = CalculatepHfCO2fromTATC(TAc(F)-PengCorrection(F), TCc(F));
end
F=Icase==13; % input TA, pH
if any(F)
TCc(F) = CalculateTCfromTApH(TAc(F)-PengCorrection(F), PHic(F));
FCic(F) = CalculatefCO2fromTCpH(TCc(F), PHic(F));
end
F=Icase==14 | Icase==15; % input TA, (pCO2 or fCO2)
if any(F)
PHic(F) = CalculatepHfromTAfCO2(TAc(F)-PengCorrection(F), FCic(F));
TCc(F) = CalculateTCfromTApH (TAc(F)-PengCorrection(F), PHic(F));
end
F=Icase==23; % input TC, pH
if any(F)
TAc(F) = CalculateTAfromTCpH (TCc(F), PHic(F)) + PengCorrection(F);
FCic(F) = CalculatefCO2fromTCpH(TCc(F), PHic(F));
end
F=Icase==24 | Icase==25; % input TC, (pCO2 or fCO2)
if any(F)
PHic(F) = CalculatepHfromTCfCO2(TCc(F), FCic(F));
TAc(F) = CalculateTAfromTCpH (TCc(F), PHic(F)) + PengCorrection(F);
end
F=Icase==34 | Icase==35; % input pH, (pCO2 or fCO2)
if any(F)
TCc(F) = CalculateTCfrompHfCO2(PHic(F), FCic(F));
TAc(F) = CalculateTAfromTCpH (TCc(F), PHic(F)) + PengCorrection(F);
end
% By now, an fCO2 value is available for each sample.
% Generate the associated pCO2 value:
PCic = FCic./FugFac;
% CalculateOtherParamsAtInputConditions:
[HCO3inp CO3inp BAlkinp...
OHinp PAlkinp...
SiAlkinp Hfreeinp...
HSO4inp HFinp] = CalculateAlkParts(PHic, TCc);
PAlkinp = PAlkinp+PengCorrection;
CO2inp = TCc - CO3inp - HCO3inp;
F=true(ntps,1); % i.e., do for all samples:
Revelleinp = RevelleFactor(TAc-PengCorrection, TCc);
[OmegaCainp OmegaArinp] = CaSolubility(Sal, TempCi, Pdbari, TCc, PHic);
xCO2dryinp = PCic./VPFac; % ' this assumes pTot = 1 atm
% Just for reference, convert pH at input conditions to the other scales, too.
[pHicT pHicS pHicF pHicN]=FindpHOnAllScales(PHic);
% Merge the Ks at input into an array. Ks at output will be glued to this later.
KIVEC=[K0 K1 K2 -log10(K1) -log10(K2) KW KB KF KS KP1 KP2 KP3 KSi];
% Calculate the constants for all samples at output conditions
Constants(TempCo,Pdbaro);
% Added by JM Epitalon
% For computing derivative with respect to Ks, one has to perturb the value of one K
% Requested perturbation is passed through global variables PertK and Perturb
if (~ isempty(PertK))
switch PertK
case {'K0'}
K0 = K0 + Perturb;
case {'K1'}
K1 = K1 + Perturb;
case {'K2'}
K2 = K2 + Perturb;
case {'KB'}
KB = KB + Perturb;
case {'KW'}
KW = KW + Perturb;
case {'KW'}
TB = TB + Perturb;
end
end
% Calculate, for output conditions, using conservative TA and TC, pH, fCO2 and pCO2
F=true(ntps,1); % i.e., do for all samples:
[PHoc FCoc] = CalculatepHfCO2fromTATC(TAc-PengCorrection, TCc);
PCoc = FCoc./FugFac;
% Calculate Other Stuff At Output Conditions:
[HCO3out CO3out BAlkout...
OHout PAlkout...
SiAlkout Hfreeout...
HSO4out HFout] = CalculateAlkParts(PHoc, TCc);
PAlkout = PAlkout+PengCorrection;
CO2out = TCc - CO3out - HCO3out;
Revelleout = RevelleFactor(TAc, TCc);
[OmegaCaout OmegaArout] = CaSolubility(Sal, TempCo, Pdbaro, TCc, PHoc);
xCO2dryout = PCoc./VPFac; % ' this assumes pTot = 1 atm
% Just for reference, convert pH at output conditions to the other scales, too.
[pHocT pHocS pHocF pHocN]=FindpHOnAllScales(PHoc);
KOVEC=[K0 K1 K2 -log10(K1) -log10(K2) KW KB KF KS KP1 KP2 KP3 KSi];
TVEC =[TB TF TS];
% Saving data in array, 81 columns, as many rows as samples input
DATA=[TAc*1e6 TCc*1e6 PHic PCic*1e6 FCic*1e6...
HCO3inp*1e6 CO3inp*1e6 CO2inp*1e6 BAlkinp*1e6 OHinp*1e6...
PAlkinp*1e6 SiAlkinp*1e6 Hfreeinp*1e6 Revelleinp OmegaCainp... %%% Multiplied Hfreeinp *1e6, svh20100827
OmegaArinp xCO2dryinp*1e6 PHoc PCoc*1e6 FCoc*1e6...
HCO3out*1e6 CO3out*1e6 CO2out*1e6 BAlkout*1e6 OHout*1e6...
PAlkout*1e6 SiAlkout*1e6 Hfreeout*1e6 Revelleout OmegaCaout... %%% Multiplied Hfreeout *1e6, svh20100827
OmegaArout xCO2dryout*1e6 pHicT pHicS pHicF...
pHicN pHocT pHocS pHocF pHocN...
TEMPIN TEMPOUT PRESIN PRESOUT PAR1TYPE...
PAR2TYPE K1K2CONSTANTS KSO4CONSTANTS pHSCALEIN SAL...
PO4 SI KIVEC KOVEC TVEC*1e6];
HEADERS={'TAlk';'TCO2';'pHin';'pCO2in';'fCO2in';'HCO3in';'CO3in';...
'CO2in';'BAlkin';'OHin';'PAlkin';'SiAlkin';'Hfreein';'RFin';...
'OmegaCAin';'OmegaARin';'xCO2in';'pHout';'pCO2out';'fCO2out';...
'HCO3out';'CO3out';'CO2out';'BAlkout';'OHout';'PAlkout';...
'SiAlkout';'Hfreeout';'RFout';'OmegaCAout';'OmegaARout';'xCO2out';
'pHinTOTAL';'pHinSWS';'pHinFREE';'pHinNBS';...
'pHoutTOTAL';'pHoutSWS';'pHoutFREE';'pHoutNBS';'TEMPIN';'TEMPOUT';...
'PRESIN';'PRESOUT';'PAR1TYPE';'PAR2TYPE';'K1K2CONSTANTS';'KSO4CONSTANTS';...
'pHSCALEIN';'SAL';'PO4';'SI';'K0input';'K1input';'K2input';'pK1input';...
'pK2input';'KWinput';'KBinput';'KFinput';'KSinput';'KP1input';'KP2input';...
'KP3input';'KSiinput';'K0output';'K1output';'K2output';'pK1output';...
'pK2output';'KWoutput';'KBoutput';'KFoutput';'KSoutput';'KP1output';...
'KP2output';'KP3output';'KSioutput';'TB';'TF';'TS';};
NICEHEADERS={...
'01 - TAlk (umol/kgSW) ';
'02 - TCO2 (umol/kgSW) ';
'03 - pHin () ';
'04 - pCO2in (uatm) ';
'05 - fCO2in (uatm) ';
'06 - HCO3in (umol/kgSW) ';
'07 - CO3in (umol/kgSW) ';
'08 - CO2in (umol/kgSW) ';
'09 - BAlkin (umol/kgSW) ';
'10 - OHin (umol/kgSW) ';
'11 - PAlkin (umol/kgSW) ';
'12 - SiAlkin (umol/kgSW) ';
'13 - Hfreein (umol/kgSW) ';
'14 - RevelleFactorin () ';
'15 - OmegaCain () ';
'16 - OmegaArin () ';
'17 - xCO2in (ppm) ';
'18 - pHout () ';
'19 - pCO2out (uatm) ';
'20 - fCO2out (uatm) ';
'21 - HCO3out (umol/kgSW) ';
'22 - CO3out (umol/kgSW) ';
'23 - CO2out (umol/kgSW) ';
'24 - BAlkout (umol/kgSW) ';
'25 - OHout (umol/kgSW) ';
'26 - PAlkout (umol/kgSW) ';
'27 - SiAlkout (umol/kgSW) ';
'28 - Hfreeout (umol/kgSW) '; %%% Changed 'Hfreein' to 'Hfreeout', svh20100827
'29 - RevelleFactorout () ';
'30 - OmegaCaout () ';
'31 - OmegaArout () ';
'32 - xCO2out (ppm) ';
'33 - pHin (Total) () ';
'34 - pHin (SWS) () ';
'35 - pHin (Free) () ';
'36 - pHin (NBS ) () ';
'37 - pHout(Total) () ';
'38 - pHout(SWS) () ';
'39 - pHout(Free) () ';
'40 - pHout(NBS ) () ';
'41 - TEMPIN (Deg C) ';
'42 - TEMPOUT (Deg C) ';
'43 - PRESIN (dbar) ';
'44 - PRESOUT (dbar) ';
'45 - PAR1TYPE () ';
'46 - PAR2TYPE () ';
'47 - K1K2CONSTANTS () ';
'48 - KSO4CONSTANTS () ';
'49 - pHSCALEIN () ';
'50 - SAL (umol/kgSW) ';
'51 - PO4 (umol/kgSW) ';
'52 - SI (umol/kgSW) ';
'53 - K0input () ';
'54 - K1input () ';
'55 - K2input () ';
'56 - pK1input () ';
'57 - pK2input () ';
'58 - KWinput () ';
'59 - KBinput () ';
'60 - KFinput () ';
'61 - KSinput () ';
'62 - KP1input () ';
'63 - KP2input () ';
'64 - KP3input () ';
'65 - KSiinput () ';
'66 - K0output () ';
'67 - K1output () ';
'68 - K2output () ';
'69 - pK1output () ';
'70 - pK2output () ';
'71 - KWoutput () ';
'72 - KBoutput () ';
'73 - KFoutput () ';
'74 - KSoutput () ';
'75 - KP1output () ';
'76 - KP2output () ';
'77 - KP3output () ';
'78 - KSioutput () ';
'79 - TB (umol/kgSW) ';
'80 - TF (umol/kgSW) ';
'81 - TS (umol/kgSW) '};
clear global F K2 KP3 Pdbari Sal TS VPFac ntps
clear global FugFac KB KS Pdbaro T TSi WhichKs pHScale
clear global K KF KSi PengCorrection TB TempCi WhoseKSO4 sqrSal
clear global K0 KP1 KW RGasConstant TF TempCo fH
clear global K1 KP2 Pbar RT TP TempK logTempK
end % end main function
%**************************************************************************
% Subroutines:
%**************************************************************************
function Constants(TempC,Pdbar)
global pHScale WhichKs WhoseKSO4 sqrSal Pbar RT;
global K0 fH FugFac VPFac ntps TempK logTempK;
global K1 K2 KW KB KF KS KP1 KP2 KP3 KSi;
global TB TF TS TP TSi RGasConstant Sal;
% SUB Constants, version 04.01, 10-13-97, written by Ernie Lewis.
% Inputs: pHScale%, WhichKs%, WhoseKSO4%, Sali, TempCi, Pdbar
% Outputs: K0, K(), T(), fH, FugFac, VPFac
% This finds the Constants of the CO2 system in seawater or freshwater,
% corrects them for pressure, and reports them on the chosen pH scale.
% The process is as follows: the Constants (except KS, KF which stay on the
% free scale - these are only corrected for pressure) are
% 1) evaluated as they are given in the literature
% 2) converted to the SWS scale in mol/kg-SW or to the NBS scale
% 3) corrected for pressure
% 4) converted to the SWS pH scale in mol/kg-SW
% 5) converted to the chosen pH scale
%
% PROGRAMMER'S NOTE: all logs are log base e
% PROGRAMMER'S NOTE: all Constants are converted to the pH scale
% pHScale% (the chosen one) in units of mol/kg-SW
% except KS and KF are on the free scale
% and KW is in units of (mol/kg-SW)^2
TempK = TempC + 273.15;
RT = RGasConstant.*TempK;
logTempK = log(TempK);
Pbar = Pdbar ./ 10;
% Generate empty vectors for holding results
TB = nan(ntps,1);
TF = nan(ntps,1);
TS = nan(ntps,1);
% CalculateTB - Total Borate:
F=(WhichKs==8); % Pure water case.
if any(F)
TB(F) = 0;
end
F=(WhichKs==6 | WhichKs==7);
if any(F)
TB(F) = 0.0004106.*Sal(F)./35; % in mol/kg-SW
% this is .00001173.*Sali
% this is about 1% lower than Uppstrom's value
% Culkin, F., in Chemical Oceanography,
% ed. Riley and Skirrow, 1965:
% GEOSECS references this, but this value is not explicitly
% given here
end
F=(WhichKs~=6 & WhichKs~=7 & WhichKs~=8); % All other cases
if any(F)
FF=F&(WhoseKSO4==1|WhoseKSO4==2); % If user opted for Uppstrom's values:
if any(FF)
% Uppstrom, L., Deep-Sea Research 21:161-162, 1974:
% this is .000416.*Sali./35. = .0000119.*Sali
% TB(FF) = (0.000232./10.811).*(Sal(FF)./1.80655); % in mol/kg-SW
TB(FF) = 0.0004157.*Sal(FF)./35; % in mol/kg-SW
end
FF=F&(WhoseKSO4==3|WhoseKSO4==4); % If user opted for the new Lee values:
if any(FF)
% Lee, Kim, Byrne, Millero, Feely, Yong-Ming Liu. 2010.
% Geochimica Et Cosmochimica Acta 74 (6): 1801–1811.
TB(FF) = 0.0004326.*Sal(FF)./35; % in mol/kg-SW
end
end
% CalculateTF;
% Riley, J. P., Deep-Sea Research 12:219-220, 1965:
% this is .000068.*Sali./35. = .00000195.*Sali
TF = (0.000067./18.998).*(Sal./1.80655); % in mol/kg-SW
% CalculateTS ;
% Morris, A. W., and Riley, J. P., Deep-Sea Research 13:699-705, 1966:
% this is .02824.*Sali./35. = .0008067.*Sali
TS = (0.14./96.062).*(Sal./1.80655); % in mol/kg-SW
% CalculateK0:
% Weiss, R. F., Marine Chemistry 2:203-215, 1974.
TempK100 = TempK./100;
lnK0 = -60.2409 + 93.4517 ./ TempK100 + 23.3585 .* log(TempK100) + Sal .*...
(0.023517 - 0.023656 .* TempK100 + 0.0047036 .* TempK100 .^2);
K0 = exp(lnK0); % this is in mol/kg-SW/atm
% CalculateIonS:
% This is from the DOE handbook, Chapter 5, p. 13/22, eq. 7.2.4:
IonS = 19.924 .* Sal ./ (1000 - 1.005 .* Sal);
% CalculateKS:
lnKS = nan(ntps,1); pKS = nan(ntps,1); KS = nan(ntps,1);
F=(WhoseKSO4==1|WhoseKSO4==3);
if any(F)
% Dickson, A. G., J. Chemical Thermodynamics, 22:113-127, 1990
% The goodness of fit is .021.
% It was given in mol/kg-H2O. I convert it to mol/kg-SW.
% TYPO on p. 121: the constant e9 should be e8.
% This is from eqs 22 and 23 on p. 123, and Table 4 on p 121:
lnKS(F) = -4276.1./TempK(F) + 141.328 - 23.093.*logTempK(F) +...
(-13856./TempK(F) + 324.57 - 47.986.*logTempK(F)).*sqrt(IonS(F)) +...
(35474./TempK(F) - 771.54 + 114.723.*logTempK(F)).*IonS(F) +...
(-2698./TempK(F)).*sqrt(IonS(F)).*IonS(F) + (1776./TempK(F)).*IonS(F).^2;
KS(F) = exp(lnKS(F))... % this is on the free pH scale in mol/kg-H2O
.* (1 - 0.001005 .* Sal(F)); % convert to mol/kg-SW
end
F=(WhoseKSO4==2|WhoseKSO4==4);
if any(F)
% Khoo, et al, Analytical Chemistry, 49(1):29-34, 1977
% KS was found by titrations with a hydrogen electrode
% of artificial seawater containing sulfate (but without F)
% at 3 salinities from 20 to 45 and artificial seawater NOT
% containing sulfate (nor F) at 16 salinities from 15 to 45,
% both at temperatures from 5 to 40 deg C.
% KS is on the Free pH scale (inherently so).
% It was given in mol/kg-H2O. I convert it to mol/kg-SW.
% He finds log(beta) which = my pKS;
% his beta is an association constant.
% The rms error is .0021 in pKS, or about .5% in KS.
% This is equation 20 on p. 33:
pKS(F) = 647.59 ./ TempK(F) - 6.3451 + 0.019085.*TempK(F) - 0.5208.*sqrt(IonS(F));
KS(F) = 10.^(-pKS(F))... % this is on the free pH scale in mol/kg-H2O
.* (1 - 0.001005.*Sal(F)); % convert to mol/kg-SW
end
% CalculateKF:
% Dickson, A. G. and Riley, J. P., Marine Chemistry 7:89-99, 1979:
lnKF = 1590.2./TempK - 12.641 + 1.525.*IonS.^0.5;
KF = exp(lnKF)... % this is on the free pH scale in mol/kg-H2O
.*(1 - 0.001005.*Sal); % convert to mol/kg-SW
% Another expression exists for KF: Perez and Fraga 1987. Not used here since ill defined for low salinity. (to be used for S: 10-40, T: 9-33)
% Nonetheless, P&F87 might actually be better than the fit of D&R79 above, which is based on only three salinities: [0 26.7 34.6]
% lnKF = 874./TempK - 9.68 + 0.111.*Sal.^0.5;
% KF = exp(lnKF); % this is on the free pH scale in mol/kg-SW
% CalculatepHScaleConversionFactors:
% These are NOT pressure-corrected.
SWStoTOT = (1 + TS./KS)./(1 + TS./KS + TF./KF);
FREEtoTOT = 1 + TS./KS;
% CalculatefH
fH = nan(ntps,1);
% Use GEOSECS's value for cases 1,2,3,4,5 (and 6) to convert pH scales.
F=(WhichKs==8);
if any(F)
fH(F) = 1; % this shouldn't occur in the program for this case
end
F=(WhichKs==7);
if any(F)
fH(F) = 1.29 - 0.00204.* TempK(F) + (0.00046 -...
0.00000148.*TempK(F)).*Sal(F).*Sal(F);
% Peng et al, Tellus 39B:439-458, 1987:
% They reference the GEOSECS report, but round the value
% given there off so that it is about .008 (1%) lower. It
% doesn't agree with the check value they give on p. 456.
end
F=(WhichKs~=7 & WhichKs~=8);
if any(F)
fH(F) = 1.2948 - 0.002036.*TempK(F) + (0.0004607 -...
0.000001475.*TempK(F)).*Sal(F).^2;
% Takahashi et al, Chapter 3 in GEOSECS Pacific Expedition,
% v. 3, 1982 (p. 80);
end
% CalculateKB:
KB = nan(ntps,1); logKB = nan(ntps,1);
lnKBtop = nan(ntps,1); lnKB = nan(ntps,1);
F=(WhichKs==8); % Pure water case
if any(F)
KB(F) = 0;
end
F=(WhichKs==6 | WhichKs==7);
if any(F)
% This is for GEOSECS and Peng et al.
% Lyman, John, UCLA Thesis, 1957
% fit by Li et al, JGR 74:5507-5525, 1969:
logKB(F) = -9.26 + 0.00886.*Sal(F) + 0.01.*TempC(F);
KB(F) = 10.^(logKB(F))... % this is on the NBS scale
./fH(F); % convert to the SWS scale
end
F=(WhichKs~=6 & WhichKs~=7 & WhichKs~=8);
if any(F)
% Dickson, A. G., Deep-Sea Research 37:755-766, 1990:
lnKBtop(F) = -8966.9 - 2890.53.*sqrSal(F) - 77.942.*Sal(F) +...
1.728.*sqrSal(F).*Sal(F) - 0.0996.*Sal(F).^2;
lnKB(F) = lnKBtop(F)./TempK(F) + 148.0248 + 137.1942.*sqrSal(F) +...
1.62142.*Sal(F) + (-24.4344 - 25.085.*sqrSal(F) - 0.2474.*...
Sal(F)).*logTempK(F) + 0.053105.*sqrSal(F).*TempK(F);
KB(F) = exp(lnKB(F))... % this is on the total pH scale in mol/kg-SW
./SWStoTOT(F); % convert to SWS pH scale
end
% CalculateKW:
lnKW = nan(ntps,1); KW = nan(ntps,1);
F=(WhichKs==7);
if any(F)
% Millero, Geochemica et Cosmochemica Acta 43:1651-1661, 1979
lnKW(F) = 148.9802 - 13847.26./TempK(F) - 23.6521.*logTempK(F) +...
(-79.2447 + 3298.72./TempK(F) + 12.0408.*logTempK(F)).*...
sqrSal(F) - 0.019813.*Sal(F);
end
F=(WhichKs==8);
if any(F)
% Millero, Geochemica et Cosmochemica Acta 43:1651-1661, 1979
% refit data of Harned and Owen, The Physical Chemistry of
% Electrolyte Solutions, 1958
lnKW(F) = 148.9802 - 13847.26./TempK(F) - 23.6521.*logTempK(F);
end
F=(WhichKs~=6 & WhichKs~=7 & WhichKs~=8);
if any(F)
% Millero, Geochemica et Cosmochemica Acta 59:661-677, 1995.
% his check value of 1.6 umol/kg-SW should be 6.2
lnKW(F) = 148.9802 - 13847.26./TempK(F) - 23.6521.*logTempK(F) +...
(-5.977 + 118.67./TempK(F) + 1.0495.*logTempK(F)).*...
sqrSal(F) - 0.01615.*Sal(F);
end
KW = exp(lnKW); % this is on the SWS pH scale in (mol/kg-SW)^2
F=(WhichKs==6);
if any(F)
KW(F) = 0; % GEOSECS doesn't include OH effects
end
% CalculateKP1KP2KP3KSi:
KP1 = nan(ntps,1); KP2 = nan(ntps,1);
KP3 = nan(ntps,1); KSi = nan(ntps,1);
lnKP1 = nan(ntps,1); lnKP2 = nan(ntps,1);
lnKP3 = nan(ntps,1); lnKSi = nan(ntps,1);
F=(WhichKs==7);
if any(F)
KP1(F) = 0.02;
% Peng et al don't include the contribution from this term,
% but it is so small it doesn't contribute. It needs to be
% kept so that the routines work ok.
% KP2, KP3 from Kester, D. R., and Pytkowicz, R. M.,
% Limnology and Oceanography 12:243-252, 1967:
% these are only for sals 33 to 36 and are on the NBS scale
KP2(F) = exp(-9.039 - 1450./TempK(F))... % this is on the NBS scale
./fH(F); % convert to SWS scale
KP3(F) = exp(4.466 - 7276./TempK(F))... % this is on the NBS scale
./fH(F); % convert to SWS scale
% Sillen, Martell, and Bjerrum, Stability Constants of metal-ion complexes,
% The Chemical Society (London), Special Publ. 17:751, 1964:
KSi(F) = 0.0000000004... % this is on the NBS scale
./fH(F); % convert to SWS scale
end
F=(WhichKs==6 | WhichKs==8);
if any(F)
KP1(F) = 0; KP2(F) = 0; KP3(F) = 0; KSi(F) = 0;
% Neither the GEOSECS choice nor the freshwater choice
% include contributions from phosphate or silicate.
end
F=(WhichKs~=6 & WhichKs~=7 & WhichKs~=8);
if any(F)
% Yao and Millero, Aquatic Geochemistry 1:53-88, 1995
% KP1, KP2, KP3 are on the SWS pH scale in mol/kg-SW.
% KSi was given on the SWS pH scale in molal units.
lnKP1(F) = -4576.752./TempK(F) + 115.54 - 18.453.*logTempK(F) + (-106.736./TempK(F) +...
0.69171).*sqrSal(F) + (-0.65643./TempK(F) - 0.01844).*Sal(F);
KP1(F) = exp(lnKP1(F));
lnKP2(F) = -8814.715./TempK(F) + 172.1033 - 27.927.*logTempK(F) + (-160.34./TempK(F) +...
1.3566).*sqrSal(F) + (0.37335./TempK(F) - 0.05778).*Sal(F);
KP2(F) = exp(lnKP2(F));
lnKP3(F) = -3070.75./TempK(F) - 18.126 + (17.27039./TempK(F) + 2.81197).*sqrSal(F) +...
(-44.99486./TempK(F) - 0.09984).*Sal(F);
KP3(F) = exp(lnKP3(F));
lnKSi(F) = -8904.2./TempK(F) + 117.4 - 19.334.*logTempK(F) + (-458.79./TempK(F) +...
3.5913).*sqrt(IonS(F)) + (188.74./TempK(F) - 1.5998).*IonS(F) +...
(-12.1652./TempK(F) + 0.07871).*IonS(F).^2;
KSi(F) = exp(lnKSi(F))... % this is on the SWS pH scale in mol/kg-H2O
.*(1 - 0.001005.*Sal(F)); % convert to mol/kg-SW
end
% CalculateK1K2:
logK1 = nan(ntps,1); lnK1 = nan(ntps,1);
pK1 = nan(ntps,1); K1 = nan(ntps,1);
logK2 = nan(ntps,1); lnK2 = nan(ntps,1);
pK2 = nan(ntps,1); K2 = nan(ntps,1);
F=(WhichKs==1);
if any(F)
% ROY et al, Marine Chemistry, 44:249-267, 1993
% (see also: Erratum, Marine Chemistry 45:337, 1994
% and Erratum, Marine Chemistry 52:183, 1996)
% Typo: in the abstract on p. 249: in the eq. for lnK1* the
% last term should have S raised to the power 1.5.
% They claim standard deviations (p. 254) of the fits as
% .0048 for lnK1 (.5% in K1) and .007 in lnK2 (.7% in K2).
% They also claim (p. 258) 2s precisions of .004 in pK1 and
% .006 in pK2. These are consistent, but Andrew Dickson
% (personal communication) obtained an rms deviation of about
% .004 in pK1 and .003 in pK2. This would be a 2s precision
% of about 2% in K1 and 1.5% in K2.
% T: 0-45 S: 5-45. Total Scale. Artificial sewater.
% This is eq. 29 on p. 254 and what they use in their abstract:
lnK1(F) = 2.83655 - 2307.1266./TempK(F) - 1.5529413.*logTempK(F) +...
(-0.20760841 - 4.0484./TempK(F)).*sqrSal(F) + 0.08468345.*Sal(F) -...
0.00654208.*sqrSal(F).*Sal(F);
K1(F) = exp(lnK1(F))... % this is on the total pH scale in mol/kg-H2O
.*(1 - 0.001005.*Sal(F))... % convert to mol/kg-SW
./SWStoTOT(F); % convert to SWS pH scale
% This is eq. 30 on p. 254 and what they use in their abstract:
lnK2(F) = -9.226508 - 3351.6106./TempK(F) - 0.2005743.*logTempK(F) +...
(-0.106901773 - 23.9722./TempK(F)).*sqrSal(F) + 0.1130822.*Sal(F) -...
0.00846934.*sqrSal(F).*Sal(F);
K2(F) = exp(lnK2(F))... % this is on the total pH scale in mol/kg-H2O
.*(1 - 0.001005.*Sal(F))... % convert to mol/kg-SW
./SWStoTOT(F); % convert to SWS pH scale
end
F=(WhichKs==2);
if any(F)
% GOYET AND POISSON, Deep-Sea Research, 36(11):1635-1654, 1989
% The 2s precision in pK1 is .011, or 2.5% in K1.
% The 2s precision in pK2 is .02, or 4.5% in K2.
% This is in Table 5 on p. 1652 and what they use in the abstract:
pK1(F) = 812.27./TempK(F) + 3.356 - 0.00171.*Sal(F).*logTempK(F)...
+ 0.000091.*Sal(F).^2;
K1(F) = 10.^(-pK1(F)); % this is on the SWS pH scale in mol/kg-SW
% This is in Table 5 on p. 1652 and what they use in the abstract:
pK2(F) = 1450.87./TempK(F) + 4.604 - 0.00385.*Sal(F).*logTempK(F)...
+ 0.000182.*Sal(F).^2;
K2(F) = 10.^(-pK2(F)); % this is on the SWS pH scale in mol/kg-SW
end
F=(WhichKs==3);
if any(F)
% HANSSON refit BY DICKSON AND MILLERO
% Dickson and Millero, Deep-Sea Research, 34(10):1733-1743, 1987
% (see also Corrigenda, Deep-Sea Research, 36:983, 1989)
% refit data of Hansson, Deep-Sea Research, 20:461-478, 1973
% and Hansson, Acta Chemica Scandanavia, 27:931-944, 1973.
% on the SWS pH scale in mol/kg-SW.
% Hansson gave his results on the Total scale (he called it
% the seawater scale) and in mol/kg-SW.
% Typo in DM on p. 1739 in Table 4: the equation for pK2*
% for Hansson should have a .000132 *S^2
% instead of a .000116 *S^2.
% The 2s precision in pK1 is .013, or 3% in K1.
% The 2s precision in pK2 is .017, or 4.1% in K2.
% This is from Table 4 on p. 1739.
pK1(F) = 851.4./TempK(F) + 3.237 - 0.0106.*Sal(F) + 0.000105.*Sal(F).^2;
K1(F) = 10.^(-pK1(F)); % this is on the SWS pH scale in mol/kg-SW
% This is from Table 4 on p. 1739.
pK2(F) = -3885.4./TempK(F) + 125.844 - 18.141.*logTempK(F)...
- 0.0192.*Sal(F) + 0.000132.*Sal(F).^2;
K2(F) = 10.^(-pK2(F)); % this is on the SWS pH scale in mol/kg-SW
end
F=(WhichKs==4);
if any(F)
% MEHRBACH refit BY DICKSON AND MILLERO
% Dickson and Millero, Deep-Sea Research, 34(10):1733-1743, 1987
% (see also Corrigenda, Deep-Sea Research, 36:983, 1989)
% refit data of Mehrbach et al, Limn Oc, 18(6):897-907, 1973
% on the SWS pH scale in mol/kg-SW.
% Mehrbach et al gave results on the NBS scale.
% The 2s precision in pK1 is .011, or 2.6% in K1.
% The 2s precision in pK2 is .020, or 4.6% in K2.
% Valid for salinity 20-40.
% This is in Table 4 on p. 1739.
pK1(F) = 3670.7./TempK(F) - 62.008 + 9.7944.*logTempK(F)...
- 0.0118.*Sal(F) + 0.000116.*Sal(F).^2;
K1(F) = 10.^(-pK1(F)); % this is on the SWS pH scale in mol/kg-SW
% This is in Table 4 on p. 1739.
pK2(F) = 1394.7./TempK(F) + 4.777 - 0.0184.*Sal(F) + 0.000118.*Sal(F).^2;
K2(F) = 10.^(-pK2(F)); % this is on the SWS pH scale in mol/kg-SW
end