#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ carbonate-solver.py -- Science Journaling Club, Volume 1, Issue 3, Spring 2025. THE QUESTION ------------ The seawater carbonate system is four coupled equilibria (carbon dioxide dissolution, the first and second dissociations of carbonic acid, and the dissociation of water) tied together by one conservation statement, total alkalinity. Almost nobody solves it correctly by hand, because the honest version is a quartic in the hydrogen ion concentration and the usual shortcut is to drop the second dissociation. What happens to surface seawater pH, to the carbonate ion concentration, and to the aragonite saturation state as atmospheric carbon dioxide rises from 180 to 1000 parts per million, when the system is solved properly rather than approximated? And how large is the error the common approximation introduces? ONE RESULT UP FRONT, BECAUSE IT IS A DISAGREEMENT ------------------------------------------------- Solved at full atmospheric equilibrium, the rise in hydrogen ion concentration from 280 to 420 ppm is 39.80 percent, not the 30 percent that is usually quoted. Section 2d of the output debugs that gap in four steps. It is not the solver and it is not the uncertainty in the constants. Asked the question the quoted figure was originally answering, 280 to 370 ppm, this solver returns a pH decline of 0.0993 units against the published observational estimate of 0.11 +/- 0.03. The quoted figure is anchored to a 370 ppm atmosphere and has not been updated. The number is reported as it comes out. WHAT THIS PROGRAM IS -------------------- This is a computation, not an observation. The Science Journaling Club has no ship, no bottle sampler, no spectrophotometric pH cell and no ocean. Nothing below was measured in seawater by us. Every number this file prints is the output of a numerical solver applied to published thermodynamic constants. Where we write "measured" we mean "computed from the model". The experiment is the code. THE MODEL --------- State variables: total alkalinity TA (mol/kg-SW), temperature T (deg C), practical salinity S, and the dry-air mole fraction of carbon dioxide xCO2 (ppm) in equilibrium with the surface. From those we solve for [H+]. fugacity fCO2 = xCO2 * 1e-6 * (P - pH2O) * phi(T) Henry's law [CO2*] = K0 * fCO2 first acid [HCO3-] = K1 [CO2*] / [H+] second acid [CO3^2-] = K1 K2 [CO2*] / [H+]^2 water [OH-] = Kw / [H+] borate [B(OH)4-] = KB * BT / (KB + [H+]) alkalinity TA = [HCO3-] + 2[CO3^2-] + [B(OH)4-] + [OH-] - [H+] Substituting the speciation into the alkalinity definition gives one nonlinear equation in [H+], strictly decreasing on (0, inf), which we solve by a safeguarded Newton iteration and check against bisection and against the roots of the equivalent quartic polynomial. Constants, all on the total hydrogen ion pH scale, mol/(kg-SW): K0 Weiss (1974), CO2 solubility, mol/(kg atm), defined on fugacity K1,K2 Lueker, Dickson & Keeling (2000), the Mehrbach refit Kw Millero (1995) as given in the Dickson best-practice guide KB Dickson (1990) Ksp aragonite and calcite, Mucci (1983) BT Uppstrom (1974) boron-to-salinity ratio Ca Riley & Tongudai (1967) calcium-to-salinity ratio pH2O Weiss & Price (1980) phi virial fugacity coefficient, Weiss (1974) ASSUMPTIONS ----------- 1. Surface water at one atmosphere total pressure. No pressure correction is applied to any constant, so every result is a surface result. 2. Thermodynamic equilibrium with the atmosphere. Real surface water is almost never in equilibrium; air-sea gas exchange takes months to a year for a mixed layer, and the model grants it instantly. 3. Total alkalinity is held fixed while carbon dioxide is added. This is the correct bookkeeping for dissolving CO2 into seawater, since CO2 is a neutral species and does not change the charge balance, but it is wrong on timescales long enough for carbonate dissolution, river input or calcification to change alkalinity. 4. Alkalinity contains only carbonate, borate, water and the free proton. Phosphate, silicate, ammonia, sulfide, fluoride and the bisulfate proton are all omitted. In open-ocean surface water they together account for well under half a percent of TA. 5. A single well-mixed parcel. No biology, no vertical mixing, no upwelling, no seasonality, no spatial structure of any kind. 6. The published constants are treated as exact in the deterministic run. The Monte Carlo section relaxes that assumption. LIMITATIONS ----------- The aragonite saturation state computed here is a thermodynamic quantity, not a biological one. Omega below one means aragonite is thermodynamically unstable in that water; it does not mean a named organism dies, and many calcifiers maintain internal chemistry quite different from the water around them. The model cannot say anything about the deep ocean, about the saturation horizon, or about any real place, because it has no pressure term and no geography. The constants are calibrated over roughly 2-35 degrees C and salinity 19-43; results at 0 degrees C are a modest extrapolation of K1 and K2 and are flagged where they appear. SEED ---- Master seed 20250320, hard-coded below and used for every random draw through numpy's SeedSequence. The whole output is deterministic. Run: python carbonate-solver.py > carbonate-solver-output.txt """ import math import platform import sys import time import numpy as np SEED = 20250320 # ---------------------------------------------------------------------------- # physical constants and ratios # ---------------------------------------------------------------------------- RGAS = 82.05736 # cm^3 atm / (mol K), for the fugacity coefficient TK0 = 273.15 P_TOT = 1.0 # atm, total pressure at the sea surface def kelvin(t_c): return np.asarray(t_c, dtype=float) + TK0 # --- CO2 solubility, Weiss (1974) eq. 12, mol/(kg-SW atm), fugacity based ---- def K0_weiss(t_c, sal): T = kelvin(t_c) t100 = T / 100.0 ln_k0 = (93.4517 / t100 - 60.2409 + 23.3585 * np.log(t100) + sal * (0.023517 - 0.023656 * t100 + 0.0047036 * t100 * t100)) return np.exp(ln_k0) # --- first and second dissociation, Lueker, Dickson & Keeling (2000) --------- def K1_lueker(t_c, sal): T = kelvin(t_c) pk1 = (3633.86 / T - 61.2172 + 9.6777 * np.log(T) - 0.011555 * sal + 0.0001152 * sal * sal) return np.power(10.0, -pk1) def K2_lueker(t_c, sal): T = kelvin(t_c) pk2 = (471.78 / T + 25.9290 - 3.16967 * np.log(T) - 0.01781 * sal + 0.0001122 * sal * sal) return np.power(10.0, -pk2) # --- ion product of water, Millero (1995), total scale ----------------------- def Kw_millero(t_c, sal): T = kelvin(t_c) s = np.asarray(sal, dtype=float) lnkw = (148.9802 - 13847.26 / T - 23.6521 * np.log(T) + (-5.977 + 118.67 / T + 1.0495 * np.log(T)) * np.sqrt(s) - 0.01615 * s) return np.exp(lnkw) # --- boric acid, Dickson (1990), total scale -------------------------------- def KB_dickson(t_c, sal): T = kelvin(t_c) s = np.asarray(sal, dtype=float) sq = np.sqrt(s) lnkb = ((-8966.90 - 2890.53 * sq - 77.942 * s + 1.728 * s ** 1.5 - 0.0996 * s * s) / T + 148.0248 + 137.1942 * sq + 1.62142 * s + (-24.4344 - 25.085 * sq - 0.2474 * s) * np.log(T) + 0.053105 * sq * T) return np.exp(lnkb) # --- aragonite solubility product, Mucci (1983) ----------------------------- def Ksp_aragonite(t_c, sal): T = kelvin(t_c) s = np.asarray(sal, dtype=float) log10ksp = (-171.945 - 0.077993 * T + 2903.293 / T + 71.595 * np.log10(T) + (-0.068393 + 0.0017276 * T + 88.135 / T) * np.sqrt(s) - 0.10018 * s + 0.0059415 * s ** 1.5) return np.power(10.0, log10ksp) # --- calcite solubility product, Mucci (1983), for the comparison column ---- def Ksp_calcite(t_c, sal): T = kelvin(t_c) s = np.asarray(sal, dtype=float) log10ksp = (-171.9065 - 0.077993 * T + 2839.319 / T + 71.595 * np.log10(T) + (-0.77712 + 0.0028426 * T + 178.34 / T) * np.sqrt(s) - 0.07711 * s + 0.0041249 * s ** 1.5) return np.power(10.0, log10ksp) # --- conservative ratios ---------------------------------------------------- def boron_total(sal): """Uppstrom (1974): 0.000232 kg-B per kg-SW at chlorinity 19.374.""" return 0.000232 / 10.811 * (np.asarray(sal, dtype=float) / 1.80655) def calcium_total(sal): """Riley & Tongudai (1967), the ratio Mucci (1983) used.""" return 0.02128 / 40.087 * (np.asarray(sal, dtype=float) / 1.80655) # --- water vapour pressure over seawater, Weiss & Price (1980), atm --------- def p_h2o(t_c, sal): T = kelvin(t_c) s = np.asarray(sal, dtype=float) return np.exp(24.4543 - 67.4509 * (100.0 / T) - 4.8489 * np.log(T / 100.0) - 0.000544 * s) # --- fugacity coefficient for CO2 in air, Weiss (1974) ---------------------- def fugacity_coefficient(t_c, x_co2_ppm): T = kelvin(t_c) B = (-1636.75 + 12.0408 * T - 3.27957e-2 * T * T + 3.16528e-5 * T ** 3) delta = 57.7 - 0.118 * T x = np.asarray(x_co2_ppm, dtype=float) * 1e-6 return np.exp(P_TOT * (B + 2.0 * (1.0 - x) ** 2 * delta) / (RGAS * T)) def fco2_from_xco2(x_co2_ppm, t_c, sal): """Dry-air mole fraction (ppm) to fugacity (atm) at the sea surface.""" pw = p_h2o(t_c, sal) pco2 = np.asarray(x_co2_ppm, dtype=float) * 1e-6 * (P_TOT - pw) return pco2 * fugacity_coefficient(t_c, x_co2_ppm), pco2, pw # ---------------------------------------------------------------------------- # the solver # ---------------------------------------------------------------------------- def alk_residual(H, C, TA, K1, K2, KB, BT, Kw): """TA(H) minus TA_input. Strictly decreasing in H on (0, inf).""" return (K1 * C / H + 2.0 * K1 * K2 * C / (H * H) + KB * BT / (KB + H) + Kw / H - H - TA) def alk_derivative(H, C, TA, K1, K2, KB, BT, Kw): return (-K1 * C / (H * H) - 4.0 * K1 * K2 * C / (H ** 3) - KB * BT / ((KB + H) ** 2) - Kw / (H * H) - 1.0) def solve_H(C, TA, K1, K2, KB, BT, Kw, lo=1e-12, hi=1e-2, itmax=200): """Safeguarded Newton on [H+], vectorised. Returns (H, iterations).""" arrs = np.broadcast_arrays( *[np.asarray(v, dtype=float) for v in (C, TA, K1, K2, KB, BT, Kw)]) C, TA, K1, K2, KB, BT, Kw = arrs shape = C.shape lo = np.full(shape, lo, dtype=float) hi = np.full(shape, hi, dtype=float) H = np.sqrt(lo * hi) used = 0 for it in range(itmax): used = it + 1 F = alk_residual(H, C, TA, K1, K2, KB, BT, Kw) lo = np.where(F > 0.0, H, lo) # F decreasing: F>0 means H too small hi = np.where(F < 0.0, H, hi) dF = alk_derivative(H, C, TA, K1, K2, KB, BT, Kw) Hn = H - F / dF bad = ~np.isfinite(Hn) | (Hn <= lo) | (Hn >= hi) Hn = np.where(bad, np.sqrt(lo * hi), Hn) step = np.max(np.abs(Hn - H) / H) H = Hn if step < 1e-16: break return H, used def solve_H_bisect(C, TA, K1, K2, KB, BT, Kw, lo=1e-12, hi=1e-2, itmax=400): """Pure bisection in log space, scalar, kept as an independent check.""" def f(h): return alk_residual(h, C, TA, K1, K2, KB, BT, Kw) a, b = lo, hi if f(a) <= 0 or f(b) >= 0: raise ValueError("bracket does not contain a sign change") steps = itmax for i in range(itmax): m = math.sqrt(a * b) fm = f(m) if fm > 0: a = m elif fm < 0: b = m else: return m, i + 1 if (b - a) <= 4.0 * np.finfo(float).eps * m: steps = i + 1 break return math.sqrt(a * b), steps def solve_H_polynomial(C, TA, K1, K2, KB, BT, Kw): """The same equation cleared of denominators, solved by numpy.roots. H^4 + (KB+TA) H^3 + (TA*KB - K1*C - KB*BT - Kw) H^2 - (K1*C*KB + 2*K1*K2*C + Kw*KB) H - 2*K1*K2*C*KB = 0 """ coef = [1.0, KB + TA, TA * KB - K1 * C - KB * BT - Kw, -(K1 * C * KB + 2.0 * K1 * K2 * C + Kw * KB), -2.0 * K1 * K2 * C * KB] r = np.roots(coef) real = r[np.abs(r.imag) < 1e-18].real real = real[real > 0] if real.size == 0: real = r.real[r.real > 0] return float(np.max(real)) def speciate(x_co2_ppm, t_c, sal, TA, k=None): """Full solve. Concentrations in mol/kg-SW.""" if k is None: k = dict(K0=K0_weiss(t_c, sal), K1=K1_lueker(t_c, sal), K2=K2_lueker(t_c, sal), KB=KB_dickson(t_c, sal), Kw=Kw_millero(t_c, sal), Ksp=Ksp_aragonite(t_c, sal)) BT = boron_total(sal) Ca = calcium_total(sal) fco2, pco2, pw = fco2_from_xco2(x_co2_ppm, t_c, sal) C = k["K0"] * fco2 H, iters = solve_H(C, TA, k["K1"], k["K2"], k["KB"], BT, k["Kw"]) hco3 = k["K1"] * C / H co3 = k["K1"] * k["K2"] * C / (H * H) oh = k["Kw"] / H boh4 = k["KB"] * BT / (k["KB"] + H) dic = C + hco3 + co3 omega = Ca * co3 / k["Ksp"] resid = alk_residual(H, C, TA, k["K1"], k["K2"], k["KB"], BT, k["Kw"]) out = dict(H=H, pH=-np.log10(H), co2=C, hco3=hco3, co3=co3, dic=dic, oh=oh, boh4=boh4, omega=omega, fco2=fco2, pco2=pco2, pH2O=pw, resid=resid, iters=iters, BT=BT, Ca=Ca) out.update(k) return out # --- the approximations we are here to indict ------------------------------- def solve_H_no_K2(C, TA, K1, KB, BT, Kw): """Approximation A: drop the 2[CO3] term from the alkalinity balance. TA = K1 C / H + KB BT/(KB+H) + Kw/H - H. Borate and water kept. """ def f(h): return K1 * C / h + KB * BT / (KB + h) + Kw / h - h - TA a, b = 1e-12, 1e-2 for _ in range(400): m = math.sqrt(a * b) if f(m) > 0: a = m else: b = m if (b - a) <= 4.0 * np.finfo(float).eps * m: break return math.sqrt(a * b) def solve_H_carbonic_only(C, TA, K1): """Approximation B: the schoolroom version, TA = [HCO3-] alone. Then H = K1 C / TA exactly, one line of algebra and no iteration. """ return K1 * C / TA def banner(text): print() print("=" * 78) print(text) print("=" * 78) def main(): t_start = time.time() seq = np.random.SeedSequence(SEED) print("=" * 78) print("SOLVING SEAWATER: a carbonate equilibrium solver") print("Science Journaling Club, Volume 1 Issue 3, Spring 2025") print("=" * 78) print("python : %s" % sys.version.split()[0]) print("numpy : %s" % np.__version__) print("platform : %s" % platform.platform()) print("master seed : %d" % SEED) print("float eps : %.6e" % np.finfo(float).eps) print() print("This is a computation. No seawater was sampled. Every number below") print("comes from a numerical solver applied to published constants.") TREF = 18.0 # deg C, close to the area-weighted mean surface ocean SREF = 35.0 # practical salinity TAREF = 2300e-6 # mol/kg-SW, representative open-ocean surface alkalinity # ------------------------------------------------------------------------ banner("1. EQUILIBRIUM CONSTANTS AT THE REFERENCE CONDITION") print("Reference seawater: T = %.1f degC, S = %.1f, TA = %.1f umol/kg" % (TREF, SREF, TAREF * 1e6)) print() kref = dict(K0=float(K0_weiss(TREF, SREF)), K1=float(K1_lueker(TREF, SREF)), K2=float(K2_lueker(TREF, SREF)), KB=float(KB_dickson(TREF, SREF)), Kw=float(Kw_millero(TREF, SREF)), Ksp=float(Ksp_aragonite(TREF, SREF))) print(" K0 (mol/kg/atm) = %.6e pK0 = %.4f" % (kref["K0"], -math.log10(kref["K0"]))) print(" K1 (mol/kg) = %.6e pK1 = %.4f" % (kref["K1"], -math.log10(kref["K1"]))) print(" K2 (mol/kg) = %.6e pK2 = %.4f" % (kref["K2"], -math.log10(kref["K2"]))) print(" KB (mol/kg) = %.6e pKB = %.4f" % (kref["KB"], -math.log10(kref["KB"]))) print(" Kw (mol/kg)^2 = %.6e pKw = %.4f" % (kref["Kw"], -math.log10(kref["Kw"]))) print(" Ksp aragonite = %.6e pKsp = %.4f" % (kref["Ksp"], -math.log10(kref["Ksp"]))) print(" Ksp calcite = %.6e" % float(Ksp_calcite(TREF, SREF))) print(" total boron BT = %.6e mol/kg" % float(boron_total(SREF))) print(" total calcium Ca = %.6e mol/kg" % float(calcium_total(SREF))) print(" pH2O at ref = %.6e atm" % float(p_h2o(TREF, SREF))) print(" fugacity coeff. = %.6f (dimensionless)" % float(fugacity_coefficient(TREF, 420.0))) print() print("At 25 degC and S = 35, where the literature usually tabulates them:") print(" pK1(25,35) = %.4f" % (-math.log10(float(K1_lueker(25.0, 35.0))))) print(" pK2(25,35) = %.4f" % (-math.log10(float(K2_lueker(25.0, 35.0))))) print(" pKB(25,35) = %.4f" % (-math.log10(float(KB_dickson(25.0, 35.0))))) print(" pKw(25,35) = %.4f" % (-math.log10(float(Kw_millero(25.0, 35.0))))) print(" K0 (25,35) = %.6e mol/kg/atm" % float(K0_weiss(25.0, 35.0))) print(" Ksp_ar(25,35) = %.6e" % float(Ksp_aragonite(25.0, 35.0))) print() print(" A note on what xCO2 means here. The atmospheric figure quoted in") print(" the news is a dry-air mole fraction. Getting from it to the") print(" fugacity the solubility constant wants costs two corrections:") print(" subtract the water vapour pressure, then apply the virial") print(" fugacity coefficient. At %.0f degC those two together take" % TREF) fc, pc, pw = fco2_from_xco2(420.0, TREF, SREF) print(" 420 ppm to fCO2 = %.3f uatm, a reduction of %.2f percent." % (float(fc) * 1e6, 100.0 * (1.0 - float(fc) / 420e-6))) # ------------------------------------------------------------------------ banner("2. VALIDATION") print("2a. Reference pH values, club solver beside the accepted figures") print() accepted = [(280.0, 8.17, "pre-industrial surface ocean"), (420.0, 8.05, "present day surface ocean")] print(" %-7s %-32s %10s %10s %10s" % ("xCO2", "condition", "club pH", "accepted", "diff")) vdiff = [] for x, acc, label in accepted: rr = speciate(x, TREF, SREF, TAREF, kref) ph = float(rr["pH"]) print(" %-7.0f %-32s %10.4f %10.2f %+10.4f" % (x, label, ph, acc, ph - acc)) vdiff.append(ph - acc) print() print(" Largest disagreement: %.4f pH units, at the present-day point." % max(abs(d) for d in vdiff)) print(" The two signs are opposite, so our solver spreads the two") print(" reference points %.4f units apart while the accepted pair is" % (float(speciate(280.0, TREF, SREF, TAREF, kref)["pH"]) - float(speciate(420.0, TREF, SREF, TAREF, kref)["pH"]))) print(" 0.12 apart. That gap is not rounding, and section 2d takes it") print(" apart rather than averaging it away.") print() print(" A note on pH scales. K1, K2 and KB above are on the total") print(" hydrogen ion scale. Millero's Kw as published is on the seawater") print(" scale, and the conversion would move pKw by under 0.001. At the") print(" reference condition [OH-] is %.2f umol/kg in a 2300 umol/kg" % (float(kref["Kw"]) / float(speciate(420.0, TREF, SREF, TAREF, kref)["H"]) * 1e6)) print(" alkalinity balance, so a 0.001 shift in pKw moves TA by about") print(" %.1e umol/kg. We do not apply the conversion, and say so." % (float(kref["Kw"]) / float(speciate(420.0, TREF, SREF, TAREF, kref)["H"]) * 1e6 * (10 ** 0.001 - 1))) print() print("2b. Alkalinity definition and charge balance at the solution") print() print(" In this model the charge balance and the alkalinity definition are") print(" the same equation. Total alkalinity is defined as the excess of") print(" conservative cations over conservative anions, so writing charge") print(" balance out and substituting that definition returns exactly") print(" TA - [HCO3-] - 2[CO3=] - [B(OH)4-] - [OH-] + [H+] = 0.") print(" Here is that residual at every point of a check grid.") print() xs_chk = np.array([180.0, 280.0, 350.0, 420.0, 560.0, 700.0, 850.0, 1000.0]) ts_chk = np.array([0.0, 5.0, 10.0, 15.0, 18.0, 20.0, 25.0, 30.0]) XX, TT = np.meshgrid(xs_chk, ts_chk) Rg = speciate(XX, TT, SREF, TAREF) max_abs = float(np.max(np.abs(Rg["resid"]))) max_rel = float(np.max(np.abs(Rg["resid"]) / TAREF)) eps = float(np.finfo(float).eps) print(" grid points checked : %d" % XX.size) print(" max |residual| : %.6e mol/kg" % max_abs) print(" max |residual| / TA : %.6e" % max_rel) print(" that is : %.3f ulp of TA (1 ulp = %.3e)" % (max_rel / eps, eps)) print(" Newton iterations used : %d" % int(Rg["iters"])) print(" VERDICT: the alkalinity definition is satisfied to machine") print(" precision. The residual cannot be driven below one ulp of the") print(" largest term in the sum, and it is not.") print() print("2c. Three independent solvers on the same equation") print() r_ref = speciate(420.0, TREF, SREF, TAREF, kref) Cref = float(r_ref["co2"]) BTref = float(boron_total(SREF)) h_newton = float(r_ref["H"]) h_bis, n_bis = solve_H_bisect(Cref, TAREF, kref["K1"], kref["K2"], kref["KB"], BTref, kref["Kw"]) h_poly = solve_H_polynomial(Cref, TAREF, kref["K1"], kref["K2"], kref["KB"], BTref, kref["Kw"]) print(" safeguarded Newton [H+] = %.12e pH = %.10f" % (h_newton, -math.log10(h_newton))) print(" log-space bisection [H+] = %.12e pH = %.10f (%d steps)" % (h_bis, -math.log10(h_bis), n_bis)) print(" quartic via np.roots [H+] = %.12e pH = %.10f" % (h_poly, -math.log10(h_poly))) print(" Newton vs bisection rel diff = %.3e" % (abs(h_newton - h_bis) / h_newton)) print(" Newton vs quartic rel diff = %.3e" % (abs(h_newton - h_poly) / h_newton)) print(" The quartic is a genuinely different algorithm: it clears every") print(" denominator and hands the polynomial to an eigenvalue routine.") print() print("2d. The widely quoted 30 percent rise in hydrogen ion concentration") print() r280 = speciate(280.0, TREF, SREF, TAREF, kref) r420 = speciate(420.0, TREF, SREF, TAREF, kref) r410 = speciate(410.0, TREF, SREF, TAREF, kref) r400 = speciate(400.0, TREF, SREF, TAREF, kref) h280, h420 = float(r280["H"]), float(r420["H"]) pct = 100.0 * (h420 / h280 - 1.0) print(" [H+] at 280 ppm : %.6e mol/kg pH = %.4f" % (h280, float(r280["pH"]))) print(" [H+] at 420 ppm : %.6e mol/kg pH = %.4f" % (h420, float(r420["pH"]))) print(" ratio : %.6f" % (h420 / h280)) print(" club value : %+.2f percent" % pct) print(" quoted figure : about +30 percent") print(" difference : %+.2f percentage points" % (pct - 30.0)) print() print(" The same calculation to other end points:") for rr, xx in ((r400, 400.0), (r410, 410.0), (r420, 420.0)): pp = 100.0 * (float(rr["H"]) / h280 - 1.0) print(" 280 -> %4.0f ppm : dpH = %+.4f d[H+] = %+.2f %%" % (xx, float(rr["pH"]) - float(r280["pH"]), pp)) print() print(" WE DO NOT REPRODUCE THE QUOTED FIGURE, and we are not going to") print(" round our way to it. Solved at full atmospheric equilibrium the") print(" answer is %.2f percent, not 30. What follows is the debugging." % pct) print() print("2d-i. Debugging step one: is it the solver?") print() print(" Three algorithms agree to %.0e relative (section 2c). The" % (abs(h_newton - h_poly) / h_newton)) print(" alkalinity residual sits at %.2f ulp (section 2b). The calcite" % (max_rel / eps)) print(" solubility product our Mucci code returns at 25 degC and S = 35") print(" is %.4e, against the value tabulated from that paper of about" % float(Ksp_calcite(25.0, 35.0))) print(" 4.27e-07, and pKB at the same condition is %.4f against Dickson's" % (-math.log10(float(KB_dickson(25.0, 35.0))))) print(" published 8.5975. The arithmetic is not the problem.") print() print("2d-ii. Debugging step two: what change DOES give 30 percent?") print() print(" Invert the solver. Hold the 280 ppm baseline and bisect on the") print(" second end point until the rise in [H+] is exactly 30.00 percent.") print() def h_rise(x): return 100.0 * (float(speciate(float(x), TREF, SREF, TAREF, kref)["H"]) / h280 - 1.0) def invert(target, fn, lo_x=200.0, hi_x=900.0): for _ in range(200): mid = math.sqrt(lo_x * hi_x) if fn(mid) < target: lo_x = mid else: hi_x = mid if hi_x - lo_x < 1e-9 * mid: break return math.sqrt(lo_x * hi_x) x30 = invert(30.0, h_rise) x26 = invert(100.0 * (10.0 ** 0.1 - 1.0), h_rise) def ph_drop(x): return float(r280["pH"]) - float(speciate(float(x), TREF, SREF, TAREF, kref)["pH"]) x805 = invert(float(r280["pH"]) - 8.05, ph_drop) x817 = invert(float(speciate(280.0, TREF, SREF, TAREF, kref)["pH"]) - 8.17, ph_drop, 100.0, 600.0) print(" xCO2 giving exactly +30.00 %% rise in [H+] : %.1f ppm" % x30) print(" xCO2 giving the classic 0.10 pH decline : %.1f ppm (= %+.2f %%)" % (x26, h_rise(x26))) print(" xCO2 at which our solver reads pH 8.05 : %.1f ppm" % x805) print(" xCO2 at which our solver reads pH 8.17 : %.1f ppm" % x817) print() print(" So the quoted 30 percent corresponds, in a solver held at full") print(" equilibrium, to an atmosphere of about %.0f ppm rather than 420." % x30) print(" And the accepted pH pair (8.17, 8.05) is itself a %.1f percent" % (100.0 * (10.0 ** 0.12 - 1.0))) print(" rise, since 10^0.12 - 1 = %.4f. The commonest form of the claim," % (10.0 ** 0.12 - 1.0)) print(" a decline of 0.1 pH units, is only %.1f percent." % (100.0 * (10.0 ** 0.1 - 1.0))) print() print("2d-iii. Debugging step three: test the model against an") print(" observational synthesis at the CO2 level the quoted figure") print(" was actually computed for") print() print(" The '0.1 pH units, about 30 percent' claim entered circulation in") print(" the 2000s, when the atmosphere was near 370 ppm, not 420. Jiang") print(" and colleagues (2019) report the global surface decline from 1770") print(" to 2000 as 0.11 +/- 0.03 pH units. So we ask our solver the") print(" question that claim was answering: 280 -> 370 ppm.") print() r370 = speciate(370.0, TREF, SREF, TAREF, kref) d370 = float(r280["pH"]) - float(r370["pH"]) p370 = 100.0 * (float(r370["H"]) / h280 - 1.0) print(" club dpH, 280 -> 370 ppm : %.4f units" % d370) print(" published dpH, 1770 -> 2000 : 0.1100 +/- 0.0300 units") print(" difference : %+.4f units" % (d370 - 0.11)) print(" in units of the published s.d.: %.2f" % (abs(d370 - 0.11) / 0.03)) print(" club rise in [H+], 280 -> 370 : %+.2f percent" % p370) print(" the same claim restated : 10^0.11 - 1 = %+.2f percent" % (100.0 * (10.0 ** 0.11 - 1.0))) print() print(" At the CO2 level the claim belongs to, our solver agrees with the") print(" published decline to %.2f standard deviations, and our rise in" % (abs(d370 - 0.11) / 0.03)) print(" [H+] is %.2f percent against the claim's %.2f percent." % (p370, 100.0 * (10.0 ** 0.11 - 1.0))) print() print("2d-iv. What we conclude") print() print(" The disagreement is not in the solver and, as section 8 shows, it") print(" is not in the constants either: propagating every published") print(" uncertainty jointly moves the answer by a fifth of a percent, and") print(" the gap to be explained is %.1f percent. Two things account for" % (pct - 30.0)) print(" it, and the first is much the larger.") print() print(" One. The quoted figure is old. It was computed for an atmosphere") print(" near 370 ppm and has been repeated unchanged while the atmosphere") print(" went to 420. Asked the 370 ppm question our solver returns") print(" %.2f percent, inside the published observational interval. Asked" % p370) print(" the 420 ppm question it returns %.2f percent. The solver is doing" % pct) print(" what a solver should do. The quoted number has stopped doing what") print(" a number should do.") print() print(" Two. Our parcel is in instantaneous equilibrium with the air") print(" above it. The real surface ocean is not; a mixed layer takes") print(" months to a year to equilibrate, and the atmosphere has been") print(" rising the whole time. That pushes the observed value below the") print(" equilibrium one. Our inversion puts the effective mixing ratio") print(" behind a reading of pH 8.05 at %.0f ppm." % x805) print() print(" We report %.2f percent. Anyone quoting 30 percent for a 420 ppm" % pct) print(" atmosphere is quoting a number from a smaller atmosphere.") print() print("2e. Sensitivity of the headline number to the reference condition") print(" (the quoted 30 percent is not a constant of nature)") print() print(" %-8s %-7s %-11s %10s %10s %10s" % ("T degC", "S", "TA umol/kg", "pH(280)", "pH(420)", "d[H+] %")) sens_rows = [] for tt, ss, ta in ((0.0, 34.0, 2300.0), (5.0, 34.0, 2290.0), (10.0, 34.5, 2295.0), (18.0, 35.0, 2300.0), (20.0, 35.5, 2320.0), (25.0, 35.0, 2300.0), (29.0, 34.0, 2310.0), (18.0, 35.0, 2200.0), (18.0, 35.0, 2400.0)): a = speciate(280.0, tt, ss, ta * 1e-6) b = speciate(420.0, tt, ss, ta * 1e-6) pp = 100.0 * (float(b["H"]) / float(a["H"]) - 1.0) print(" %-8.1f %-7.1f %-11.0f %10.4f %10.4f %10.2f" % (tt, ss, ta, float(a["pH"]), float(b["pH"]), pp)) sens_rows.append((tt, ss, ta, float(a["pH"]), float(b["pH"]), pp)) lo_p = min(r[5] for r in sens_rows) hi_p = max(r[5] for r in sens_rows) print() print(" Across this block the rise in [H+] runs from %.2f to %.2f percent," % (lo_p, hi_p)) print(" while the absolute pH at 420 ppm runs from %.4f to %.4f." % (min(r[4] for r in sens_rows), max(r[4] for r in sens_rows))) # ------------------------------------------------------------------------ banner("3. THE MAIN GRID: xCO2 FROM 180 TO 1000 ppm") xs = np.array([180.0, 200.0, 240.0, 280.0, 320.0, 350.0, 380.0, 420.0, 450.0, 500.0, 560.0, 650.0, 750.0, 850.0, 1000.0]) print("Salinity %.1f, TA %.0f umol/kg throughout. Speciation in umol/kg." % (SREF, TAREF * 1e6)) grid_store = {} for tt in (5.0, 18.0, 25.0): R = speciate(xs, tt, SREF, TAREF) grid_store[tt] = R print() print(" T = %.0f degC" % tt) print(" %6s %8s %9s %10s %9s %10s %8s %8s" % ("xCO2", "pH", "CO2*", "HCO3-", "CO3=", "DIC", "Om_ar", "Om_ca")) Kspc = float(Ksp_calcite(tt, SREF)) for i, x in enumerate(xs): om_ca = float(R["Ca"]) * float(R["co3"][i]) / Kspc print(" %6.0f %8.4f %9.2f %10.2f %9.2f %10.2f %8.3f %8.3f" % (x, R["pH"][i], R["co2"][i] * 1e6, R["hco3"][i] * 1e6, R["co3"][i] * 1e6, R["dic"][i] * 1e6, R["omega"][i], om_ca)) R18 = grid_store[18.0] i180 = 0 i1000 = len(xs) - 1 print() print("Summary over the full sweep at T = 18 degC:") print(" pH : %.4f at 180 ppm -> %.4f at 1000 ppm (change %+.4f)" % (R18["pH"][i180], R18["pH"][i1000], R18["pH"][i1000] - R18["pH"][i180])) print(" [H+] : %.4e -> %.4e mol/kg (factor %.3f)" % (R18["H"][i180], R18["H"][i1000], R18["H"][i1000] / R18["H"][i180])) print(" [CO2*] : %8.2f -> %8.2f umol/kg (factor %.3f)" % (R18["co2"][i180] * 1e6, R18["co2"][i1000] * 1e6, R18["co2"][i1000] / R18["co2"][i180])) print(" [HCO3-] : %8.2f -> %8.2f umol/kg (factor %.3f)" % (R18["hco3"][i180] * 1e6, R18["hco3"][i1000] * 1e6, R18["hco3"][i1000] / R18["hco3"][i180])) print(" [CO3=] : %8.2f -> %8.2f umol/kg (factor %.3f)" % (R18["co3"][i180] * 1e6, R18["co3"][i1000] * 1e6, R18["co3"][i1000] / R18["co3"][i180])) print(" DIC : %8.2f -> %8.2f umol/kg (factor %.3f)" % (R18["dic"][i180] * 1e6, R18["dic"][i1000] * 1e6, R18["dic"][i1000] / R18["dic"][i180])) print(" Om_ar : %8.3f -> %8.3f (fraction remaining %.3f)" % (R18["omega"][i180], R18["omega"][i1000], R18["omega"][i1000] / R18["omega"][i180])) print() print(" Note the shape of it. Dissolved CO2 rises by the same factor as") print(" the atmosphere, %.2f. Bicarbonate barely moves, %+.1f percent." % (R18["co2"][i1000] / R18["co2"][i180], 100.0 * (R18["hco3"][i1000] / R18["hco3"][i180] - 1.0))) print(" Carbonate is the species that pays: %.1f percent of it is gone," % (100.0 * (1.0 - R18["co3"][i1000] / R18["co3"][i180]))) print(" and DIC has risen only %.1f percent to absorb all of that carbon." % (100.0 * (R18["dic"][i1000] / R18["dic"][i180] - 1.0))) # ------------------------------------------------------------------------ banner("4. BJERRUM FRACTIONS AND THE CROSSOVER POINTS") print("Fraction of DIC in each species as a function of pH, T = %.0f degC," % TREF) print("S = %.0f. These depend only on K1, K2 and [H+], not on alkalinity." % SREF) print() print(" %6s %10s %10s %10s" % ("pH", "CO2*", "HCO3-", "CO3=")) ph_grid = np.arange(5.0, 11.001, 0.25) Hg = 10.0 ** (-ph_grid) den = Hg * Hg + kref["K1"] * Hg + kref["K1"] * kref["K2"] a0 = Hg * Hg / den a1 = kref["K1"] * Hg / den a2 = kref["K1"] * kref["K2"] / den for i in range(len(ph_grid)): print(" %6.2f %10.6f %10.6f %10.6f" % (ph_grid[i], a0[i], a1[i], a2[i])) print() print(" crossover CO2* = HCO3- at pH = pK1 = %.4f" % (-math.log10(kref["K1"]))) print(" crossover HCO3- = CO3= at pH = pK2 = %.4f" % (-math.log10(kref["K2"]))) ph420 = float(r420["pH"]) den0 = 10.0 ** (-2 * ph420) + kref["K1"] * 10.0 ** (-ph420) + kref["K1"] * kref["K2"] print(" present-day surface pH = %.4f, which sits %.4f units below pK2." % (ph420, -math.log10(kref["K2"]) - ph420)) print(" At that pH the DIC split is CO2* %.4f, HCO3- %.4f, CO3= %.4f." % (10.0 ** (-2 * ph420) / den0, kref["K1"] * 10.0 ** (-ph420) / den0, kref["K1"] * kref["K2"] / den0)) print(" Ocean pH sits in the flat middle of the bicarbonate plateau, which") print(" is why bicarbonate hardly notices and carbonate does all the") print(" moving.") # ------------------------------------------------------------------------ banner("5. THE COST OF NEGLECTING THE SECOND DISSOCIATION") print("Approximation A: solve the alkalinity balance with the 2[CO3=] term") print(" deleted, keeping borate and water. Carbonate is") print(" then back-calculated from the H that comes out.") print("Approximation B: the schoolroom version, TA = [HCO3-] alone, which") print(" gives H = K1 [CO2*] / TA in closed form.") print() print(" %6s %9s %9s %8s %9s %8s %9s %9s" % ("xCO2", "pH exact", "pH A", "err A", "pH B", "err B", "Om errA%", "Om errB%")) approx_rows = [] for x in xs: rx = speciate(float(x), TREF, SREF, TAREF, kref) Cx = float(rx["co2"]) hA = solve_H_no_K2(Cx, TAREF, kref["K1"], kref["KB"], BTref, kref["Kw"]) hB = solve_H_carbonic_only(Cx, TAREF, kref["K1"]) co3A = kref["K1"] * kref["K2"] * Cx / (hA * hA) co3B = kref["K1"] * kref["K2"] * Cx / (hB * hB) Ca = float(rx["Ca"]) omX = float(rx["omega"]) omA = Ca * co3A / kref["Ksp"] omB = Ca * co3B / kref["Ksp"] phX = float(rx["pH"]) phA, phB = -math.log10(hA), -math.log10(hB) print(" %6.0f %9.4f %9.4f %+8.4f %9.4f %+8.4f %+9.2f %+9.2f" % (x, phX, phA, phA - phX, phB, phB - phX, 100.0 * (omA / omX - 1.0), 100.0 * (omB / omX - 1.0))) approx_rows.append((float(x), phX, phA, phB, omX, omA, omB)) eA = [abs(r[2] - r[1]) for r in approx_rows] eB = [abs(r[3] - r[1]) for r in approx_rows] oA = [abs(100.0 * (r[5] / r[4] - 1.0)) for r in approx_rows] oB = [abs(100.0 * (r[6] / r[4] - 1.0)) for r in approx_rows] print() print(" Approximation A: mean |pH error| = %.4f, worst = %.4f units" % (float(np.mean(eA)), max(eA))) print(" mean |Omega error| = %.2f %%, worst = %.2f %%" % (float(np.mean(oA)), max(oA))) print(" Approximation B: mean |pH error| = %.4f, worst = %.4f units" % (float(np.mean(eB)), max(eB))) print(" mean |Omega error| = %.2f %%, worst = %.2f %%" % (float(np.mean(oB)), max(oB))) print() print(" Read the Omega columns twice. Approximation A shifts pH by a few") print(" hundredths of a unit and still misstates the aragonite saturation") print(" state by tens of percent, because Omega runs on [H+] squared and") print(" the neglected term was carrying exactly the quantity Omega needs.") print() dh = [100.0 * (10.0 ** (-(r[2])) / 10.0 ** (-(r[1])) - 1.0) for r in approx_rows] print(" Error in [H+] from approximation A, percent, across the sweep:") print(" ", " ".join("%+.2f" % v for v in dh)) print() print(" And the same question asked the other way round. If you use") print(" approximation A to work out how much the hydrogen ion rises from") print(" 280 to 420 ppm, you get:") cA280 = float(speciate(280.0, TREF, SREF, TAREF, kref)["co2"]) cA420 = float(speciate(420.0, TREF, SREF, TAREF, kref)["co2"]) hA280 = solve_H_no_K2(cA280, TAREF, kref["K1"], kref["KB"], BTref, kref["Kw"]) hA420 = solve_H_no_K2(cA420, TAREF, kref["K1"], kref["KB"], BTref, kref["Kw"]) hB280 = solve_H_carbonic_only(cA280, TAREF, kref["K1"]) hB420 = solve_H_carbonic_only(cA420, TAREF, kref["K1"]) pctA = 100.0 * (hA420 / hA280 - 1.0) pctB = 100.0 * (hB420 / hB280 - 1.0) print(" exact solver : %+.2f %%" % pct) print(" approximation A : %+.2f %% (off by %+.2f points)" % (pctA, pctA - pct)) print(" approximation B : %+.2f %% (off by %+.2f points)" % (pctB, pctB - pct)) print(" The ratio is less wrong than the absolute values, but it is still") print(" wrong by %.1f and %.1f percentage points. An approximation that" % (pctA - pct, pctB - pct)) print(" misses the headline number by a fifth of its own size is not a") print(" shortcut. It is a different answer.") # ------------------------------------------------------------------------ banner("6. WHERE THE ARAGONITE SATURATION STATE CROSSES ONE") print("Bisection on xCO2 for Omega_aragonite = 1, at fixed TA = %.0f umol/kg." % (TAREF * 1e6)) print("Below the printed xCO2 the water is supersaturated with respect to") print("aragonite; above it, aragonite is thermodynamically unstable.") print() def omega_at(x, tt, ss, ta): return float(speciate(float(x), tt, ss, ta)["omega"]) print(" %8s %6s %14s %10s %11s" % ("T degC", "S", "xCO2 at Om=1", "pH there", "CO3= there")) cross_rows = [] for tt, ss in ((0.0, 34.0), (5.0, 34.0), (10.0, 34.5), (15.0, 35.0), (18.0, 35.0), (20.0, 35.0), (25.0, 35.0), (30.0, 35.0)): lo_x, hi_x = 50.0, 40000.0 if omega_at(hi_x, tt, ss, TAREF) > 1.0: print(" %8.1f %6.1f %14s %10s %11s" % (tt, ss, "> 40000", "-", "-")) continue if omega_at(lo_x, tt, ss, TAREF) < 1.0: print(" %8.1f %6.1f %14s %10s %11s" % (tt, ss, "< 50", "-", "-")) continue mid = lo_x for _ in range(300): mid = math.sqrt(lo_x * hi_x) if omega_at(mid, tt, ss, TAREF) > 1.0: lo_x = mid else: hi_x = mid if hi_x - lo_x < 1e-10 * mid: break xc = math.sqrt(lo_x * hi_x) rc = speciate(xc, tt, ss, TAREF) print(" %8.1f %6.1f %14.1f %10.4f %11.2f" % (tt, ss, xc, float(rc["pH"]), float(rc["co3"]) * 1e6)) cross_rows.append((tt, ss, xc, float(rc["pH"]), float(rc["co3"]) * 1e6)) print() print(" Cold water crosses first, and by a long way. Carbon dioxide is more") print(" soluble in cold water and the dissociation constants shift with") print(" temperature, so a polar parcel at the same alkalinity runs out of") print(" carbonate at a far lower atmospheric mixing ratio than a tropical") print(" one. The 0 degC row is outside the calibration range of K1 and K2") print(" and should be read as an extrapolation.") print() print() print(" A check worth making. Published model studies put surface") print(" aragonite undersaturation in the Southern Ocean at roughly 550 to") print(" 600 ppm. Our coldest row, %.1f degC at S = %.1f, crosses at" % (cross_rows[0][0], cross_rows[0][1])) print(" %.1f ppm. That is an independent result our model did not see" % cross_rows[0][2]) print(" during construction, and it lands inside the published window.") print() print(" Carbonate concentration at the crossing, for reference:") for row in cross_rows: print(" T = %5.1f degC : [CO3=] = %6.2f umol/kg at Omega = 1" % (row[0], row[4])) # ------------------------------------------------------------------------ banner("7. THE REVELLE FACTOR, COMPUTED NUMERICALLY") print("Revelle factor R = (d fCO2 / fCO2) / (d DIC / DIC) at constant TA,") print("by central difference on a 0.1 percent perturbation of xCO2.") print() print(" %8s %8s %10s %12s %10s" % ("xCO2", "T degC", "DIC", "Revelle R", "pH")) rev_rows = [] for tt in (5.0, 18.0, 25.0): for x in (180.0, 280.0, 420.0, 560.0, 1000.0): ex = 1e-3 rp = speciate(x * (1 + ex), tt, SREF, TAREF) rm = speciate(x * (1 - ex), tt, SREF, TAREF) r0 = speciate(x, tt, SREF, TAREF) dlnf = math.log(float(rp["fco2"]) / float(rm["fco2"])) dlnc = math.log(float(rp["dic"]) / float(rm["dic"])) Rv = dlnf / dlnc print(" %8.0f %8.1f %10.2f %12.3f %10.4f" % (x, tt, float(r0["dic"]) * 1e6, Rv, float(r0["pH"]))) rev_rows.append((x, tt, float(r0["dic"]) * 1e6, Rv, float(r0["pH"]))) r420_18 = [r[3] for r in rev_rows if r[0] == 420.0 and r[1] == 18.0][0] print() print(" The published range for present-day surface water is roughly 8 to") print(" 15, lowest in the warm subtropics and highest in cold") print(" high-latitude water. Our block runs %.2f to %.2f, which is wider" % (min(r[3] for r in rev_rows), max(r[3] for r in rev_rows))) print(" at both ends, and it should be: the block also contains glacial") print(" air at 180 ppm and an atmosphere at 1000 ppm that has never") print(" existed in the Quaternary. Restricting to the present-day column") r_now = [r[3] for r in rev_rows if r[0] == 420.0] print(" at 420 ppm gives %.2f to %.2f, inside the published range." % (min(r_now), max(r_now))) print(" R = %.2f at 420 ppm and 18 degC means a one percent rise in DIC" % r420_18) print(" buys a %.1f percent rise in fCO2. The buffer gets worse as CO2" % r420_18) print(" rises: R climbs from %.2f at 180 ppm to %.2f at 1000 ppm at this" % ([r[3] for r in rev_rows if r[0] == 180.0 and r[1] == 18.0][0], [r[3] for r in rev_rows if r[0] == 1000.0 and r[1] == 18.0][0])) print(" temperature, so each additional tonne is absorbed less willingly.") # ------------------------------------------------------------------------ banner("8. MONTE CARLO: WHAT SURVIVES THE UNCERTAINTY IN THE CONSTANTS") N_MC = 200000 child = seq.spawn(4) rng = np.random.default_rng(child[0]) print("Trials : %d" % N_MC) print("Generator : PCG64 from SeedSequence(%d).spawn(4)[0]" % SEED) print() print("Standard uncertainties applied, one independent normal draw each,") print("held common between the 280 ppm and 420 ppm end members, because it") print("is the same water and the same constants on both sides:") U = dict(pK0=0.002, pK1=0.0075, pK2=0.015, pKB=0.010, pKw=0.010, pKsp=0.020, TA=2.0e-6, T=0.10, S=0.010) for kk in ("pK0", "pK1", "pK2", "pKB", "pKw", "pKsp"): print(" u(%-4s) = %.4f (log10 units)" % (kk, U[kk])) print(" u(TA) = %.1f umol/kg" % (U["TA"] * 1e6)) print(" u(T) = %.2f degC" % U["T"]) print(" u(S) = %.3f" % U["S"]) print() d_pK0 = rng.normal(0.0, U["pK0"], N_MC) d_pK1 = rng.normal(0.0, U["pK1"], N_MC) d_pK2 = rng.normal(0.0, U["pK2"], N_MC) d_pKB = rng.normal(0.0, U["pKB"], N_MC) d_pKw = rng.normal(0.0, U["pKw"], N_MC) d_pKsp = rng.normal(0.0, U["pKsp"], N_MC) d_TA = rng.normal(0.0, U["TA"], N_MC) d_T = rng.normal(0.0, U["T"], N_MC) d_S = rng.normal(0.0, U["S"], N_MC) Tm = TREF + d_T Sm = SREF + d_S TAm = TAREF + d_TA km = dict( K0=K0_weiss(Tm, Sm) * 10.0 ** (-d_pK0), K1=K1_lueker(Tm, Sm) * 10.0 ** (-d_pK1), K2=K2_lueker(Tm, Sm) * 10.0 ** (-d_pK2), KB=KB_dickson(Tm, Sm) * 10.0 ** (-d_pKB), Kw=Kw_millero(Tm, Sm) * 10.0 ** (-d_pKw), Ksp=Ksp_aragonite(Tm, Sm) * 10.0 ** (-d_pKsp)) m280 = speciate(280.0, Tm, Sm, TAm, km) m420 = speciate(420.0, Tm, Sm, TAm, km) pct_mc = 100.0 * (m420["H"] / m280["H"] - 1.0) ph_mc = m420["pH"] om_mc = m420["omega"] dph_mc = m420["pH"] - m280["pH"] co3_mc = m420["co3"] * 1e6 def report(name, arr, det, unit=""): mean = float(np.mean(arr)) sd = float(np.std(arr, ddof=1)) se = sd / math.sqrt(len(arr)) lo_q, hi_q = np.percentile(arr, [2.5, 97.5]) print(" %s" % name) print(" deterministic : %12.5f %s" % (det, unit)) print(" Monte Carlo mean : %12.5f %s" % (mean, unit)) print(" standard error of mean %12.5f %s" % (se, unit)) print(" spread (1 s.d.) : %12.5f %s" % (sd, unit)) print(" 95%% interval : %12.5f to %.5f %s" % (lo_q, hi_q, unit)) print(" mean - deterministic : %+12.5f %s (%.2f SE)" % (mean - det, unit, abs(mean - det) / se if se > 0 else float("nan"))) return mean, sd, se, float(lo_q), float(hi_q) print("Results. SE is the standard error of the Monte Carlo mean, computed") print("from the trials themselves as s.d./sqrt(N).") print() mc_pct = report("rise in [H+], 280->420 ppm", pct_mc, pct, "%") print() mc_ph = report("pH at 420 ppm", ph_mc, float(r420["pH"]), "") print() mc_om = report("Omega aragonite at 420 ppm", om_mc, float(r420["omega"]), "") print() mc_dph = report("pH change 280->420 ppm", dph_mc, float(r420["pH"]) - float(r280["pH"]), "") print() mc_co3 = report("[CO3=] at 420 ppm", co3_mc, float(r420["co3"]) * 1e6, "umol/kg") print() print(" The rise in [H+] is the durable number here: its 95 percent") print(" interval is %.2f to %.2f percent, a width of only %.2f points," % (mc_pct[3], mc_pct[4], mc_pct[4] - mc_pct[3])) print(" because the same constants act on both end members and most of") print(" the error cancels. The absolute pH is softer, spread %.4f, and" % mc_ph[1]) print(" Omega is softest of all, spread %.4f on a value of %.3f, because" % (mc_om[1], mc_om[0])) print(" pKsp enters Omega directly and nothing cancels it.") print() print("8a-ii. How far out is the 30 percent figure?") print() print(" Our value : %.5f %%" % mc_pct[0]) print(" Quoted value : 30 %") print(" Gap : %.2f percentage points" % (mc_pct[0] - 30.0)) print(" In units of the Monte Carlo standard error of the mean : %.0f SE" % ((mc_pct[0] - 30.0) / mc_pct[2])) print(" In units of the Monte Carlo spread (1 s.d.) : %.1f s.d." % ((mc_pct[0] - 30.0) / mc_pct[1])) print() print(" Those two numbers say the same thing twice. Every published") print(" uncertainty in every constant, propagated jointly, moves the") print(" answer by about %.2f percent. The gap to be explained is %.1f" % (mc_pct[1], mc_pct[0] - 30.0)) print(" percent. The constants cannot account for it and we are not") print(" going to pretend they can. The explanation is in section 2d-iii:") print(" the quoted figure describes an ocean that has not finished") print(" equilibrating, and our solver describes one that has.") print() print("8b. Convergence of the Monte Carlo mean") print() print(" %10s %14s %14s %16s" % ("trials", "running mean", "running SE", "|mean-det|/SE")) checks = [100, 200, 500, 1000, 2000, 5000, 10000, 20000, 50000, 100000, 150000, 200000] conv_rows = [] csum = np.cumsum(pct_mc) csum2 = np.cumsum(pct_mc * pct_mc) for n in checks: m = csum[n - 1] / n var = (csum2[n - 1] - n * m * m) / (n - 1) se = math.sqrt(var / n) print(" %10d %14.5f %14.5f %16.2f" % (n, m, se, abs(m - pct) / se)) conv_rows.append((n, float(m), se)) print() print(" The standard error falls as 1/sqrt(N), as it must. From %d trials" % checks[0]) print(" to %d it drops by a factor of %.1f against a predicted %.1f." % (checks[-1], conv_rows[0][2] / conv_rows[-1][2], math.sqrt(checks[-1] / float(checks[0])))) # ------------------------------------------------------------------------ banner("9. FIGURE DATA (machine readable, for the article and the LaTeX)") print("FIG1 pH and Omega_aragonite against xCO2, three temperatures") print("# T_degC xCO2 pH Omega_ar") for tt in (5.0, 18.0, 25.0): R = grid_store[tt] for i, x in enumerate(xs): print("FIG1 %4.0f %6.0f %8.4f %8.4f" % (tt, x, R["pH"][i], R["omega"][i])) print() print("FIG2 speciation in umol/kg against xCO2 at T = 18 degC") print("# xCO2 CO2star HCO3 CO3 DIC") for i, x in enumerate(xs): print("FIG2 %6.0f %9.2f %9.2f %9.2f %9.2f" % (x, R18["co2"][i] * 1e6, R18["hco3"][i] * 1e6, R18["co3"][i] * 1e6, R18["dic"][i] * 1e6)) print() print("FIG3 Bjerrum fractions against pH at T = 18 degC, S = 35") print("# pH f_CO2 f_HCO3 f_CO3") for i in range(len(ph_grid)): print("FIG3 %5.2f %9.6f %9.6f %9.6f" % (ph_grid[i], a0[i], a1[i], a2[i])) print() print("FIG4 approximation error against xCO2 at T = 18 degC") print("# xCO2 pH_exact pH_A pH_B Omega_exact Omega_A Omega_B") for row in approx_rows: print("FIG4 %6.0f %9.4f %9.4f %9.4f %9.4f %9.4f %9.4f" % row) print() print("FIG5 Monte Carlo convergence of the rise in [H+] (percent)") print("# trials running_mean running_SE") trace_n = np.unique(np.round(np.logspace(2, math.log10(N_MC), 60)).astype(int)) for n in trace_n: m = csum[n - 1] / n var = (csum2[n - 1] - n * m * m) / (n - 1) se = math.sqrt(var / n) print("FIG5 %8d %10.5f %10.5f" % (n, m, se)) print() print("FIG6 Omega = 1 crossing against temperature") print("# T_degC S xCO2_cross pH CO3_umol") for row in cross_rows: print("FIG6 %5.1f %5.1f %9.1f %8.4f %8.2f" % row) print() print("FIG7 Revelle factor against xCO2") print("# xCO2 T_degC DIC RevelleR pH") for row in rev_rows: print("FIG7 %6.0f %6.1f %9.2f %8.3f %8.4f" % row) # ------------------------------------------------------------------------ banner("10. CLOSING NUMBERS") r1000 = speciate(1000.0, TREF, SREF, TAREF, kref) print("Reference condition : T = %.1f degC, S = %.1f, TA = %.0f umol/kg" % (TREF, SREF, TAREF * 1e6)) print("pH, pre-industrial 280 ppm : %.4f (accepted about 8.17)" % float(r280["pH"])) print("pH, present day 420 ppm : %.4f (accepted about 8.05)" % float(r420["pH"])) print("pH change 280 -> 420 : %+.4f units" % (float(r420["pH"]) - float(r280["pH"]))) print("rise in [H+] 280 -> 420 : %+.2f %% (quoted about 30 %%)" % pct) print(" gap to the quoted figure : %+.2f points, %.0f SE and %.1f s.d. out" % (pct - 30.0, (mc_pct[0] - 30.0) / mc_pct[2], (mc_pct[0] - 30.0) / mc_pct[1])) print(" xCO2 that would give 30%% : %.1f ppm at full equilibrium" % x30) print(" xCO2 that reads pH 8.05 : %.1f ppm at full equilibrium" % x805) print("Monte Carlo on that rise : %.2f %% +/- %.4f %% (SE), 95%% %.2f to %.2f" % (mc_pct[0], mc_pct[2], mc_pct[3], mc_pct[4])) print("Omega aragonite at 280 ppm : %.3f" % float(r280["omega"])) print("Omega aragonite at 420 ppm : %.3f" % float(r420["omega"])) print("Omega aragonite at 1000 ppm: %.3f" % float(r1000["omega"])) print("[CO3=] at 280 / 420 / 1000 : %.2f / %.2f / %.2f umol/kg" % (float(r280["co3"]) * 1e6, float(r420["co3"]) * 1e6, float(r1000["co3"]) * 1e6)) print("Carbonate lost 180 -> 1000 : %.1f %%" % (100.0 * (1.0 - R18["co3"][i1000] / R18["co3"][i180]))) print("Max alkalinity residual : %.3e mol/kg (%.3f ulp of TA)" % (max_abs, max_rel / eps)) print("Newton vs quartic agreement: %.3e relative" % (abs(h_newton - h_poly) / h_newton)) print("Worst Omega error, approx A: %.2f %%" % max(oA)) print("Worst Omega error, approx B: %.2f %%" % max(oB)) print("Revelle factor at 420, 18C : %.3f" % r420_18) print("Master seed : %d" % SEED) print("Wall clock : %.2f s" % (time.time() - t_start)) print() print("END OF OUTPUT") if __name__ == "__main__": main()