Computes, in closed form, the four partial derivatives of the Bolin sensitivity \(B_e = \partial DIC / \partial [CO_2^*]\) (at constant total alkalinity, temperature and salinity) with respect to temperature, salinity, DIC and total alkalinity.
buffderiv(flag, var1, var2, S = 35, T = 25, Patm = 1, Pt = 0, Sit = 0,
NH4t = 0, HSt = 0, k1k2 = "l", kf = "dg", ks = "d",
pHscale = "T", b = "u74", warn = "y", npolish = 3)A data frame containing:
Bolin sensitivity, dimensionless
derivative of Be with respect to temperature, per degree Celsius
derivative of Be with respect to salinity, per practical salinity unit
derivative of Be with respect to DIC, per (mol/kg)
derivative of Be with respect to total alkalinity, per (mol/kg)
derivative of [H+] with respect to temperature, (mol/kg) per degree Celsius
derivative of [H+] with respect to salinity, (mol/kg) per psu
derivative of [H+] with respect to DIC, dimensionless
derivative of [H+] with respect to total alkalinity, dimensionless
Ac(h) + Anc(h) - ALK (mol/kg). A diagnostic. It must vanish to machine precision; if it does not, the alkalinity being differentiated is not the one seacarb solved, and every derivative above is suspect.
select the pair of variables available, as in carb
value of the first variable in mol/kg, except for pH, and for pCO2 in \(\mu\)atm
value of the second variable in mol/kg, except for pH
salinity
temperature in degrees Celsius
surface atmospheric pressure in atm, default is 1 atm
concentration of total dissolved inorganic phosphorus (mol/kg)
concentration of total dissolved inorganic silicon (mol/kg)
concentration of total dissolved inorganic ammonia (mol/kg)
concentration of total dissolved inorganic hydrogen sulfide (mol/kg)
formulation of K1 and K2. Must be "l" (Lueker et al., 2000). See Details.
formulation of Kf. Must be "dg" (Dickson and Riley, 1979). See Details.
formulation of Ks. Must be "d" (Dickson, 1990).
choice of pH scale. Must be "T" (total scale).
formulation of total boron. Must be "u74" (Uppstrom, 1974).
"y" to show warnings when T or S go beyond the valid range
number of Newton refinement steps applied to [H+] (default 3). See Details.
James Orr james.orr@lsce.ipsl.fr
Surface only (P = 0).
Which variables are held constant. Two distinct sets are involved, and it is worth
separating them. \(B_e\) is itself defined as a partial derivative at constant
\(A_T\), \(T\) and \(S\); that clause is part of the definition of the quantity,
and fixes which perturbation \(B_e\) measures. Once defined, however,
\(B_e\) is simply a function of the state, \(B_e(C_T,A_T,T,S)\),
and each of the four derivatives returned here is an ordinary partial derivative of that
function, holding the other three state variables constant. So dBe_dT is taken
at constant \(C_T\), \(A_T\) and \(S\); dBe_dAT at constant
\(C_T\), \(T\) and \(S\); and so on. This is the usual situation for a
second-order thermodynamic quantity, and is no different from differentiating the heat
capacity \(C_p = (\partial H/\partial T)_p\), itself defined at constant pressure, with
respect to pressure.
The choice of state variables does matter, because \(C_T\) and \([CO_2^*]\)
are not independent: \(\partial B_e/\partial T\) at constant \(C_T\) is not the same
number as \(\partial B_e/\partial T\) at constant \([CO_2^*]\). This routine uses
\((C_T, A_T, T, S)\) throughout, these being the variables carried by
ocean models and reported by carb.
Method. The proton concentration \(h\) is fixed implicitly by the total
alkalinity constraint
$$F(h; T,S,C_T,A_T) = A_c(h,C_T) + A_{nc}(h) - A_T = 0$$
and \(B_e\) is then an explicit algebraic function \(G(h,C_T,T,S)\),
the same one evaluated by buffsun. For any \(Y\) in
\(\{T,S,C_T,A_T\}\), the implicit function theorem gives
$$\frac{dB_e}{dY} = G_Y + G_h \left(-\frac{F_Y}{F_h}\right)$$
Every term is closed form. The only nonlinear solve is the single call to carb
(or carbfull) that supplies \(h\). This makes the routine roughly two orders
of magnitude cheaper than differencing carb() numerically, and free of step-size error.
Acid-base systems. Carbonate, borate, water, phosphate, silicate, fluoride, and
(on the carbfull path) ammonia and sulfide. Sulfate and the free proton are carried
exactly as in SolveSAPHE. The alkalinity is precisely the one that seacarb itself inverts;
alk_residual is returned so this can be checked at every point.
Silicate follows the solver. When NH4t = HSt = 0 the routine calls
carb, whose alkalinity treats silicate as monoprotic (K1si only).
Otherwise it calls carbfull, which is diprotic. K2si is set to zero on
the first path, which reduces the diprotic expressions to the monoprotic ones exactly.
Constants are not optional. The analytic dK/dT and dK/dS (in dlnK.R) are the
exact derivatives of specific formulations, so buffderiv stops rather than silently
returning derivatives of the wrong function. Two of the guards matter in practice:
k1k2 = "x" silently switches to Waters et al. (2014) when T < 2 or T > 35 or S < 19 or
S > 43, which is exactly where polar surface waters lie; and kf = "x" silently switches
between Perez and Fraga and Dickson and Riley by T and S range, changing the SWS-to-total
conversion and hence Kw, Ksi, K1p, K2p and K3p by 0.5 to 0.9 percent. Neither switch warns.
Newton refinement. carb() returns \(h\) converged to roughly 1e-10
relative. \(B_e\) is steep in \(h\) and its derivatives steeper still, so
buffderiv applies npolish Newton steps to \(F(h) = 0\) using \(F_h\), which
it needs anyway. Speciation, \(B_e\) and all four derivatives then follow from the same
\(h\) to machine precision. Set npolish = 0 to use carb()'s \(h\) untouched.
Bolin, B. and Eriksson, E., 1959 Changes in the carbon dioxide content of the atmosphere and sea due to fossil fuel combustion. The Atmosphere and the Sea in Motion, Rossby Memorial Volume, Rockefeller Institute Press, 130-142.
Munhoven, G., 2013 Mathematics of the total alkalinity-pH equation - pathway to robust and universal solution algorithms: the SolveSAPHE package v1.0.1. Geoscientific Model Development 6, 1367-1388, tools:::Rd_expr_doi("10.5194/gmd-6-1367-2013").
Orr, J. C., Epitalon, J.-M. and Gattuso, J.-P., 2015 Comparison of ten packages that compute ocean carbonate chemistry. Biogeosciences 12, 1483-1510, tools:::Rd_expr_doi("10.5194/bg-12-1483-2015").
buffsun, buffer, buffesm, buffzwg,
carb, carbfull, derivnum.
## Bolin sensitivity and its four partial derivatives,
## from total alkalinity (2300 umol/kg) and DIC (2000 umol/kg)
buffderiv(flag = 15, var1 = 2300e-6, var2 = 2000e-6, S = 35, T = 20,
Pt = 0.2e-6, Sit = 3e-6,
k1k2 = "l", kf = "dg", ks = "d", pHscale = "T", b = "u74")
## Vectorised over a field. The single carb() call is the only nonlinear solve.
n <- 1000
S <- runif(n, 30, 37)
T <- runif(n, -1.8, 30)
AT <- runif(n, 2100, 2400) * 1e-6
CT <- AT * runif(n, 0.85, 0.97)
d <- buffderiv(15, AT, CT, S = S, T = T, k1k2 = "l", kf = "dg", ks = "d",
pHscale = "T", b = "u74", warn = "n")
summary(d[, c("Be", "dBe_dT", "dBe_dS", "dBe_dCT", "dBe_dAT")])
## The alkalinity residual must be at machine precision
max(abs(d$alk_residual) / AT)
Run the code above in your browser using DataLab