Solving Seawater: A Carbonate Equilibrium Solver and What It Says About Acidification

Four Equations and One Unknown

Seawater holds its carbon in three forms, and they are in conversation with each other and with the air. There is dissolved carbon dioxide, there is bicarbonate, and there is carbonate, which is the one this article is about. Add a proton to carbonate and you get bicarbonate; add another and you get carbonic acid, which falls apart into water and dissolved gas. Each of those steps has an equilibrium constant, water itself has one, and holding the whole arrangement in place is a conservation law called total alkalinity, which is a statement about charge that does not change when you dissolve a neutral gas into the ocean.

Write those down together and something nice happens. Every concentration becomes a function of the single unknown \([\mathrm{H}^+]\). Substitute them all into the alkalinity definition, clear the denominators, and four coupled equilibria fall out as one quartic polynomial with exactly one positive root. The system is over before it starts. There is nothing to iterate between and nothing to relax, only a root to find on an interval where the function is strictly decreasing.

We wrote the solver because we wanted to see that happen. Then we ran it across the range of atmospheric carbon dioxide that the Earth has had and is going to have, from the 180 parts per million of the last glacial maximum to a thousand, and the elegance stopped being the interesting part.

The third of those numbers is the one that keeps us awake. Carbonate ion is the raw material every shell and every coral skeleton is assembled from. It is also the minority species, about eight per cent of the dissolved inorganic carbon at present-day pH, and it is the species that absorbs almost the entire chemical shock of adding carbon dioxide to the ocean. Bicarbonate barely notices. Carbonate pays.

What This Is, and What It Is Not

The club has no research vessel. We have no Niskin bottles, no coulometer, no spectrophotometric pH cell, and no access to the ocean beyond a train ride. What we have is a room with computers in it, so our research is computation. Every number in this article is the printed output of a Python file we wrote, applied to equilibrium constants that other people measured in laboratories over the past fifty years.

That distinction matters more here than it usually does, because the output of a carbonate solver looks exactly like the output of a measurement. A pH of 8.0361 reads as though somebody put an electrode in the water. Nobody did. It is what the thermodynamics say the water would do if it sat in equilibrium with 420 ppm of air at 18 degrees and salinity 35 for long enough, with no biology in it and nothing mixing it. When we write that a quantity was "measured" in this article, we mean measured from the solver.

What is genuinely empirical here belongs to other people. Weiss measured the solubility of carbon dioxide in seawater [1]. Mehrbach and colleagues measured the apparent dissociation constants of carbonic acid, in real seawater rather than an artificial substitute, across 2 to 35 degrees and salinity 19 to 43 [2], and Lueker, Dickson and Keeling refit those measurements onto the total hydrogen ion scale and validated the refit against gas-solution equilibrations [3]. Dickson did the boric acid [4]. Millero assembled the water dissociation [5]. Mucci measured the solubility of calcite and aragonite [6]. The boron and calcium contents of seawater come from Uppström's chlorinity ratio [7] and from Riley and Tongudai's major-cation ratios [8], and the vapour pressure over seawater from Weiss and Price [9]. Our contribution is arithmetic on top of their work, and an account of what that arithmetic says.

The System, Written Out

Four equilibria. Carbon dioxide crossing the sea surface and dissolving, governed by Henry's law with a solubility constant \(K_0\) defined on fugacity rather than partial pressure:

Then the two dissociations of carbonic acid, which between them do all the work:

Then water, \(K_w = [\mathrm{H}^+][\mathrm{OH}^-]\), and boric acid, which is not part of the carbonate system at all but contributes a tenth of the alkalinity and cannot be left out. Total alkalinity is defined as the excess of proton acceptors over proton donors relative to a zero level of protons, and in surface seawater it is, to better than half a per cent,

Every term on the right is a known function of \([\mathrm{H}^+]\) once the constants and \(f\mathrm{CO_2}\) are fixed. Write \(C = K_0 f\mathrm{CO_2}\) and \(h = [\mathrm{H}^+]\), then the equation to solve is

Every term in \(F'(h)\) is negative. The function is strictly decreasing on \((0, \infty)\), which means the root is unique and any bracketing method will find it. Multiply through by \(h^2(K_B + h)\) and the same statement becomes a quartic:

Two forms of the same statement. We solve the first by safeguarded Newton iteration, which takes a Newton step when it lands inside the current bracket and a geometric bisection when it does not, and we solve the second by handing the coefficients to an eigenvalue routine. They are completely different algorithms and they have to agree. Section 7 reports by how much.

The aragonite saturation state falls out afterwards as a ratio, not as a solve:

Calcium is conservative in seawater, so it tracks salinity and does essentially nothing. The numerator moves because carbonate moves. \(\Omega\) above one means aragonite is thermodynamically stable in that water; below one it dissolves.

Working Notes From the Club Table

The Arithmetic

Here is the whole calculation at the reference condition, by hand, so that nothing in the rest of the article has to be taken on trust. Temperature 18 °C, salinity 35, total alkalinity 2300 µmol per kilogram of seawater. The constants the solver prints at that condition:

Now do it again at 280 ppm, which gives \(h = 6.581793 \times 10^{-9}\) and pH 8.1817. The ratio of the two hydrogen ion concentrations is 9.201560/6.581793 = 1.398032. Subtract one, multiply by a hundred: 39.80%. That is the number, and it is arithmetic anybody can repeat with a calculator once they have the two roots.

The pH difference is 0.1455 units. Those two statements are the same statement, because pH is a logarithm, and \(10^{0.1455} = 1.3980\). Any time you see a pH decline quoted, you can convert it yourself: a decline of \(x\) units is a \((10^x - 1)\times 100\) per cent rise in hydrogen ion concentration. A decline of 0.1 units is 25.9%. A decline of 0.12 is 31.8%.

Where the Carbonate Goes

Figure 1 is the whole result in two panels. Left, pH against atmospheric carbon dioxide on a log axis, at three temperatures. Right, the aragonite saturation state on the same axis, with the \(\Omega = 1\) line drawn across it.

The pH curves are almost straight on a log axis and almost parallel, which is a real feature rather than an artefact of the plotting: over this range the buffer system is nearly logarithmic in carbon dioxide and nearly temperature-independent in slope. The saturation curves are not parallel at all. Cold water starts lower and crosses one sooner. A 5 °C parcel at alkalinity 2300 crosses \(\Omega_{ar} = 1\) at 740 ppm; an 18 °C parcel crosses at 1431 ppm; a 30 °C parcel not until 2602 ppm. Push the cold parcel down to 0 °C and it crosses at 577 ppm, which is a number nobody should be comfortable with.

Figure 2 shows where the carbon actually goes, and it is the figure we would put on the wall.

Read the right panel first. Between 180 and 1000 ppm, dissolved carbon dioxide rises by a factor of 5.556, which is very nearly the factor the atmosphere rose by, exactly as Henry's law requires. Bicarbonate rises 30.7%. Carbonate falls by 69.2%, from 285.41 to 87.80 µmol/kg. Total dissolved inorganic carbon rises only 17.0%, because the ocean is a buffer and that is what buffers do.

The mechanism is visible in the alkalinity balance. Alkalinity is fixed. Adding carbon dioxide adds protons. Those protons have to be absorbed by something, and the only thing in the balance with a coefficient of 2 is carbonate. Each proton converts one carbonate ion to one bicarbonate ion, which keeps alkalinity constant while raising dissolved inorganic carbon and lowering carbonate. The ocean buffers its pH by spending its carbonate, and it has a finite amount.

The full grid

The full grid, three temperatures by fifteen atmospheres, all at salinity 35 and alkalinity 2300 µmol/kg. Rows where the aragonite saturation state has fallen below one are marked. Every figure is printed exactly as the solver produced it.

The atmosphere at which \(\Omega_{ar}\) reaches one, by temperature, at this alkalinity: 577 ppm at 0 °C and salinity 34, 740 at 5 °C, 956 at 10 °C, 1233 at 15 °C, 1431 at 18 °C, 1580 at 20 °C, 2025 at 25 °C, 2602 at 30 °C. The carbonate concentration at the crossing is remarkably stable across that whole range, 66.18 µmol/kg at the cold end and 61.42 at the warm end, because the solubility product and the calcium concentration both move only slightly. It is the carbon dioxide needed to get there that changes.

The Bjerrum plot behind all of it

Strip out alkalinity entirely and plot the fraction of dissolved inorganic carbon in each species against pH. That plot, named after Niels Bjerrum, depends on nothing but \(K_1\), \(K_2\) and the hydrogen ion concentration, and it is the single most useful picture in aqueous carbonate chemistry.

At 18 °C and salinity 35 our solver puts pK₁ at 5.9103 and pK₂ at 9.0793. Those are the crossover points: below pH 5.91 dissolved carbon dioxide dominates, between 5.91 and 9.08 bicarbonate dominates, above 9.08 carbonate dominates. Surface seawater sits at 8.04, which is 1.0432 units below pK₂, on the shoulder of the bicarbonate plateau where the bicarbonate fraction is flat and the carbonate fraction is steep.

At present-day pH the split is 0.0068 dissolved carbon dioxide, 0.9107 bicarbonate, 0.0825 carbonate. Now nudge pH down by 0.1455 units, the change from 280 to 420 ppm. The bicarbonate fraction is near its maximum, so its derivative is near zero and it barely responds. The carbonate fraction is on a steep part of its curve and it falls hard. This is the whole mechanism described above, and the Bjerrum plot shows it with no alkalinity bookkeeping at all.

It also explains why the ocean's buffering degrades. The Revelle factor, the ratio of the fractional change in carbon dioxide fugacity to the fractional change in dissolved inorganic carbon at fixed alkalinity [10], measures how reluctantly the ocean takes up another tonne. We compute it by central difference on a 0.1% perturbation. At 18 °C it rises from 8.184 at 180 ppm to 11.006 at 420 ppm to 15.690 at 1000 ppm. The present-day column across our three temperatures runs 9.682 to 14.404, inside the published range of roughly 8 to 15 [14]. As carbonate is spent, the buffer weakens, and the ocean absorbs a smaller share of what we emit.

The Number That Did Not Come Out Right

We expected 30%. Everyone quotes 30%. We got 39.80%, and the honest thing to do with a disagreement of that size is to publish the debugging rather than the rounding.

Step one was the solver. Three algorithms, sharing no code: the safeguarded Newton, a log-space bisection that took 54 halvings to run out of floating-point room, and the roots of the cleared quartic. Newton and bisection returned bit-identical values. Newton and the quartic agreed to 1.870 × 10⁻¹⁴ relative, which is about a hundred ulp and is what you expect from a quartic root-finder working on coefficients that span twenty orders of magnitude. The alkalinity residual at the solution, evaluated at all 64 points of a temperature-by-carbon-dioxide check grid, never exceeded 4.34 × 10⁻¹⁹ mol/kg, which is 0.849 ulp of the alkalinity it is checking. The solver is not wrong.

Step two was the constants. Our calcite solubility product reads 4.2724 × 10⁻⁷ at 25 °C and salinity 35, against the value from Mucci's paper of about 4.27 × 10⁻⁷. Our pK_B reads 8.5975 at the same condition, which is Dickson's published figure to four decimal places. Our pK₁ and pK₂ reproduce the Lueker coefficients as printed. The constants are not wrong either.

Step three was to ask what change does give 30%. Hold the 280 ppm baseline, bisect on the other end point, and the answer is 385.2 ppm. Ask instead which atmosphere produces the classic decline of exactly 0.1 pH units: 370.7 ppm, giving a rise of 25.89%. Ask which atmosphere produces a reading of pH 8.05: 404.4 ppm.

Which is when somebody looked up when the claim was written.

The "0.1 pH units, about 30 per cent" statement entered circulation in the middle 2000s, off the back of Caldeira and Wickett [11] and Feely and colleagues [12]. Atmospheric carbon dioxide then was near 370 ppm. It is now 420. The claim was correct when it was made and has been repeated unchanged for twenty years while the quantity it describes kept moving.

So we asked the solver the question the claim was originally answering. Jiang and colleagues, combining the Surface Ocean CO₂ Atlas with Earth system model output, put the global surface decline from 1770 to 2000 at 0.11 ± 0.03 pH units [14]. Our solver, run from 280 to 370 ppm at the reference condition, returns 0.0993 units. That is 0.36 standard deviations from the published estimate, which is agreement.

There is a second effect underneath the first, smaller and harder to pin down. Our parcel is in instantaneous equilibrium with the air above it. A real mixed layer takes months to a year to equilibrate through gas exchange, and the atmosphere has been rising the entire time, so observed surface pH tracks an effective mixing ratio somewhat below the contemporaneous atmospheric one. Our inversion puts the effective value behind a reading of pH 8.05 at 404.4 ppm rather than 420. We cannot separate the two effects with this model, and we say so.

The conclusion we draw is narrow and we want to keep it narrow. We are not claiming the ocean has acidified more than the literature says. We are claiming that a specific round number, computed for a 370 ppm atmosphere, is still being quoted for a 420 ppm one, and that a solver asked about 420 ppm returns 39.80%. If you want the number for today's atmosphere at full equilibrium, that is it.

What the Shortcut Costs

Textbooks that introduce the carbonate system very often drop the second dissociation. The reasoning is reasonable on its face: carbonate is only eight per cent of the dissolved carbon, the algebra without it is a cubic instead of a quartic, and in the crudest version the whole thing collapses to \(h = K_1 C / A_T\), which is one line and no iteration at all.

We implemented two versions of the shortcut. Approximation A deletes the \(2[\mathrm{CO_3^{2-}}]\) term from the alkalinity balance and keeps borate and water. Approximation B is the schoolroom version, \(A_T = [\mathrm{HCO_3^-}]\) and nothing else. In both cases we then back-calculate carbonate and \(\Omega\) from the hydrogen ion concentration that comes out, which is what anybody using the shortcut would do.

Approximation A is wrong by 0.0739 pH units on average and by 0.1270 units at worst. Approximation B is wrong by 0.0925 on average and 0.1594 at worst. Put those beside the entire acidification signal we have been discussing, which is 0.1455 units, and the shortcut's error is the same size as the phenomenon. A student who computes pre-industrial pH with approximation A and present-day pH with the exact solver would conclude that the ocean had become more alkaline.

The saturation state is worse, and it is worse for a reason worth stating. \(\Omega\) depends on \([\mathrm{H}^+]^{-2}\), so a fractional error in \(h\) is doubled. Approximation A misstates \(\Omega_{ar}\) by 41.66% on average and by 79.46% at the glacial end. Approximation B reaches 108.36%, which is to say it more than doubles the saturation state. An approximation that halves the acidity and doubles the carbonate availability is not a simplification of the problem. It is a different problem with a reassuring answer.

The ratio survives better than the absolute values, which is a genuine and slightly awkward finding: approximation A gives a rise in hydrogen ion of 47.83% against the exact 39.80%, and approximation B gives 50.00%. Still wrong by eight and ten percentage points respectively, but far less wrong than the pH values that produced them, because the two errors partly cancel in the ratio.

Two Hundred Thousand Oceans

The constants are measurements and measurements have uncertainties. We ran 200,000 trials in which every equilibrium constant was perturbed by an independent normal draw with the standard uncertainty recommended by Orr and colleagues [16] and implemented as the default in the seacarb package: 0.0020 in pK₀, 0.0075 in pK₁, 0.0150 in pK₂, 0.0100 in pK_B and pK_w, 0.0200 in pK_sp, all in log₁₀ units, plus 2 µmol/kg on alkalinity, 0.1 °C on temperature and 0.010 on salinity. Each trial's draw was held common between the 280 ppm and 420 ppm end members, because it is the same water and the same constants on both sides and pretending otherwise would inflate the spread dishonestly.

The rise in hydrogen ion concentration comes out at 39.80254% with a standard error of 0.00043%, computed from the trials themselves as \(s/\sqrt{N}\). The trial-to-trial spread is 0.19032% and the 95% interval runs 39.42786% to 40.17320%. The standard error fell by a factor of 45.9 between 100 trials and 200,000, against the \(\sqrt{N}\) prediction of 44.7.

The reason that interval is so narrow is the common draw. Most of the error in the constants shifts both end members in the same direction and cancels in the ratio. The quantities where nothing cancels are much softer. The absolute pH at 420 ppm has a spread of 0.0068 units, with a 95% interval of 8.02271 to 8.04954. The aragonite saturation state has a spread of 0.14286 on a mean of 2.63073, a 95% interval of 2.36165 to 2.92155, because pK_sp enters \(\Omega\) directly and has the largest stated uncertainty of any constant in the set.

That matters for how the \(\Omega = 1\) crossings in §6 should be read. A five per cent uncertainty in \(\Omega\) translates into a real uncertainty in where the crossing sits, and our crossing values are printed to one decimal place because the solver is deterministic, not because the ocean is.

The gap to the 30% figure, in these units, is 51.5 trial-to-trial standard deviations and about 23,000 standard errors. The constants cannot close it, which is exactly what §7 concluded by a different route.

The Strongest Objection We Could Make

Here is the objection, stated as hard as we can state it.

This study contains no new information about the ocean. It is a reimplementation of CO2SYS, a program that has existed since 1998 and has since been rebuilt in ten languages and cross-validated between them to 0.0002 pH units [15]. Everything our solver does, seacarb and mocsy and PyCO2SYS already did, better, with the minor acid systems included and the pressure corrections applied. A club that wanted the numbers could have installed a package in thirty seconds. The only quantity we produced that is not already in the literature is 39.80%, and that is 39.80% for a model condition we chose, at a temperature and alkalinity we selected to be representative, in a state of equilibrium the surface ocean is never in. Change the reference condition and the number changes: our own sensitivity block spans 37.25% to 44.00% across a plausible range of surface waters. Publishing a two-decimal figure that moves by seven percentage points when you change the temperature is not a measurement. It is a demonstration.

We think most of that is correct, and we want to be precise about which parts we accept.

We accept that the solver is a reimplementation and we say so in §3. The value in writing it ourselves is pedagogical and verificational, not novel. A club that installs a package gets the right answer and learns nothing about why the quartic has one positive root. A club that writes the solver and checks it three ways against itself and against published constants can then read the packages' documentation and understand what it is reading. We would make the same argument about writing a least-squares routine once before using one.

We accept that 39.80% is condition-dependent and we print the dependence rather than hiding it. The sensitivity block is in the output and its range is in the abstract of this section. What we do not accept is that condition-dependence makes the number meaningless. The quoted 30% is condition-dependent in exactly the same way and is never presented with a range at all. Our contribution is to show the range, and to show that the centre of it at 420 ppm is nowhere near 30.

The part of the objection we cannot answer is the equilibrium assumption. It is the largest approximation in the study by a distance, it is not quantifiable within the model, and it pushes the answer in a known direction. A solver that granted instantaneous equilibrium will always report more acidification than an ocean that is still catching up. We have said where that matters and we cannot say by how much.

What a Shell Has To Do

An oyster larva has about forty-eight hours. In that window it has to build a shell from nothing, out of aragonite, using carbonate ions pulled from the water around it, and if it cannot finish it does not get a second attempt. The energy cost of that construction depends on how far the surrounding water is from saturation, because the further below saturation the water sits, the harder the organism has to work to concentrate carbonate at the site of calcification against the gradient.

Our solver says the surface ocean at 18 °C has gone from \(\Omega_{ar} = 3.423\) at 280 ppm to 2.627 at 420 ppm. Still supersaturated. Still, thermodynamically, a place where aragonite is stable. But the carbonate concentration has fallen from 221.70 to 170.15 µmol/kg, which is 23% of the raw material gone in two centuries, and at 1000 ppm it is 87.80 µmol/kg, which is 69% gone. Cold water is already much closer to the line, and our coldest parcel crosses it at 577 ppm.

The meta-analytic evidence on what this does to organisms is mixed by design, because the organisms are not alike. Kroeker and colleagues synthesised 228 studies and found reduced survival, calcification, growth and development when the whole range of marine species is pooled, with enormous variation between taxa and between studies [17]. Some calcifiers are barely affected. Some are affected badly. Our model cannot distinguish between them and does not try.

What the model can say is the part that does not depend on biology at all. The carbonate ion concentration in surface seawater is falling, the arithmetic of why is four equations long, and the quantity that is falling is the one that shells are made of.

Where a different modelling choice would change the answer

Equilibrium against disequilibrium. This is the big one and it is worth repeating. Granting our parcel instantaneous equilibrium with the atmosphere overstates present-day acidification relative to an ocean whose mixed layer is still catching up. Our inversion puts the effective mixing ratio behind an observed pH of 8.05 at 404 ppm rather than 420, which is a hint at the size of the effect and not a measurement of it. A model with a gas-exchange timescale in it would report a smaller number than 39.80%.

Fixed alkalinity against responsive alkalinity. We hold total alkalinity constant while adding carbon dioxide, which is correct bookkeeping for dissolving a neutral gas. Over centuries it stops being correct, because alkalinity is moved by carbonate sediment dissolution, by river input, by changes in biological calcification and by the weathering of rock on land. Dissolving carbonate sediment raises alkalinity and partly restores the carbonate ion, which is the mechanism that eventually neutralises fossil carbon over tens of thousands of years [12]. None of that is in our model, and including it would make the long-term picture less severe and the short-term picture unchanged.

Which set of constants. We used the Mehrbach measurements as refit by Lueker, Dickson and Keeling [2][3], which is the recommended set for seawater in the range we work in. Other published sets exist and differ. Lueker's own paper reports that formulations based on other studies gave calculated fugacities about 10% lower than their equilibrations. A 10% error in fugacity is roughly a 0.04 pH unit shift, which is a third of our entire acidification signal. This is the choice most likely to move our absolute pH values, and it is the reason the solver checks in §7 are worth the space they take.

Surface against depth. Every constant here is at one atmosphere. The solubility products of both carbonate minerals increase substantially with pressure, which is why the saturation horizon exists and sits at a few hundred to a few thousand metres depending on basin [12]. A pressure-corrected version of this solver would produce a completely different and considerably more alarming figure for the deep ocean. We did not write one, and no number in this article applies below the surface.

Total scale against seawater scale. Our K₁, K₂ and K_B are on the total hydrogen ion scale. Millero's K_w as published is on the seawater scale and we do not convert it. The conversion would move pK_w by under 0.001, and since hydroxide is 3.51 µmol/kg in a 2300 µmol/kg balance, the effect on alkalinity is about 0.008 µmol/kg. That one genuinely does not matter, and we checked rather than assuming.

Minor acid systems. Phosphate, silicate, ammonia, sulfide and the fluoride and bisulfate protons are all absent from our alkalinity definition. In open-ocean surface water they are together well under half a per cent of alkalinity. In a coastal estuary, in an upwelling zone, or anywhere with organic alkalinity, they are not, and this solver should not be pointed at such water.

Reproducing this

One file, one command, no packages beyond numpy.

python carbonate-solver.py > carbonate-solver-output.txt

Expected runtime is a couple of seconds. Ours was 1.33 s on a laptop under Python 3.12.3 and numpy 2.4.2. The dominant cost is the 200,000-trial Monte Carlo, which is fully vectorised; the deterministic grids are effectively free. The master seed is 20250320, hard-coded at the top of the file, and the Monte Carlo draws its generator from it through numpy's SeedSequence.spawn, so the entire output is deterministic and byte-reproducible. If your numbers differ from ours by more than the printed standard errors, we would like to hear about it. The interactive model runs the same solver in your browser and reproduces the headline values at its default settings.

References