How Statisticians Found 15% Less Uncertainty in the Ocean's Hidden Heat

Roughly 90 percent of the excess heat trapped by rising greenhouse gases hasn't warmed the atmosphere. It hasn't melted ice or dried out soils. It has sunk, quietly and invisibly, into the world's oceans—accumulating over two decades into a reservoir of thermal energy that now dwarfs anything humanity has ever managed to put to practical use. This hidden heat shapes everything from sea level rise to hurricane intensity, from the pace of climate change itself to the survival of coral reefs and fisheries. And for years, scientists have struggled to measure it precisely enough to know whether their calculations of Earth's energy budget were trustworthy.

That may be changing. A new statistical framework, developed by researchers at Carnegie Mellon University and the University of Colorado Boulder, offers a more rigorous way to estimate how much heat the oceans have absorbed—and crucially, to quantify just how uncertain those estimates are. The approach doesn't rely on better instruments or more data. It relies on smarter math. By modeling the relationship between different depths in the ocean rather than treating them as separate problems, the team achieved up to a 15 percent reduction in uncertainty for global ocean heat content estimates. The work appears in a recent preprint by Thea Sukianto, Donata Giglio, and Mikael Kuusela, building on a univariate framework the group published earlier this year.

The 15 percent figure sounds modest. In isolation, it's hard to appreciate what it means. But uncertainty matters enormously in climate science, where small differences can determine whether a trend qualifies as statistically significant—whether scientists can confidently say the ocean is warming faster than expected, or whether a spike in heat content reflects a real shift or just noise in the data. "We've shown that modeling the vertical dependence reduces uncertainties on downstream quantity estimates by as much as 15 percent," the researchers write. "These estimated uncertainties are essential to analyze the statistical significance of OHC anomalies on both regional and global scales."

The Science

To understand what the researchers did, it helps to understand what they're measuring—and why it's so hard to measure.

Ocean heat content, or OHC, is the total amount of thermal energy stored in the water column. The basic formula involves integrating temperature measurements over depth and geographic area: you take the temperature at each point, multiply by the density of seawater and its heat capacity, then sum across the entire ocean basin. In practice, this means converting millions of individual temperature readings from Argo floats into gridded maps of heat content, then adding up all those grid cells to get a global total.

The problem is that the ocean is vast, observations are sparse, and the stakes are high. The Argo program, established in the early 2000s, has deployed roughly 4,100 autonomous profiling floats across the world's oceans. Each float descends to about 2,000 decibars of pressure—roughly 2,000 meters depth—collecting temperature and salinity profiles as it rises to the surface every ten days or so. Since the program reached near-global coverage around 2004, these floats have accumulated over three million vertical profiles. That's an extraordinary dataset, but it's still just a fraction of what would be needed to observe the ocean completely.

{ "image": 4 }

The Argo floats sample at roughly a 3-degree by 3-degree nominal resolution—grid cells about 330 kilometers on a side. Within each cell, you might have a handful of profiles per month. Between cells, especially in the deep ocean or near the poles, coverage can be spotty or nonexistent. And here's the complication that drives the new paper: not all profiles reach the same depth. Some floats don't descend as far due to technical limitations. Others are constrained by bathymetry—undersea mountains, continental shelves, and ridges that block deeper sampling. The result is that observations become increasingly sparse with depth.

For the upper ocean layer (15 to 975 decibars, roughly the top 1,000 meters), the researchers identified 1,540,593 quality-controlled profiles available for analysis. For the midocean layer (975 to 1,850 decibars, roughly 1,000 to 2,000 meters), that number drops to 1,128,932—a difference of about 27 percent. This isn't a trivial gap. When you try to map temperatures in the midocean using only midocean observations, you're working with fewer data points and noisier estimates. The midocean matters too: it accounts for roughly half of the ocean's heat storage capacity and plays a key role in long-term climate dynamics.

Previous approaches have dealt with this uneven sampling by partitioning the ocean into layers and mapping each one separately. The Roemmich-Gilson climatology, one of the most widely used Argo products, divides the water column into 58 pressure layers and interpolates each month independently. It works well for producing temperature fields, but it doesn't provide reliable uncertainty estimates—it's essentially a best-guess map without a measure of confidence.

Other methods have attempted to incorporate the vertical dimension directly into the statistical model, treating temperature or salinity as a continuous function of pressure rather than a set of discrete layers. Functional kriging approaches, for example, can predict temperature at any depth. But they run into computational walls when you try to extract uncertainties for global totals. "Obtaining an uncertainty on the global OHC anomaly would require the predictive covariance between every spatial location in the Argo-sampled global ocean," the researchers note, "which is not readily available from the local models."

The challenge, then, is twofold: produce accurate maps of ocean heat content at different depths, and attach reliable uncertainty estimates to those maps—especially when you add them together to get a global total. The Sukianto et al. framework tackles both problems using a technique called bivariate locally stationary Gaussian process regression.

{ "image": 5 }

The core idea is to model the upper and midocean layers jointly rather than separately, allowing information from the better-sampled upper layer to inform predictions in the midocean. If the two layers are correlated—and physical oceanography suggests they should be, since water masses at different depths interact through mixing and circulation—then knowing something about the upper ocean should tell you something about the deeper ocean. By modeling this correlation explicitly, the method can borrow strength across layers.

The approach builds on a univariate framework the same group published in 2026, which used locally stationary Gaussian processes to map ocean heat content within single pressure layers. "Locally stationary" means the statistical properties of the ocean vary depending on where you are—the correlation structure in the tropics differs from that in the subpolar gyres, and the method accounts for this by fitting separate models in local windows across the globe. The computational trick that makes this feasible is a white noise convolution representation, which allows the researchers to simulate from the predictive distributions without having to invert enormous covariance matrices.

The bivariate extension adds a second layer of complexity: instead of modeling one layer at a time, the method models both layers simultaneously, with a covariance structure that captures both the spatial dependence within each layer and the vertical dependence between them. "We jointly model the upper and midocean OHC anomalies," the researchers write. "Since the midocean layer has fewer observations, this improves the mapped anomaly fields by allowing the midocean to borrow predictive strength from the upper ocean."

To validate the approach, the team used Argo data from January 2004 through December 2022—a period of nearly two decades of measurements. They trained the model on quality-controlled profiles, mapped the heat content anomalies onto a 1-degree by 1-degree monthly grid, and then ran a series of experiments to see whether the bivariate approach outperformed the univariate one.

What They Found

The results are striking in their consistency. Across every test the researchers ran, the bivariate model outperformed the univariate model.

The primary validation involved cross-validation—a statistical technique where you hide a portion of the data, train the model on the rest, and then check how well the predictions match the hidden observations. If a model can accurately predict observations it hasn't seen, that's a sign it's capturing real patterns rather than overfitting to noise.

The researchers used two different cross-validation schemes. In leave-one-year-out validation (LOOO), they trained on all years except one, then predicted the held-out year. This tests how well the model generalizes across time. In leave-one-float-out validation (LOFO), they trained on all profiles except those from a single float, then predicted that float's observations. This tests how well the model generalizes across space.

In both schemes, for both the upper ocean and midocean layers, the bivariate model had lower root mean squared error than the univariate model. The improvements were modest but consistent—typically in the range of 5 to 15 percent depending on the layer and validation scheme. The midocean showed larger relative gains than the upper ocean, which makes sense: the midocean has fewer observations to begin with, so it has more to gain from borrowing information from the upper layer.

{ "image": 0 }

{ "image": 1 }

But the most consequential finding relates to uncertainty quantification. This is where the 15 percent figure comes in—and where the practical value of the bivariate approach becomes clearest.

When you add together the heat content of the upper ocean and midocean to get a global total, you need to account for the fact that these two estimates are not independent. They're correlated: if the upper ocean is warmer than expected, the midocean probably is too, and vice versa. If you ignore this correlation and simply add the uncertainties, you get a conservative (i.e., overestimate) the total uncertainty. The previous univariate framework produced just such an upper bound.

The bivariate approach allows the researchers to model this correlation explicitly, using the full covariance structure between the layers. The result is tighter confidence intervals for global ocean heat content. "Modeling this correlation results in improved OHC anomaly mapping and up to a 15 percent reduction of global OHC anomaly uncertainties in comparison to mapping the two layers separately without accounting for their dependence," they report.

{ "image": 6 }

{ "image": 7 }

This isn't just a statistical nicety. Uncertainty bounds matter for interpreting trends. If the uncertainty is too large, you can't tell whether a change in ocean heat content from one year to the next represents a real shift or just random variation. With tighter uncertainties, scientists can detect trends earlier and distinguish signal from noise more reliably.

The researchers demonstrated this with three climatological case studies using the 2004–2022 Argo data.

The first was the total ocean heat content trend—the long-term rate at which the ocean is gaining heat. Estimating this trend requires fitting a line through nearly two decades of monthly data and quantifying how precisely the slope is known. With the bivariate model, the confidence interval on the trend estimate is narrower than with the univariate model.

The second case study involved the total ocean heat uptake and its comparison with the top-of-atmosphere radiative flux. The Earth receives a certain amount of energy from the sun each year; some is reflected, some is absorbed by the atmosphere, and some is stored in the ocean. The ocean heat uptake is an indirect measure of Earth's energy imbalance—one of the most fundamental quantities in climate science. With improved uncertainty estimates, the researchers could better assess how closely the ocean heat content changes track independent satellite measurements of the energy budget.

The third case study was the cross-correlation between total ocean heat content anomalies and the Oceanic Niño Index (ONI)—a measure of El Niño and La Niña conditions in the tropical Pacific. This is perhaps the most demanding test of the method, because calculating a correlation between two time series requires knowing the uncertainty in both series and how they co-vary. "A rigorous uncertainty on the latter estimate would not have been possible to obtain without incorporating vertical dependence into the local conditional simulation realizations," the researchers note. "Without the bivariate model, you simply couldn't produce the kind of statistically principled confidence intervals that climate scientists need to test hypotheses."

The comparison between univariate and bivariate mapped anomaly fields reveals why. In the upper ocean, the two models produce similar patterns—the upper ocean is well-enough sampled that the additional information from the midocean doesn't dramatically change the predictions. But in the midocean, the bivariate model produces smoother, more coherent anomaly fields. Isolated patches of warmth or cool that appeared in the univariate model disappear in the bivariate model, because the midocean predictions are now informed by what the upper ocean is doing.

Why This Changes Things

Understanding the ocean's heat budget isn't an abstract scientific exercise. It's fundamental to understanding how fast the planet is warming, how much sea level will rise, and how extreme weather events will change in a warming world.

The ocean's role in climate is paradoxical. Water has an enormous capacity to store heat—far greater than air or land. This means the ocean acts as a buffer, absorbing most of the excess energy from greenhouse gases and moderating atmospheric warming. Without this thermal sink, surface temperatures would rise much faster than they have. But the ocean's capacity isn't infinite, and its uptake of heat has consequences of its own.

Warm water expands. As the ocean absorbs heat, it thermally expands, contributing to sea level rise. This "steric" expansion accounts for roughly half of global mean sea level rise over the past century, with the other half coming from melting ice sheets and glaciers. The IPCC's Sixth Assessment Report projects that sea level will continue rising for centuries, even if emissions are drastically reduced, because the ocean has already absorbed so much heat that some rise is now committed.

Heat also affects ocean circulation. The Atlantic Meridional Overturning Circulation (AMOC)—the conveyor belt of ocean currents that carries warm water northward and returns cold water southward—is partly driven by differences in water density. As the ocean warms and freshens (from melting ice), these density differences change, potentially weakening the circulation. A weaker AMOC would have profound effects on European climate, Atlantic hurricane activity, and global heat distribution. Scientists have been debating whether the AMOC has already weakened, but the uncertainties in ocean heat content estimates make it difficult to answer with confidence.

{ "chart": 0 }

Then there's the matter of marine ecosystems. Corals, in particular, are exquisitely sensitive to temperature. Even a degree or two of warming above historical norms can trigger bleaching events, where corals expel the symbiotic algae that give them color and nutrients. Severe or prolonged bleaching can kill corals outright. The Great Barrier Reef has experienced unprecedented mass bleaching events in recent years, and scientists attribute these directly to ocean warming. More precise estimates of ocean heat content help researchers understand when and where bleaching risks are highest.

Fish populations are also affected. Many commercially important species prefer specific temperature ranges and follow those ranges as oceans warm. This redistributes fish stocks, disrupting ecosystems and fisheries. Studies have linked warming oceans to declines in sardine populations off California, shifts in tuna distributions in the Pacific, and changes in recruitment success for cod in the North Atlantic. Better ocean heat estimates improve the models that forecast these changes.

Tropical cyclones offer another connection. Warmer ocean temperatures provide more energy to fuel storms, and research suggests that the most intense hurricanes are already becoming more frequent as the ocean warms. But it's not just the surface temperature that matters—the temperature of the upper ocean layer determines how much heat is available to a storm as it churns the water column. Warm core eddies and ocean heat islands can intensify hurricane rapid intensification, the process by which weak storms suddenly become major hurricanes within hours. More accurate ocean heat mapping could improve hurricane forecasting.

The researchers' work matters precisely because it improves the reliability of ocean heat estimates. "Recent intercomparison efforts have showed a wide spread in OHC estimates based on varying data analysis pipelines and modeling assumptions," they note, citing a 2024 study by Hakuba et al. and a 2023 study by von Schuckmann et al. Different research groups, using different methods, arrive at somewhat different values for how much heat the ocean has accumulated. This spread creates uncertainty not just in the estimates themselves, but in everything that depends on them—including Earth's energy imbalance, sea level projections, and climate sensitivity estimates.

Some of this spread comes from real differences in how groups process Argo data—their quality control procedures, their mapping methods, their choices about which regions to include. But some of it comes from insufficient accounting for uncertainty. If one group reports an ocean heat content estimate without a rigorous confidence interval, and another reports an estimate with an overly conservative interval, the spread between them partly reflects methodological choices rather than true differences in the ocean.

The bivariate framework helps address this by providing a statistically rigorous way to quantify uncertainty. The confidence intervals it produces aren't ad hoc or conservative—they're derived from the data and the model structure. When the researchers sum the upper ocean and midocean heat content to get a global total, their uncertainty estimate properly accounts for the covariance between the layers. This makes the estimate more honest: if there's substantial correlation between the layers, the total uncertainty is smaller than you might expect from simply adding the individual uncertainties.

This has practical implications for how scientists interpret trends. Consider a question like: is the ocean warming faster in recent years than in earlier decades? To answer this, you'd fit separate trends to different time periods and ask whether the difference is statistically significant—that is, larger than you'd expect from random variation alone. If your uncertainty estimates are too large, you might conclude that an apparent acceleration is just noise. With tighter, more accurate uncertainties, you have more power to detect real changes.

The researchers demonstrate this with their three case studies. The ONI correlation is particularly illustrative. El Niño and La Niña events redistribute heat across the Pacific, affecting atmospheric circulation, rainfall patterns, and temperatures worldwide. Understanding how ocean heat content relates to these events helps scientists attribute observed changes to the right causes. But quantifying this relationship requires knowing not just the correlation point estimate, but its uncertainty—which is difficult to obtain without the kind of rigorous simulation-based inference that the bivariate model enables.

From a computational perspective, the method is also notable for its feasibility. "Most of the floats in the Argo array descend to approximately 20,000 dbar and collect T/S profiles roughly every 10 days as they ascend the water column to the surface," the researchers note. With over three million profiles collected since the early 2000s, producing estimates using traditional spatio-temporal methods would be computationally prohibitive. The white noise convolution trick in the locally stationary Gaussian process framework makes the calculations tractable on standard hardware, while the bivariate extension doesn't substantially increase the computational burden.

The modular, open-source implementation matters too. The researchers have released their code under a permissive license at GitHub (github.com/ttsukianto/LocalGP_OHC), allowing other groups to apply the method to different problems or extend it in new directions. "While we demonstrate our improved framework by jointly mapping OHC anomaly fields in two pressure layers," they write, "this work will facilitate further oceanographic applications with rigorously quantified uncertainties such as joint mapping of OHC and sea surface temperature (SST) or sea surface height (SSH)."

Sea surface height is particularly interesting. As the ocean warms, it expands; as it gets saltier, it contracts. Measuring sea surface height from satellites therefore provides another window into ocean density and heat content changes. But connecting satellite SSH measurements to subsurface ocean properties requires models that account for how the ocean stratifies—and that means models that can jointly map temperature, salinity, and their vertical structure. The bivariate framework developed here provides a template for how to do this with proper uncertainty quantification.

What's Next

The paper leaves several questions open, and the researchers are candid about them.

First, the analysis considers just two pressure layers. The ocean has more structure than that—above 1,000 meters, there's the warm surface mixed layer and the seasonal thermocline; below, there's the permanent thermocline and the deep ocean. A more complete treatment would map a dozen or more layers, or better yet, treat pressure as a continuous variable as in functional kriging approaches. The challenge is that each additional layer adds computational and modeling complexity, and the researchers don't yet know whether the gains from modeling more layers would justify the added complexity.

Second, the method assumes stationarity within local windows but allows the covariance parameters to vary across space. This is a reasonable approximation, but it's not perfect. The ocean has strong regional variations in its dynamics—the Gulf Stream behaves very differently from the subtropical gyres, which behave very differently from the Southern Ocean. A global model with spatially varying parameters is a compromise between computational tractability and physical realism. Future work might explore whether more flexible nonstationary models, or models that incorporate physical constraints like known circulation patterns, could improve predictions further.

Third, the cross-validation experiments show consistent improvements from the bivariate model, but they're based on in-sample residual analysis. Testing the uncertainty estimates themselves—whether the confidence intervals actually contain the true values 95 percent of the time—would require external validation against independent data or simulation studies. The researchers acknowledge this limitation and suggest it as a direction for future work.

Fourth, the Argo float network itself has limitations. Coverage is sparse in some regions, particularly the Arctic Ocean (where Argo floats struggle under sea ice), the ocean surrounding Antarctica, and shallow coastal zones. The method's performance in these poorly-sampled regions is less certain, and improving Argo coverage there remains an operational challenge.

Looking further ahead, the framework opens several avenues for oceanographic research. Joint mapping of ocean heat content and sea surface temperature could help bridge satellite and in-situ observations, improving the consistency of climate data records. Joint mapping of ocean heat content and sea surface height could improve estimates of steric sea level rise and help attribute observed sea level changes to thermal expansion versus ice melt. Extensions to salinity—mapping ocean freshening from ice melt alongside warming—could illuminate changes in ocean circulation.

The researchers also point to applications in understanding Earth's energy imbalance more directly. The energy imbalance at the top of the atmosphere—the difference between incoming solar radiation and outgoing thermal radiation—is a small number (order 1 watt per square meter) compared to the fluxes themselves (roughly 340 watts per square meter of solar input). Detecting changes in this imbalance, and attributing them to natural variability versus human forcing, requires precise accounting of where excess energy is going. The ocean is the dominant sink, so better OHC estimates directly translate to better EEI estimates.

For the broader climate science community, the paper offers a template for how to do uncertainty quantification rigorously without sacrificing computational feasibility. The gap between what statisticians know how to do in principle and what can actually be computed on large datasets is often large. The white noise convolution approach, combined with local stationarity assumptions, is a practical solution that others could adapt to different problems.

The ocean has absorbed more than a thousand times humanity's annual energy consumption over the past century. That heat is invisible from the surface, but its consequences—rising seas, shifting currents, intensifying storms—are becoming increasingly visible. Better measurements of ocean heat content won't stop the warming, but they'll help society understand what we're already committed to, and what choices remain.

{ "chart": 1 }