Surface only (P = 0).
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.