CO2 Headspace equilibration.jsl

//
// THIS SCRIPT WAS WRITTEN BY YVES PRAIRIE (UNIVERSITÉ DU QUEBEC À MONTRÉAL) JULY 2020. 
// FOR EASE OF USE,THE SAME CODE IS ALSO AVAILABLE AS A JMP ADD-IN
//  
//

Names Default To Here( 1 );
Clear Globals();

dt = Current Data Table();
filename = dt << get name;
If( Is Empty( dt ),
	Try( dt = Open(), Throw( "No data table found" ) ),
	filename = dt << get name
);

dlg = New Window( "CO2 equilibration calculator",
	<<modal,
	H List Box(
		Panel Box( "Select Columns", MDSColList = Col List Box( All, width( 140 ), nlines( 10 ) ) ),
		V List Box(
			Panel Box( "Required variables",
				Lineup Box( N Col( 2 ), Spacing( 3 ), 
					
					bbc = Button Box( "HS CO2 before equilibration (ppmv)", CGPCO2 << Append( MDSColList << GetSelected ) ),
					CGPCO2 = Col List Box( width( 100 ), nLines( 1 ), MinItems( 1 ) ),
					Button Box( "HS CO2 after equilibration (ppmv)", CPCO2 << Append( MDSColList << GetSelected ) ),
					CPCO2 = Col List Box( width( 100 ), nLines( 1 ), MinItems( 1 ) ),
					Mbb = Button Box( "Water Temperature (C) during sampling", CTw << Append( MDSColList << GetSelected ) ),
					CTw = Col List Box( width( 100 ), nLines( 1 ), MinItems( 1 ) ),
					Mbb = Button Box( "Water Temperature (C) during equil. ", CT << Append( MDSColList << GetSelected ) ),
					CT = Col List Box( width( 100 ), nLines( 1 ), MinItems( 1 ) ),
					Button Box( "Alkalinity (uequiv/L)", Calk << Append( MDSColList << GetSelected ) ),
					CAlk = Col List Box( width( 100 ), nLines( 1 ), MinItems( 1 ) ),
					bbc2 = Button Box( "Headspace ratio (gas-to-liquid volume ratio)", CHSR << Append( MDSColList << GetSelected ) ),
					H List Box( CHSR = Col List Box( width( 140 ), nLines( 1 ), MinItems( 1 ) ) ), 
					bbc3 = Button Box( "Barometric Pressure (kPa)", CBP << Append( MDSColList << GetSelected ) ),
					H List Box( CBP = Col List Box( width( 140 ), nLines( 1 ), MinItems( 1 ) ) ),
					bbc4 = Button Box( "Salinity (PSU)", CSAL << Append( MDSColList << GetSelected ) ),
					H List Box( CSAL = Col List Box( width( 140 ), nLines( 1 ), MinItems( 1 ) ) ),bbc4<<Enable(0),CSAL<<Enable(0),
				),
				H List Box(
					cb1 = Check Box(
						{"Constant HS gas pCO2 (ppmv)  "}, 
				
						If( (cb1 << Get()) == 1,
							CGPCO2 << Enabled( 0 );
							bbc << Enabled( 0 );
							nb << Enabled( 1 );
						,
							CGPCO2 << Enabled( 1 );
							bbc << Enabled( 1 );
						)
					),
					,
					nb = Number Edit Box( 400, 5 )
				),
				H List Box(
					cb2 = Check Box(
						{"Constant HS ratio (Vgas:Vliquid)"},
						If( (cb2 << Get()) == 1,
							CHSR << Enabled( 0 );
							bbc2 << Enabled( 0 );
							nb << Enabled( 1 );
						,
							CHSR << Enabled( 1 );
							bbc2 << Enabled( 1 );
						)
					),
					,
					nb2 = Number Edit Box( 1.0, 5 )
				),
				H List Box(
					cb3 = Check Box(
						{"Constant Bar. Pres (kPa)"},
						If( (cb3 << Get()) == 1,
							CBP << Enabled( 0 );
							bbc3 << Enabled( 0 );
							nb3 << Enabled( 1 );
						,
							CBP << Enabled( 1 );
							bbc3 << Enabled( 1 );
						)
					),
					,
					nb3 = Number Edit Box( 101.325, 7 )
				),
				H List Box(
					cb4 = Check Box(
						{"Constant Salinity (g/kg)"},
						If( (cb4 << Get()) == 1,
							CSAL << Enabled( 0 );
							bbc4 << Enabled( 0 );
							nb4 << Enabled( 1 );
						,
							CSAL << Enabled( 1 );
							bbc4 << Enabled( 1 );
						)
					),
					,
					nb4 = Number Edit Box( 0, 4 )
				)
			),
			V List Box(
				Panel Box( "Carbonate equilibrium equation set",
					rb6 = Radio Box( {"Freshwater (Millero 1979)", "Estuarine (Millero et al. 2010)","Marine (Dickson et al. 2007)"},
					If( (rb6 << Get()) == 1, bbc4 << Enabled( 0 );CSAL<<Enable(0),bbc4 << Enabled( 1 );CSAL<<Enable(1))	
					)			
				)
			), 
			V List Box(
				Panel Box( "Solution method",
					rb5 = Radio Box( {"Analytical solutions (faster)", "Iterative solutions (slower but more stable in extreme situations)"} )
				
				)
			), 
	
		),
		H List Box(
			Spacer Box(),
			Button Box( "Remove",
				CPCO2 << RemoveSelected;
				CGPCO2 << RemoveSelected;
				CT << RemoveSelected;
				CTw << RemoveSelected;
				CAlk << RemoveSelected;
				CHSR << RemoveSelected;
				CSAL << RemoveSelected;
				CBP << RemoveSelected;		
					
			),
			Button Box( "OK",
				CPCO2var = CPCO2 << Get Items;
				CTvar = CT << Get Items;
				CTeqvar = CTw << Get Items;
				CAlkvar = CAlk << Get Items;
				trig1 = cb1 << Get();
				trig2 = cb2 << Get();
				trig3 = cb3 << Get();
				trig4 = cb4 << Get();
				Eqset=rb6<<get();
				If( trig1 == 1,
					cHSGpCO2 = nb << get(),
					HSGpCO2var = CGPCO2 << get items
				);
				If( trig2 == 1,
					cHSRc = nb2 << get(),
					HSRvar = CHSR << get items
				);
				If( trig3 == 1,
					cBPc = nb3 << get(),
					BPvar = CBP << get items
				);
				If( (trig4 == 1|Eqset==1),
					cSALc = nb4 << get(),
					SALvar = CSAL << get items
				);
				soltype = rb5 << get();
				
					
			)
		), 

	)
); 

OTemp = Column( CTvar ) << Get As Matrix;
OTeq = Column( CTeqvar ) << Get As Matrix;
OpCO2 = Column( CPCO2var ) << Get As Matrix;
OAlk = Column( CAlkvar ) << Get As Matrix;
If( trig1 == 0,
	OHSGpCO2 = Column( HSGpCO2var ) << get as matrix
);
If( trig2 == 0,
	OHSR = Column( HSRvar ) << get as matrix
);
If( trig3 == 0,
	OBP = Column( BPvar ) << get as matrix
);
If( trig4 == 0&Eqset>1,
	OSAL = Column( SALvar ) << get as matrix
);
no = N Row( Otemp );
OpCO2orig = J( no, 1, . );
OCO2orig = J( no, 1, . );
OpCO2wrong = J( no, 1, . );
nm = Loc Nonmissing( OTemp || OTeq || OpCO2 || OAlk );
nnm = N Row( nm );
Temp = OTemp[nm];
Teq = OTeq[nm];
pCO2 = OpCO2[nm];
ANC = OAlk[nm] * 1e-6;
If( trig1 == 0,
	HSGpCO2 = OHSGpCO2[nm],
	HSGpCO2 = J( nnm, 1, cHSGpCO2 )
);
If( trig2 == 0,
	HSR = OHSR[nm],
	HSR = J( nnm, 1, cHSRc )
);
If( trig4 == 0&Eqset>1,
	SAL = OSAL[nm],
	SAL = J( nnm, 1, cSALc )
);
If( trig3 == 0,
	BP = OBP[nm],
	BP = J( nnm, 1, cBPc )
);


MolVol = (0.082057 * (Temp + 273.15));
//
// Solubility coefficients from Weiss 1974 with Sal=0 for freshwater option
// Dissociation of water from Dickson and RIley (1979)
//
//

if(
Eqset==1,
//
// Millero 1979  Temp=20; Sal=0;
//
K1=10^-(-126.34048+6320.813/(Temp+273.15)+19.568224*ln(Temp+273.15));
K2=10^-(-90.18333+5143.692/(Temp+273.15)+14.613358*ln(Temp+273.15));
Kh = exp(-60.2409 + 93.4517 * (100 / (273.15 + Temp)) + 23.3585 * Ln( (273.15 + Temp) / 100 ));
Kh2 = exp(-60.2409 + 93.4517 * (100 / (273.15 + Teq)) + 23.3585 * Ln( (273.15 + Teq) / 100 ));,
kw=exp(148.9652-13847.26/(Temp+273.15)-23.6521*ln(273.15+Temp));
Eqset==3,
//
// Dickson 2007
//
K1=10^(-3633.86/ (Temp + 273.15)+61.2172-9.67770*ln(Temp+273.15)+0.011555*Sal-0.0001152*Sal^2);
K2 = 10^(-417.78/ (Temp + 273.15) - 25.9290 + 3.16967*ln(Temp+273.15)+0.01781*Sal-0.0001112*Sal^2);
kw=exp(148.9652-13847.26/(Temp+273.15)-23.6521*ln(273.15+Temp)+sqrt(Sal):*(118.67/(Temp+273.15)-5.977+1.0495*ln(273.15+Temp))-0.01615*Sal);
Kh = exp(-60.2409 + 93.4517 * (100 / (273.15 + Temp)) + 23.3585 * Ln( (273.15 + Temp) / 100 )+Sal:*(0.023517-0.023656*(273.15+Temp)/100+0.0047036*((273.15+Temp)/100)^2));
Kh2 = exp(-60.2409 + 93.4517 * (100 / (273.15 + Teq)) + 23.3585 * Ln( (273.15 + Teq) / 100 )+Sal:*(0.023517-0.023656*(273.15+Teq)/100+0.0047036*((273.15+Teq)/100)^2)),
Eqset==2,
//
// Millero et al. 2010 as amended by Orr et al. 2015
//
kw=exp(148.9652-13847.26/(Temp+273.15)-23.6521*ln(273.15+Temp)+sqrt(Sal):*(118.67/(Temp+273.15)-5.977+1.0495*ln(273.15+Temp))-0.01615*Sal);
Kh = exp(-60.2409 + 93.4517 * (100 / (273.15 + Temp)) + 23.3585 * Ln( (273.15 + Temp) / 100 )+Sal:*(0.023517-0.023656*(273.15+Temp)/100+0.0047036*((273.15+Temp)/100)^2));
Kh2 = exp(-60.2409 + 93.4517 * (100 / (273.15 + Teq)) + 23.3585 * Ln( (273.15 + Teq) / 100 )+Sal:*(0.023517-0.023656*(273.15+Teq)/100+0.0047036*((273.15+Teq)/100)^2));
pK10=(-126.34048+6320.813/(Temp+273.15)+19.568224*ln(Temp+273.15));
A1 = 13.4038*Sal^0.5 + 0.03206*Sal - 5.242e-5*Sal^2;
B1 = -530.659*Sal^0.5 - 5.8210*Sal;
C1 = -2.0664*Sal^0.5;
pK1 = pK10 + A1 + B1:/(Temp+273.15) + C1:*log(Temp+273.15);
K1 = 10^-pK1;
pK20=(-90.18333+5143.692/(Temp+273.15)+14.613358*ln(Temp+273.15));
A2 = 21.3728*Sal^0.5 + 0.1218*Sal - 3.688e-4*Sal^2;	
B2 = -788.289*Sal^0.5 - 19.189*Sal;
C2 = -3.3738*Sal^0.5;
pK2 = pK20 + A2 + B2:/(Temp+273.15) + C2:*log(Temp+273.15);
K2 = 10^-pK2;
);



nnm = N Rows( ANC );
pH = J( nnm, 1, 7 );
pHB = J( nnm, 1, . );
CO2 = Kh :* pCO2 * 1e-6;
	
//
// Chemical equilibrium calculation after equilibration
//


If(
	soltype == 2,
		f = Expr(
			Abs( tANC - (tCO2 * tk1 / (10 ^ -tpH) + 2 * tCO2 * tK1 * tK2 / ((10 ^ -tpH) ^ 2) + tKw / (10 ^ -tpH) - (10 ^ -tpH)) )
		);
		For( i = 1, i <= nnm, i++,
			tANC = ANC[i];
			tCO2 = CO2[i];
			tk1 = k1[i];
			tk2 = k2[i];
			tpH = pH[i];
			tKw = Kw[i];
			sol = Minimize( f, {tpH( 6, 10 )}, tolerance( 1e-12 ), Method( SR1 ) );
			pH[i] = tpH;
		);
		Haft = 10 ^ -pH;,
	soltype == 1, 	
	
		a = 1;
		b = ANC;
		CO2 = Kh :* pCO2 * 1e-6;
		c = -(CO2 :* K1 + Kw);
		d = -2 * CO2 :* K1 :* K2;
		p = (3 :* a :* c - b ^ 2) / (3 * a ^ 2);
		q = (2 :* b ^ 3 - 9 * a :* b :* c + 27 * d :* a ^ 2) / (27 * a ^ 3);
		qs = c / 3 - (b ^ 2) / 9;
		rs = (b :* c - 3 * d) / 6 - (b ^ 3) / 27;
		ttt = qs ^ 3 + rs ^ 2;
		s1 = Root( rs + Sqrt( qs ^ 3 + rs ^ 2 ), 3 );
		s2 = Root( rs - Sqrt( qs ^ 3 + rs ^ 2 ), 3 );
		r1 = (s1 + s2) - b / 3;
		Haft = 2 :* Root( -p / 3 ) :* Cos( ArcCosine( 3 * q :* Root( -3 :/ p ) :/ (2 * p) ) :/ 3 ) - b / 3;
		pH = -Log10( Haft );
);

a0 = 1 / (1 + K1 :/ Haft + K1 :* K2 :/ (Haft ^ 2));
DIC = CO2 :/ a0;
TCaft = pCO2 :* HSR :/ MolVol * 1e-6 + DIC;
TCbef = TCaft - HSR :* HSGpCO2 * 1e-6 :/ MolVol;
	//
//
// Chemical equilibrium calculation before equilibration
//
//
If(
	//
	// Analytical solution
	//
	soltype == 1,
		a1 = 1;
		b1 = -(ANC :* K1 - Kw + K1 :* K2 - TCbef :* K1);
		c1 = (K1 :* K2 :* ANC - K1 :* Kw - 2 * TCbef :* K1 :* K2) :* (ANC + K1) - 4 * (-K1 :* K2 :* Kw);
		d1 = 4 * (ANC :* K1 - Kw + K1 :* K2 - TCbef :* K1) :* (-K1 :* K2 :* Kw) - (K1 :* K2 :* ANC - K1 :* Kw - 2 * TCbef :* K1 :* K2) ^ 2 - (-K1 :*
		K2 :* Kw) :* (ANC + K1) ^ 2;
		p2 = (3 * a1 :* c1 - b1 ^ 2) / (3 * a1 ^ 2);
		q2 = (2 * b1 ^ 3 - 9 * a1 :* b1 :* c1 + 27 * d1 :* a1 ^ 2) / (27 * a1 ^ 3);
		r2 = 2 * Root( -p2 / 3 ) :* Cos( ArcCosine( 3 * q2 :* Root( -3 / p2 ) :/ (2 * p2) ) / 3 ) - b1 / 3;
		RR = Root( (ANC + K1) ^ 2 / 4 - (ANC :* K1 - Kw + K1 :* K2 - TCbef :* K1) + r2 );
		SS = 3 :* (ANC + K1) ^ 2 / 4 - RR ^ 2 - 2 * (ANC :* K1 - Kw + K1 :* K2 - TCbef :* K1);
		TT = (4 * (ANC + K1) :* (ANC :* K1 - Kw + K1 :* K2 - TCbef :* K1) - 8 * (K1 :* K2 :* ANC - K1 :* Kw - 2 * TCbef :* K1 :* K2) - (ANC + K1) ^ 3
		) :/ (4 * RR);
		Hbef = -(ANC + K1) :/ 4 + RR / 2 + Root( SS + TT ) :/ 2;
		pHB = -Log10( Hbef );, 
	//
	//  Iterative solution  
	//
	soltype == 2,
		f2 = Expr(
			Abs(
				1e30 * ((10 ^ -tpHB) ^ 4 + (tANC + tk1) * (10 ^ -tpHB) ^ 3 + (tANC * tk1 - tkw + tk1 * tk2 - tTCbef * tk1) * (10 ^ -tpHB) ^ 2 + (tk1
				 * tk2 * tANC - tk1 * tkw - 2 * tTCBef * tk1 * tk2) * (10 ^ -tpHB) - tk1 * tk2 * tkw)
			)
		);
		For( i = 1, i <= nnm, i++,
			tANC = ANC[i];
			tCO2 = CO2[i];
			tk1 = k1[i];
			tk2 = k2[i];
			tpHB = pH[i] - .2;
			tKw = Kw[i];
			tTCbef = TCbef[i];
			tHB = 10 ^ -pH[i];
			sol = Minimize( f2, {tpHB( 6, 11 )}, tolerance( 1e-20 ), Method( NR ), MaxIter( 500 ) );
			pHB[i] = tpHB;
		);
		hbef = 10 ^ -pHB;
);
//
//Final calculations and result output
//
aa0 = 1 / (1 + K1 :/ Hbef + K1 :* K2 :/ (Hbef ^ 2));
pCO2orig = aa0 :* TCbef * 1e6 :/ kh2;
CO2orig = aa0 :* TCbef * 1e6;
pCO2wrong = ((pCO2 - HSGpCO2) :* HSR :/ MolVol + pCO2 :* Kh) :/ kh2;
OpCO2orig[nm] = pCO2orig:*BP/101.325;
OCO2orig[nm] = CO2orig:*BP/101.325;
OpCO2wrong[nm] = pCO2wrong:*BP/101.325;;
//dt << New Column( "pH after equil.", values( pH ) );
//dt << New Column( "pH before equil.", values( pHB ) );

dt << New Column( "uncorrected pCO2 (uatm)", values( OpCO2wrong ) );
dt << New Column( "corrected pCO2 (uatm)", values( OpCO2orig ) );
dt << New Column( "corrected CO2 (umol/L)", values( OCO2orig ) );