============================================================================== SOLVING SEAWATER: a carbonate equilibrium solver Science Journaling Club, Volume 1 Issue 3, Spring 2025 ============================================================================== python : 3.12.3 numpy : 2.4.2 platform : Windows-11-10.0.26200-SP0 master seed : 20250320 float eps : 2.220446e-16 This is a computation. No seawater was sampled. Every number below comes from a numerical solver applied to published constants. ============================================================================== 1. EQUILIBRIUM CONSTANTS AT THE REFERENCE CONDITION ============================================================================== Reference seawater: T = 18.0 degC, S = 35.0, TA = 2300.0 umol/kg K0 (mol/kg/atm) = 3.428761e-02 pK0 = 1.4649 K1 (mol/kg) = 1.229531e-06 pK1 = 5.9103 K2 (mol/kg) = 8.331038e-10 pK2 = 9.0793 KB (mol/kg) = 2.088980e-09 pKB = 8.6801 Kw (mol/kg)^2 = 3.227621e-14 pKw = 13.4911 Ksp aragonite = 6.661112e-07 pKsp = 6.1765 Ksp calcite = 4.307456e-07 total boron BT = 4.157576e-04 mol/kg total calcium Ca = 1.028457e-02 mol/kg pH2O at ref = 1.996713e-02 atm fugacity coeff. = 0.996521 (dimensionless) At 25 degC and S = 35, where the literature usually tabulates them: pK1(25,35) = 5.8472 pK2(25,35) = 8.9660 pKB(25,35) = 8.5975 pKw(25,35) = 13.2107 K0 (25,35) = 2.839188e-02 mol/kg/atm Ksp_ar(25,35) = 6.481759e-07 A note on what xCO2 means here. The atmospheric figure quoted in the news is a dry-air mole fraction. Getting from it to the fugacity the solubility constant wants costs two corrections: subtract the water vapour pressure, then apply the virial fugacity coefficient. At 18 degC those two together take 420 ppm to fCO2 = 410.182 uatm, a reduction of 2.34 percent. ============================================================================== 2. VALIDATION ============================================================================== 2a. Reference pH values, club solver beside the accepted figures xCO2 condition club pH accepted diff 280 pre-industrial surface ocean 8.1817 8.17 +0.0117 420 present day surface ocean 8.0361 8.05 -0.0139 Largest disagreement: 0.0139 pH units, at the present-day point. The two signs are opposite, so our solver spreads the two reference points 0.1455 units apart while the accepted pair is 0.12 apart. That gap is not rounding, and section 2d takes it apart rather than averaging it away. A note on pH scales. K1, K2 and KB above are on the total hydrogen ion scale. Millero's Kw as published is on the seawater scale, and the conversion would move pKw by under 0.001. At the reference condition [OH-] is 3.51 umol/kg in a 2300 umol/kg alkalinity balance, so a 0.001 shift in pKw moves TA by about 8.1e-03 umol/kg. We do not apply the conversion, and say so. 2b. Alkalinity definition and charge balance at the solution In this model the charge balance and the alkalinity definition are the same equation. Total alkalinity is defined as the excess of conservative cations over conservative anions, so writing charge balance out and substituting that definition returns exactly TA - [HCO3-] - 2[CO3=] - [B(OH)4-] - [OH-] + [H+] = 0. Here is that residual at every point of a check grid. grid points checked : 64 max |residual| : 4.336809e-19 mol/kg max |residual| / TA : 1.885569e-16 that is : 0.849 ulp of TA (1 ulp = 2.220e-16) Newton iterations used : 18 VERDICT: the alkalinity definition is satisfied to machine precision. The residual cannot be driven below one ulp of the largest term in the sum, and it is not. 2c. Three independent solvers on the same equation safeguarded Newton [H+] = 9.201560152892e-09 pH = 8.0361385304 log-space bisection [H+] = 9.201560152892e-09 pH = 8.0361385304 (54 steps) quartic via np.roots [H+] = 9.201560152893e-09 pH = 8.0361385304 Newton vs bisection rel diff = 0.000e+00 Newton vs quartic rel diff = 1.870e-14 The quartic is a genuinely different algorithm: it clears every denominator and hands the polynomial to an eigenvalue routine. 2d. The widely quoted 30 percent rise in hydrogen ion concentration [H+] at 280 ppm : 6.581793e-09 mol/kg pH = 8.1817 [H+] at 420 ppm : 9.201560e-09 mol/kg pH = 8.0361 ratio : 1.398032 club value : +39.80 percent quoted figure : about +30 percent difference : +9.80 percentage points The same calculation to other end points: 280 -> 400 ppm : dpH = -0.1277 d[H+] = +34.17 % 280 -> 410 ppm : dpH = -0.1367 d[H+] = +36.99 % 280 -> 420 ppm : dpH = -0.1455 d[H+] = +39.80 % WE DO NOT REPRODUCE THE QUOTED FIGURE, and we are not going to round our way to it. Solved at full atmospheric equilibrium the answer is 39.80 percent, not 30. What follows is the debugging. 2d-i. Debugging step one: is it the solver? Three algorithms agree to 2e-14 relative (section 2c). The alkalinity residual sits at 0.85 ulp (section 2b). The calcite solubility product our Mucci code returns at 25 degC and S = 35 is 4.2724e-07, against the value tabulated from that paper of about 4.27e-07, and pKB at the same condition is 8.5975 against Dickson's published 8.5975. The arithmetic is not the problem. 2d-ii. Debugging step two: what change DOES give 30 percent? Invert the solver. Hold the 280 ppm baseline and bisect on the second end point until the rise in [H+] is exactly 30.00 percent. xCO2 giving exactly +30.00 % rise in [H+] : 385.2 ppm xCO2 giving the classic 0.10 pH decline : 370.7 ppm (= +25.89 %) xCO2 at which our solver reads pH 8.05 : 404.4 ppm xCO2 at which our solver reads pH 8.17 : 289.4 ppm So the quoted 30 percent corresponds, in a solver held at full equilibrium, to an atmosphere of about 385 ppm rather than 420. And the accepted pH pair (8.17, 8.05) is itself a 31.8 percent rise, since 10^0.12 - 1 = 0.3183. The commonest form of the claim, a decline of 0.1 pH units, is only 25.9 percent. 2d-iii. Debugging step three: test the model against an observational synthesis at the CO2 level the quoted figure was actually computed for The '0.1 pH units, about 30 percent' claim entered circulation in the 2000s, when the atmosphere was near 370 ppm, not 420. Jiang and colleagues (2019) report the global surface decline from 1770 to 2000 as 0.11 +/- 0.03 pH units. So we ask our solver the question that claim was answering: 280 -> 370 ppm. club dpH, 280 -> 370 ppm : 0.0993 units published dpH, 1770 -> 2000 : 0.1100 +/- 0.0300 units difference : -0.0107 units in units of the published s.d.: 0.36 club rise in [H+], 280 -> 370 : +25.69 percent the same claim restated : 10^0.11 - 1 = +28.82 percent At the CO2 level the claim belongs to, our solver agrees with the published decline to 0.36 standard deviations, and our rise in [H+] is 25.69 percent against the claim's 28.82 percent. 2d-iv. What we conclude The disagreement is not in the solver and, as section 8 shows, it is not in the constants either: propagating every published uncertainty jointly moves the answer by a fifth of a percent, and the gap to be explained is 9.8 percent. Two things account for it, and the first is much the larger. One. The quoted figure is old. It was computed for an atmosphere near 370 ppm and has been repeated unchanged while the atmosphere went to 420. Asked the 370 ppm question our solver returns 25.69 percent, inside the published observational interval. Asked the 420 ppm question it returns 39.80 percent. The solver is doing what a solver should do. The quoted number has stopped doing what a number should do. Two. Our parcel is in instantaneous equilibrium with the air above it. The real surface ocean is not; a mixed layer takes months to a year to equilibrate, and the atmosphere has been rising the whole time. That pushes the observed value below the equilibrium one. Our inversion puts the effective mixing ratio behind a reading of pH 8.05 at 404 ppm. We report 39.80 percent. Anyone quoting 30 percent for a 420 ppm atmosphere is quoting a number from a smaller atmosphere. 2e. Sensitivity of the headline number to the reference condition (the quoted 30 percent is not a constant of nature) T degC S TA umol/kg pH(280) pH(420) d[H+] % 0.0 34.0 2300 8.1756 8.0173 44.00 5.0 34.0 2290 8.1797 8.0245 42.96 10.0 34.5 2295 8.1823 8.0308 41.77 18.0 35.0 2300 8.1817 8.0361 39.80 20.0 35.5 2320 8.1823 8.0386 39.22 25.0 35.0 2300 8.1771 8.0369 38.11 29.0 34.0 2310 8.1784 8.0408 37.25 18.0 35.0 2200 8.1653 8.0192 39.99 18.0 35.0 2400 8.1973 8.0523 39.62 Across this block the rise in [H+] runs from 37.25 to 44.00 percent, while the absolute pH at 420 ppm runs from 8.0173 to 8.0523. ============================================================================== 3. THE MAIN GRID: xCO2 FROM 180 TO 1000 ppm ============================================================================== Salinity 35.0, TA 2300 umol/kg throughout. Speciation in umol/kg. T = 5 degC xCO2 pH CO2* HCO3- CO3= DIC Om_ar Om_ca 180 8.3406 9.27 1802.83 198.03 2010.12 2.984 4.726 200 8.3027 10.30 1835.85 184.81 2030.96 2.785 4.411 240 8.2361 12.35 1889.88 163.21 2065.45 2.460 3.895 280 8.1788 14.41 1932.34 146.25 2093.01 2.204 3.491 320 8.1285 16.47 1966.66 132.56 2115.69 1.998 3.164 350 8.0943 18.02 1988.39 123.89 2130.29 1.867 2.957 380 8.0628 19.56 2007.42 116.30 2143.28 1.753 2.776 420 8.0240 21.62 2029.39 107.54 2158.56 1.621 2.567 450 7.9971 23.17 2043.80 101.80 2168.77 1.534 2.430 500 7.9558 25.74 2064.64 93.50 2183.88 1.409 2.232 560 7.9109 28.83 2085.54 85.18 2199.55 1.284 2.033 650 7.8514 33.46 2110.71 75.17 2219.33 1.133 1.794 750 7.7937 38.61 2132.50 66.50 2237.61 1.002 1.587 850 7.7429 43.76 2149.78 59.63 2253.17 0.899 1.423 1000 7.6763 51.48 2169.89 51.64 2273.01 0.778 1.232 T = 18 degC xCO2 pH CO2* HCO3- CO3= DIC Om_ar Om_ca 180 8.3324 6.03 1593.39 285.41 1884.83 4.407 6.814 200 8.2972 6.70 1632.55 269.65 1908.89 4.163 6.438 240 8.2352 8.04 1698.29 243.17 1949.49 3.754 5.806 280 8.1817 9.38 1751.53 221.70 1982.61 3.423 5.293 320 8.1345 10.72 1795.68 203.89 2010.29 3.148 4.868 350 8.1024 11.72 1824.20 192.39 2028.31 2.970 4.593 380 8.0727 12.72 1849.54 182.16 2044.42 2.812 4.349 420 8.0361 14.06 1879.28 170.15 2063.49 2.627 4.063 450 8.0107 15.07 1899.05 162.16 2076.28 2.504 3.872 500 7.9715 16.74 1928.07 150.44 2095.26 2.323 3.592 560 7.9289 18.75 1957.68 138.48 2114.91 2.138 3.306 650 7.8722 21.77 1994.05 123.78 2139.60 1.911 2.955 750 7.8170 25.11 2026.24 110.77 2162.12 1.710 2.645 850 7.7682 28.46 2052.23 100.26 2180.95 1.548 2.394 1000 7.7041 33.49 2083.05 87.80 2204.33 1.356 2.096 T = 25 degC xCO2 pH CO2* HCO3- CO3= DIC Om_ar Om_ca 180 8.3220 4.94 1473.80 334.60 1813.34 5.309 8.055 200 8.2882 5.49 1514.89 318.16 1838.53 5.048 7.659 240 8.2286 6.58 1584.74 290.15 1881.47 4.604 6.985 280 8.1771 7.68 1642.18 267.05 1916.92 4.237 6.429 320 8.1317 8.78 1690.45 247.61 1946.84 3.929 5.961 350 8.1008 9.60 1721.96 234.91 1966.47 3.727 5.655 380 8.0722 10.42 1750.18 223.51 1984.12 3.546 5.380 420 8.0369 11.52 1783.58 210.02 2005.12 3.332 5.056 450 8.0124 12.35 1805.96 200.97 2019.28 3.189 4.838 500 7.9745 13.72 1839.08 187.56 2040.36 2.976 4.515 560 7.9333 15.36 1873.21 173.74 2062.31 2.757 4.182 650 7.8783 17.83 1915.63 156.54 2090.01 2.484 3.768 750 7.8247 20.58 1953.66 141.11 2115.35 2.239 3.397 850 7.7771 23.32 1984.72 128.50 2136.53 2.039 3.093 1000 7.7146 27.43 2021.96 113.36 2162.75 1.799 2.729 Summary over the full sweep at T = 18 degC: pH : 8.3324 at 180 ppm -> 7.7041 at 1000 ppm (change -0.6283) [H+] : 4.6511e-09 -> 1.9765e-08 mol/kg (factor 4.250) [CO2*] : 6.03 -> 33.49 umol/kg (factor 5.556) [HCO3-] : 1593.39 -> 2083.05 umol/kg (factor 1.307) [CO3=] : 285.41 -> 87.80 umol/kg (factor 0.308) DIC : 1884.83 -> 2204.33 umol/kg (factor 1.170) Om_ar : 4.407 -> 1.356 (fraction remaining 0.308) Note the shape of it. Dissolved CO2 rises by the same factor as the atmosphere, 5.56. Bicarbonate barely moves, +30.7 percent. Carbonate is the species that pays: 69.2 percent of it is gone, and DIC has risen only 17.0 percent to absorb all of that carbon. ============================================================================== 4. BJERRUM FRACTIONS AND THE CROSSOVER POINTS ============================================================================== Fraction of DIC in each species as a function of pH, T = 18 degC, S = 35. These depend only on K1, K2 and [H+], not on alkalinity. pH CO2* HCO3- CO3= 5.00 0.890501 0.109490 0.000009 5.25 0.820562 0.179412 0.000027 5.50 0.719987 0.279939 0.000074 5.75 0.591107 0.408701 0.000191 6.00 0.448319 0.551222 0.000459 6.25 0.313510 0.685474 0.001016 6.50 0.204150 0.793759 0.002091 6.75 0.125841 0.870083 0.004076 7.00 0.074639 0.917715 0.007646 7.25 0.043125 0.942906 0.013969 7.50 0.024447 0.950512 0.025041 7.75 0.013628 0.942230 0.044142 8.00 0.007452 0.916218 0.076330 8.25 0.003968 0.867511 0.128521 8.50 0.002032 0.789875 0.208093 8.75 0.000984 0.680302 0.318714 9.00 0.000443 0.545281 0.454276 9.25 0.000184 0.402909 0.596907 9.50 0.000071 0.275121 0.724808 9.75 0.000025 0.175901 0.824074 10.00 0.000009 0.107168 0.892823 10.25 0.000003 0.063231 0.936766 10.50 0.000001 0.036570 0.963429 10.75 0.000000 0.020899 0.979101 11.00 0.000000 0.011861 0.988139 crossover CO2* = HCO3- at pH = pK1 = 5.9103 crossover HCO3- = CO3= at pH = pK2 = 9.0793 present-day surface pH = 8.0361, which sits 1.0432 units below pK2. At that pH the DIC split is CO2* 0.0068, HCO3- 0.9107, CO3= 0.0825. Ocean pH sits in the flat middle of the bicarbonate plateau, which is why bicarbonate hardly notices and carbonate does all the moving. ============================================================================== 5. THE COST OF NEGLECTING THE SECOND DISSOCIATION ============================================================================== Approximation A: solve the alkalinity balance with the 2[CO3=] term deleted, keeping borate and water. Carbonate is then back-calculated from the H that comes out. Approximation B: the schoolroom version, TA = [HCO3-] alone, which gives H = K1 [CO2*] / TA in closed form. xCO2 pH exact pH A err A pH B err B Om errA% Om errB% 180 8.3324 8.4594 +0.1270 8.4919 +0.1594 +79.46 +108.36 200 8.2972 8.4158 +0.1186 8.4461 +0.1489 +72.65 +98.48 240 8.2352 8.3401 +0.1049 8.3669 +0.1317 +62.14 +83.41 280 8.1817 8.2760 +0.0943 8.3000 +0.1183 +54.39 +72.43 320 8.1345 8.2202 +0.0857 8.2420 +0.1075 +48.41 +64.06 350 8.1024 8.1827 +0.0803 8.2031 +0.1007 +44.75 +58.97 380 8.0727 8.1482 +0.0756 8.1673 +0.0947 +41.62 +54.64 420 8.0361 8.1062 +0.0701 8.1239 +0.0877 +38.08 +49.79 450 8.0107 8.0772 +0.0665 8.0939 +0.0832 +35.81 +46.68 500 7.9715 8.0328 +0.0612 8.0482 +0.0766 +32.58 +42.30 560 7.9289 7.9849 +0.0560 7.9989 +0.0700 +29.42 +38.03 650 7.8722 7.9219 +0.0496 7.9342 +0.0620 +25.69 +33.04 750 7.8170 7.8611 +0.0441 7.8721 +0.0550 +22.53 +28.85 850 7.7682 7.8079 +0.0397 7.8177 +0.0495 +20.07 +25.60 1000 7.7041 7.7387 +0.0346 7.7471 +0.0430 +17.25 +21.92 Approximation A: mean |pH error| = 0.0739, worst = 0.1270 units mean |Omega error| = 41.66 %, worst = 79.46 % Approximation B: mean |pH error| = 0.0925, worst = 0.1594 units mean |Omega error| = 55.10 %, worst = 108.36 % Read the Omega columns twice. Approximation A shifts pH by a few hundredths of a unit and still misstates the aragonite saturation state by tens of percent, because Omega runs on [H+] squared and the neglected term was carrying exactly the quantity Omega needs. Error in [H+] from approximation A, percent, across the sweep: -25.35 -23.89 -21.47 -19.52 -17.91 -16.88 -15.97 -14.90 -14.19 -13.15 -12.10 -10.80 -9.66 -8.74 -7.65 And the same question asked the other way round. If you use approximation A to work out how much the hydrogen ion rises from 280 to 420 ppm, you get: exact solver : +39.80 % approximation A : +47.83 % (off by +8.02 points) approximation B : +50.00 % (off by +10.20 points) The ratio is less wrong than the absolute values, but it is still wrong by 8.0 and 10.2 percentage points. An approximation that misses the headline number by a fifth of its own size is not a shortcut. It is a different answer. ============================================================================== 6. WHERE THE ARAGONITE SATURATION STATE CROSSES ONE ============================================================================== Bisection on xCO2 for Omega_aragonite = 1, at fixed TA = 2300 umol/kg. Below the printed xCO2 the water is supersaturated with respect to aragonite; above it, aragonite is thermodynamically unstable. T degC S xCO2 at Om=1 pH there CO3= there 0.0 34.0 576.6 7.8902 66.18 5.0 34.0 740.1 7.8009 66.11 10.0 34.5 956.1 7.7074 65.85 15.0 35.0 1233.4 7.6141 65.33 18.0 35.0 1430.7 7.5602 64.77 20.0 35.0 1579.8 7.5243 64.33 25.0 35.0 2025.4 7.4339 63.02 30.0 35.0 2601.6 7.3429 61.42 Cold water crosses first, and by a long way. Carbon dioxide is more soluble in cold water and the dissociation constants shift with temperature, so a polar parcel at the same alkalinity runs out of carbonate at a far lower atmospheric mixing ratio than a tropical one. The 0 degC row is outside the calibration range of K1 and K2 and should be read as an extrapolation. A check worth making. Published model studies put surface aragonite undersaturation in the Southern Ocean at roughly 550 to 600 ppm. Our coldest row, 0.0 degC at S = 34.0, crosses at 576.6 ppm. That is an independent result our model did not see during construction, and it lands inside the published window. Carbonate concentration at the crossing, for reference: T = 0.0 degC : [CO3=] = 66.18 umol/kg at Omega = 1 T = 5.0 degC : [CO3=] = 66.11 umol/kg at Omega = 1 T = 10.0 degC : [CO3=] = 65.85 umol/kg at Omega = 1 T = 15.0 degC : [CO3=] = 65.33 umol/kg at Omega = 1 T = 18.0 degC : [CO3=] = 64.77 umol/kg at Omega = 1 T = 20.0 degC : [CO3=] = 64.33 umol/kg at Omega = 1 T = 25.0 degC : [CO3=] = 63.02 umol/kg at Omega = 1 T = 30.0 degC : [CO3=] = 61.42 umol/kg at Omega = 1 ============================================================================== 7. THE REVELLE FACTOR, COMPUTED NUMERICALLY ============================================================================== Revelle factor R = (d fCO2 / fCO2) / (d DIC / DIC) at constant TA, by central difference on a 0.1 percent perturbation of xCO2. xCO2 T degC DIC Revelle R pH 180 5.0 2010.12 10.008 8.3406 280 5.0 2093.01 12.030 8.1788 420 5.0 2158.56 14.404 8.0240 560 5.0 2199.55 16.206 7.9109 1000 5.0 2273.01 18.718 7.6763 180 18.0 1884.83 8.184 8.3324 280 18.0 1982.61 9.400 8.1817 420 18.0 2063.49 11.006 8.0361 560 18.0 2114.91 12.439 7.9289 1000 18.0 2204.33 15.690 7.7041 180 25.0 1813.34 7.550 8.3220 280 25.0 1916.92 8.445 8.1771 420 25.0 2005.12 9.682 8.0369 560 25.0 2062.31 10.837 7.9333 1000 25.0 2162.75 13.765 7.7146 The published range for present-day surface water is roughly 8 to 15, lowest in the warm subtropics and highest in cold high-latitude water. Our block runs 7.55 to 18.72, which is wider at both ends, and it should be: the block also contains glacial air at 180 ppm and an atmosphere at 1000 ppm that has never existed in the Quaternary. Restricting to the present-day column at 420 ppm gives 9.68 to 14.40, inside the published range. R = 11.01 at 420 ppm and 18 degC means a one percent rise in DIC buys a 11.0 percent rise in fCO2. The buffer gets worse as CO2 rises: R climbs from 8.18 at 180 ppm to 15.69 at 1000 ppm at this temperature, so each additional tonne is absorbed less willingly. ============================================================================== 8. MONTE CARLO: WHAT SURVIVES THE UNCERTAINTY IN THE CONSTANTS ============================================================================== Trials : 200000 Generator : PCG64 from SeedSequence(20250320).spawn(4)[0] Standard uncertainties applied, one independent normal draw each, held common between the 280 ppm and 420 ppm end members, because it is the same water and the same constants on both sides: u(pK0 ) = 0.0020 (log10 units) u(pK1 ) = 0.0075 (log10 units) u(pK2 ) = 0.0150 (log10 units) u(pKB ) = 0.0100 (log10 units) u(pKw ) = 0.0100 (log10 units) u(pKsp) = 0.0200 (log10 units) u(TA) = 2.0 umol/kg u(T) = 0.10 degC u(S) = 0.010 Results. SE is the standard error of the Monte Carlo mean, computed from the trials themselves as s.d./sqrt(N). rise in [H+], 280->420 ppm deterministic : 39.80324 % Monte Carlo mean : 39.80254 % standard error of mean 0.00043 % spread (1 s.d.) : 0.19032 % 95% interval : 39.42786 to 40.17320 % mean - deterministic : -0.00069 % (1.63 SE) pH at 420 ppm deterministic : 8.03614 Monte Carlo mean : 8.03611 standard error of mean 0.00002 spread (1 s.d.) : 0.00685 95% interval : 8.02271 to 8.04954 mean - deterministic : -0.00003 (1.96 SE) Omega aragonite at 420 ppm deterministic : 2.62705 Monte Carlo mean : 2.63073 standard error of mean 0.00032 spread (1 s.d.) : 0.14286 95% interval : 2.36165 to 2.92155 mean - deterministic : +0.00368 (11.53 SE) pH change 280->420 ppm deterministic : -0.14552 Monte Carlo mean : -0.14551 standard error of mean 0.00000 spread (1 s.d.) : 0.00059 95% interval : -0.14666 to -0.14435 mean - deterministic : +0.00000 (1.94 SE) [CO3=] at 420 ppm deterministic : 170.14891 umol/kg Monte Carlo mean : 170.19719 umol/kg standard error of mean 0.01095 umol/kg spread (1 s.d.) : 4.89581 umol/kg 95% interval : 160.74576 to 179.93591 umol/kg mean - deterministic : +0.04828 umol/kg (4.41 SE) The rise in [H+] is the durable number here: its 95 percent interval is 39.43 to 40.17 percent, a width of only 0.75 points, because the same constants act on both end members and most of the error cancels. The absolute pH is softer, spread 0.0068, and Omega is softest of all, spread 0.1429 on a value of 2.631, because pKsp enters Omega directly and nothing cancels it. 8a-ii. How far out is the 30 percent figure? Our value : 39.80254 % Quoted value : 30 % Gap : 9.80 percentage points In units of the Monte Carlo standard error of the mean : 23034 SE In units of the Monte Carlo spread (1 s.d.) : 51.5 s.d. Those two numbers say the same thing twice. Every published uncertainty in every constant, propagated jointly, moves the answer by about 0.19 percent. The gap to be explained is 9.8 percent. The constants cannot account for it and we are not going to pretend they can. The explanation is in section 2d-iii: the quoted figure describes an ocean that has not finished equilibrating, and our solver describes one that has. 8b. Convergence of the Monte Carlo mean trials running mean running SE |mean-det|/SE 100 39.78801 0.01951 0.78 200 39.80416 0.01360 0.07 500 39.80120 0.00850 0.24 1000 39.79980 0.00597 0.57 2000 39.80206 0.00428 0.27 5000 39.80222 0.00270 0.38 10000 39.80065 0.00191 1.36 20000 39.80013 0.00134 2.31 50000 39.80147 0.00085 2.08 100000 39.80171 0.00060 2.53 150000 39.80216 0.00049 2.20 200000 39.80254 0.00043 1.63 The standard error falls as 1/sqrt(N), as it must. From 100 trials to 200000 it drops by a factor of 45.9 against a predicted 44.7. ============================================================================== 9. FIGURE DATA (machine readable, for the article and the LaTeX) ============================================================================== FIG1 pH and Omega_aragonite against xCO2, three temperatures # T_degC xCO2 pH Omega_ar FIG1 5 180 8.3406 2.9843 FIG1 5 200 8.3027 2.7852 FIG1 5 240 8.2361 2.4596 FIG1 5 280 8.1788 2.2041 FIG1 5 320 8.1285 1.9977 FIG1 5 350 8.0943 1.8670 FIG1 5 380 8.0628 1.7527 FIG1 5 420 8.0240 1.6207 FIG1 5 450 7.9971 1.5342 FIG1 5 500 7.9558 1.4091 FIG1 5 560 7.9109 1.2837 FIG1 5 650 7.8514 1.1328 FIG1 5 750 7.7937 1.0021 FIG1 5 850 7.7429 0.8986 FIG1 5 1000 7.6763 0.7782 FIG1 18 180 8.3324 4.4066 FIG1 18 200 8.2972 4.1633 FIG1 18 240 8.2352 3.7544 FIG1 18 280 8.1817 3.4230 FIG1 18 320 8.1345 3.1481 FIG1 18 350 8.1024 2.9704 FIG1 18 380 8.0727 2.8124 FIG1 18 420 8.0361 2.6271 FIG1 18 450 8.0107 2.5038 FIG1 18 500 7.9715 2.3228 FIG1 18 560 7.9289 2.1381 FIG1 18 650 7.8722 1.9111 FIG1 18 750 7.8170 1.7102 FIG1 18 850 7.7682 1.5480 FIG1 18 1000 7.7041 1.3556 FIG1 25 180 8.3220 5.3090 FIG1 25 200 8.2882 5.0482 FIG1 25 240 8.2286 4.6038 FIG1 25 280 8.1771 4.2373 FIG1 25 320 8.1317 3.9289 FIG1 25 350 8.1008 3.7272 FIG1 25 380 8.0722 3.5464 FIG1 25 420 8.0369 3.3323 FIG1 25 450 8.0124 3.1887 FIG1 25 500 7.9745 2.9761 FIG1 25 560 7.9333 2.7567 FIG1 25 650 7.8783 2.4838 FIG1 25 750 7.8247 2.2390 FIG1 25 850 7.7771 2.0389 FIG1 25 1000 7.7146 1.7987 FIG2 speciation in umol/kg against xCO2 at T = 18 degC # xCO2 CO2star HCO3 CO3 DIC FIG2 180 6.03 1593.39 285.41 1884.83 FIG2 200 6.70 1632.55 269.65 1908.89 FIG2 240 8.04 1698.29 243.17 1949.49 FIG2 280 9.38 1751.53 221.70 1982.61 FIG2 320 10.72 1795.68 203.89 2010.29 FIG2 350 11.72 1824.20 192.39 2028.31 FIG2 380 12.72 1849.54 182.16 2044.42 FIG2 420 14.06 1879.28 170.15 2063.49 FIG2 450 15.07 1899.05 162.16 2076.28 FIG2 500 16.74 1928.07 150.44 2095.26 FIG2 560 18.75 1957.68 138.48 2114.91 FIG2 650 21.77 1994.05 123.78 2139.60 FIG2 750 25.11 2026.24 110.77 2162.12 FIG2 850 28.46 2052.23 100.26 2180.95 FIG2 1000 33.49 2083.05 87.80 2204.33 FIG3 Bjerrum fractions against pH at T = 18 degC, S = 35 # pH f_CO2 f_HCO3 f_CO3 FIG3 5.00 0.890501 0.109490 0.000009 FIG3 5.25 0.820562 0.179412 0.000027 FIG3 5.50 0.719987 0.279939 0.000074 FIG3 5.75 0.591107 0.408701 0.000191 FIG3 6.00 0.448319 0.551222 0.000459 FIG3 6.25 0.313510 0.685474 0.001016 FIG3 6.50 0.204150 0.793759 0.002091 FIG3 6.75 0.125841 0.870083 0.004076 FIG3 7.00 0.074639 0.917715 0.007646 FIG3 7.25 0.043125 0.942906 0.013969 FIG3 7.50 0.024447 0.950512 0.025041 FIG3 7.75 0.013628 0.942230 0.044142 FIG3 8.00 0.007452 0.916218 0.076330 FIG3 8.25 0.003968 0.867511 0.128521 FIG3 8.50 0.002032 0.789875 0.208093 FIG3 8.75 0.000984 0.680302 0.318714 FIG3 9.00 0.000443 0.545281 0.454276 FIG3 9.25 0.000184 0.402909 0.596907 FIG3 9.50 0.000071 0.275121 0.724808 FIG3 9.75 0.000025 0.175901 0.824074 FIG3 10.00 0.000009 0.107168 0.892823 FIG3 10.25 0.000003 0.063231 0.936766 FIG3 10.50 0.000001 0.036570 0.963429 FIG3 10.75 0.000000 0.020899 0.979101 FIG3 11.00 0.000000 0.011861 0.988139 FIG4 approximation error against xCO2 at T = 18 degC # xCO2 pH_exact pH_A pH_B Omega_exact Omega_A Omega_B FIG4 180 8.3324 8.4594 8.4919 4.4066 7.9080 9.1816 FIG4 200 8.2972 8.4158 8.4461 4.1633 7.1877 8.2634 FIG4 240 8.2352 8.3401 8.3669 3.7544 6.0875 6.8862 FIG4 280 8.1817 8.2760 8.3000 3.4230 5.2847 5.9024 FIG4 320 8.1345 8.2202 8.2420 3.1481 4.6721 5.1646 FIG4 350 8.1024 8.1827 8.2031 2.9704 4.2997 4.7220 FIG4 380 8.0727 8.1482 8.1673 2.8124 3.9830 4.3492 FIG4 420 8.0361 8.1062 8.1239 2.6271 3.6276 3.9350 FIG4 450 8.0107 8.0772 8.0939 2.5038 3.4004 3.6726 FIG4 500 7.9715 8.0328 8.0482 2.3228 3.0797 3.3054 FIG4 560 7.9289 7.9849 7.9989 2.1381 2.7671 2.9512 FIG4 650 7.8722 7.9219 7.9342 1.9111 2.4021 2.5426 FIG4 750 7.8170 7.8611 7.8721 1.7102 2.0955 2.2036 FIG4 850 7.7682 7.8079 7.8177 1.5480 1.8586 1.9443 FIG4 1000 7.7041 7.7387 7.7471 1.3556 1.5894 1.6527 FIG5 Monte Carlo convergence of the rise in [H+] (percent) # trials running_mean running_SE FIG5 100 39.78801 0.01951 FIG5 114 39.80055 0.01838 FIG5 129 39.80501 0.01703 FIG5 147 39.80880 0.01614 FIG5 167 39.80887 0.01471 FIG5 190 39.80376 0.01399 FIG5 217 39.79580 0.01321 FIG5 246 39.80502 0.01239 FIG5 280 39.80496 0.01180 FIG5 319 39.80193 0.01100 FIG5 363 39.80242 0.01036 FIG5 413 39.80495 0.00976 FIG5 469 39.80215 0.00889 FIG5 534 39.80199 0.00819 FIG5 607 39.80269 0.00769 FIG5 691 39.79877 0.00721 FIG5 786 39.80069 0.00682 FIG5 894 39.80259 0.00630 FIG5 1016 39.79887 0.00592 FIG5 1156 39.80245 0.00555 FIG5 1315 39.80314 0.00526 FIG5 1496 39.80332 0.00500 FIG5 1702 39.80203 0.00461 FIG5 1936 39.80169 0.00436 FIG5 2202 39.80205 0.00406 FIG5 2505 39.80345 0.00382 FIG5 2849 39.80353 0.00359 FIG5 3241 39.80285 0.00337 FIG5 3686 39.80295 0.00315 FIG5 4193 39.80302 0.00296 FIG5 4770 39.80314 0.00276 FIG5 5425 39.80106 0.00260 FIG5 6171 39.80114 0.00243 FIG5 7020 39.80059 0.00227 FIG5 7985 39.80031 0.00212 FIG5 9083 39.80152 0.00200 FIG5 10332 39.80096 0.00188 FIG5 11753 39.80054 0.00176 FIG5 13369 39.80065 0.00165 FIG5 15207 39.80059 0.00155 FIG5 17298 39.80051 0.00145 FIG5 19676 39.80055 0.00136 FIG5 22381 39.79992 0.00127 FIG5 25459 39.79931 0.00119 FIG5 28959 39.80003 0.00112 FIG5 32941 39.80071 0.00105 FIG5 37470 39.80074 0.00098 FIG5 42622 39.80107 0.00092 FIG5 48482 39.80131 0.00086 FIG5 55148 39.80177 0.00081 FIG5 62731 39.80210 0.00076 FIG5 71356 39.80213 0.00071 FIG5 81168 39.80189 0.00067 FIG5 92328 39.80182 0.00063 FIG5 105022 39.80179 0.00059 FIG5 119462 39.80201 0.00055 FIG5 135888 39.80197 0.00052 FIG5 154572 39.80212 0.00048 FIG5 175825 39.80212 0.00045 FIG5 200000 39.80254 0.00043 FIG6 Omega = 1 crossing against temperature # T_degC S xCO2_cross pH CO3_umol FIG6 0.0 34.0 576.6 7.8902 66.18 FIG6 5.0 34.0 740.1 7.8009 66.11 FIG6 10.0 34.5 956.1 7.7074 65.85 FIG6 15.0 35.0 1233.4 7.6141 65.33 FIG6 18.0 35.0 1430.7 7.5602 64.77 FIG6 20.0 35.0 1579.8 7.5243 64.33 FIG6 25.0 35.0 2025.4 7.4339 63.02 FIG6 30.0 35.0 2601.6 7.3429 61.42 FIG7 Revelle factor against xCO2 # xCO2 T_degC DIC RevelleR pH FIG7 180 5.0 2010.12 10.008 8.3406 FIG7 280 5.0 2093.01 12.030 8.1788 FIG7 420 5.0 2158.56 14.404 8.0240 FIG7 560 5.0 2199.55 16.206 7.9109 FIG7 1000 5.0 2273.01 18.718 7.6763 FIG7 180 18.0 1884.83 8.184 8.3324 FIG7 280 18.0 1982.61 9.400 8.1817 FIG7 420 18.0 2063.49 11.006 8.0361 FIG7 560 18.0 2114.91 12.439 7.9289 FIG7 1000 18.0 2204.33 15.690 7.7041 FIG7 180 25.0 1813.34 7.550 8.3220 FIG7 280 25.0 1916.92 8.445 8.1771 FIG7 420 25.0 2005.12 9.682 8.0369 FIG7 560 25.0 2062.31 10.837 7.9333 FIG7 1000 25.0 2162.75 13.765 7.7146 ============================================================================== 10. CLOSING NUMBERS ============================================================================== Reference condition : T = 18.0 degC, S = 35.0, TA = 2300 umol/kg pH, pre-industrial 280 ppm : 8.1817 (accepted about 8.17) pH, present day 420 ppm : 8.0361 (accepted about 8.05) pH change 280 -> 420 : -0.1455 units rise in [H+] 280 -> 420 : +39.80 % (quoted about 30 %) gap to the quoted figure : +9.80 points, 23034 SE and 51.5 s.d. out xCO2 that would give 30% : 385.2 ppm at full equilibrium xCO2 that reads pH 8.05 : 404.4 ppm at full equilibrium Monte Carlo on that rise : 39.80 % +/- 0.0004 % (SE), 95% 39.43 to 40.17 Omega aragonite at 280 ppm : 3.423 Omega aragonite at 420 ppm : 2.627 Omega aragonite at 1000 ppm: 1.356 [CO3=] at 280 / 420 / 1000 : 221.70 / 170.15 / 87.80 umol/kg Carbonate lost 180 -> 1000 : 69.2 % Max alkalinity residual : 4.337e-19 mol/kg (0.849 ulp of TA) Newton vs quartic agreement: 1.870e-14 relative Worst Omega error, approx A: 79.46 % Worst Omega error, approx B: 108.36 % Revelle factor at 420, 18C : 11.006 Master seed : 20250320 Wall clock : 1.33 s END OF OUTPUT