Ocean heat content responses to changing anthropogenic aerosol forcing strength: regional and multi decadal ‐ variability

The causes of decadal variations in global warming are poorly understood, however it is widely understood that variations in ocean heat content (OHC) are linked with variations in surface warming. To investigate the forced response of OHC to anthropogenic aerosols (AA), we use an ensemble of historical simulations, which were carried out using a range of anthropogenic aerosol forcing magnitudes in a CMIP6-era global circulation model. We find that the centennial scale linear trends in historical OHC are significantly sensitive to AA forcing magnitude (−3.0 ± 0.1 × 10 5 (J m −3 century −1 )/(W m −2 ), R 2 = 0.99), but interannual to multi-decadal variability in global OHC appear largely independent of AA forcing magnitude. Comparison with observations find consistencies in different depth ranges and at different time scales with all but the strongest aerosol forcing magnitude, at least partly due to limited observational accuracy. We find that OHC is sensitive to aerosol forcing magnitude across much of the tropics and sub-tropics, and stronger negative forcing induces more ocean cooling. The polar regions and North Atlantic show the strongest heat content trends, and also show the strongest dependence on aerosol forcing magnitude. However, the OHC response to increasing aerosol forcing magnitude in the North Atlantic and Southern Ocean is either dominated by internal variability, or strongly state dependent, showing different behavior in different time periods. Our results suggest the response to aerosols in these regions is a complex combination of influences from ocean transport, atmospheric forcings, and sea ice.


Introduction
The ocean has taken up over 90% of the excess heat trapped in the earth system since 1955 due to human-induced global warming (Cheng et al., 2017;Levitus et al., 2012).Variations in ocean heat uptake (OHU) have also been robustly linked with decadal variations in the rate of surface warming (Marshall et al., 2015;Meehl et al., 2011).In particular, changes in OHU in the Pacific have been linked with the early 2000s slow-down in global surface warming (England et al., 2014;Kosaka & Xie, 2013;Oka & Watanabe, 2017;Stolpe et al., 2021), although Chen and Tung (2014) instead link the so-called "hiatus" with changes in heat transport in the Atlantic and Southern Oceans.Yin et al. (2018) also link the Pacific with the subsequent jump in global surface warming in 2014-2016.Understanding what contributes to variability in OHC is therefore important to understand if future climate variability is to be accurately predicted.Variability in anthropogenic aerosols (AAs) have been found to be an important source of changes in OHC (Delworth et al., 2005;Fiedler & Putrasahan, 2021), along with variability in green-house gases (GHGs), volcanic aerosols, and internal variability (IV).Anthropogenic aerosols are also Abstract The causes of decadal variations in global warming are poorly understood, however it is widely understood that variations in ocean heat content (OHC) are linked with variations in surface warming.To investigate the forced response of OHC to anthropogenic aerosols (AA), we use an ensemble of historical simulations, which were carried out using a range of anthropogenic aerosol forcing magnitudes in a CMIP6era global circulation model.We find that the centennial scale linear trends in historical OHC are significantly sensitive to AA forcing magnitude (−3.0 ± 0.1 × 10 5 (J m −3 century −1 )/(W m −2 ), R 2 = 0.99), but interannual to multi-decadal variability in global OHC appear largely independent of AA forcing magnitude.Comparison with observations find consistencies in different depth ranges and at different time scales with all but the strongest aerosol forcing magnitude, at least partly due to limited observational accuracy.We find that OHC is sensitive to aerosol forcing magnitude across much of the tropics and sub-tropics, and stronger negative forcing induces more ocean cooling.The polar regions and North Atlantic show the strongest heat content trends, and also show the strongest dependence on aerosol forcing magnitude.However, the OHC response to increasing aerosol forcing magnitude in the North Atlantic and Southern Ocean is either dominated by internal variability, or strongly state dependent, showing different behavior in different time periods.Our results suggest the response to aerosols in these regions is a complex combination of influences from ocean transport, atmospheric forcings, and sea ice.

Plain Language Summary
As well as emitting greenhouse gases that warm the planet, throughout the industrial era humans have also released substances known as aerosols into the atmosphere.In general, these aerosols reflect heat arriving at the surface of the planet and cause cooling, however we don't have a good idea how changes in the amounts of these aerosols changes ocean warming.Most of the heat that has built up in the climate system has gone into the ocean, so this is important to understand.We use a large computer model to look at how the ocean would have warmed with different levels of aerosols from humans in atmosphere, keeping other things like greenhouse gases the same.We find that the amount of aerosols in the atmosphere strongly affects the warming of the global ocean on timescales of many decades to hundreds of years.Whilst most of the ocean cools when more aerosols are added, parts of the Southern and North Atlantic oceans show different behavior.What is going on in these regions is probably the result of many things, such as changes in ocean currents, in sea ice, as well as in the amount of aerosols and where they are being emitted.BOLAND ET AL.
• Climate model analysis shows centennial trends of historical global ocean heat content (OHC) depend on anthropogenic aerosol forcing strength • Increased aerosol forcing leads to general cooling of the global ocean, but with significant regional and decadal variations • The strongest responses to aerosol forcing coincide with the strongest OHC trends, in the Southern and North Atlantic oceans

Supporting Information:
Supporting Information may be found in the online version of this article.
In general, studies of CMIP5-era GCMs which compare simulations with aerosol-only forcing and GHG-only forcing (such as Cai et al. (2006), Irving et al. (2019), andShi et al. (2018)) find that aerosols induce cooling focused in the Northern Hemisphere (due to the higher concentration of historical aerosols in this region), which then leads to more heat uptake in the Southern Hemisphere compared to the Northern, inducing an increase in northward heat transport, mostly via a strengthened AMOC.Menary et al. (2020) and Robson et al. (2022) find links between AMOC strengthening and increased AA forcing within CMIP6 models, which Robson et al. (2022) link to enhanced surface heat loss in the sub-polar North Atlantic.
Previous studies have largely focused on the effect the presence or absence of AAs have on global OHU, in comparison with GHGs.Given the uncertainty in the magnitude of the ERF of AAs, there is also a need to understand the impact of that uncertainty on OHU.We utilize a climate model ensemble specifically designed (as part of the UK SMURPHS project) to sample a wide range of historical aerosol forcings.Using the SMURPHS ensemble allows us to investigate the link between aerosol forcing and decadal variations in OHU, assessing the forced response of OHC to the magnitude of AA forcing in a CMIP6 generation model, using multiple ensemble members to improve statistical robustness (Goosse et al., 2005).Comparison with observations allows us to assess the performance of the ocean component of the model (see companion papers Dittus, Hawkins, et al. (2020), Dittus et al. (2021) for analysis of the atmospheric component), as well as ascertaining the best match for the magnitude of historical aerosol forcing.
This paper is laid out as follows: we first give background on the simulations and observations used in this study in Section 2, and then describe our findings in Section 3.This is broken up into analysis of changes in volume integrated OHC on a centennial scale (Section 3.1), spatial patterns of integrated OHC on multi-decadal scales (Section 3.2), and sensitivity of these patterns to aerosol forcing magnitude (Section 3.3).We finish by discussing the implications of our results in Section 4.

Simulations and Observations
All results are from an ensemble of historical scaled aerosol emissions simulations, conducted with the HadGEM3-GC3.1-LLglobal climate model, as describe in detail in Dittus, Hawkins, et al. (2020).The SMURPHS ensemble was designed to sample a plausible range of historical aerosol forcing (Booth et al., 2018), with 2014 (the end of the historical simulation) AA ERF ranging from −0.38 to −1.50 Wm −2 , which spans most of the 95% confidence interval presented in IPCC AR5 (Boucher et al., 2013).The effective radiative forcings were calculated using a series of fixed SST (Sea Surface Temperature) runs as described in Dittus, Hawkins, et al. (2020).The negative sign of the ERF indicates the addition of further AAs tends to cool the earth's surface.
The targeted aerosol forcings were achieved by applying a constant scaling factor in space and time to the standard historical CMIP6 AA and precursor emissions (Hoesly et al., 2018), see Dittus, Hawkins, et al. (2020) for further 10.1029/2022JC018725 3 of 20 details.Five scaling factors were used (Table 1), each with five ensemble members spanning the full historical period 1850-2014, giving a total of 25 ensemble members.The 1.0 forcing case is the standard CMIP6 historical emissions scenario.Hobbs et al. (2016) demonstrate that de-drifting OHC in CMIP5 models closes the time-varying energy budget on decadal scales.Assuming similar behavior in CMIP6 models, all calculations of OHC in this work were de-drifted using the linear trend in the same OHC calculation from the first 500 years of the HadGEM3-GC3.1-LLpre-industrial control run.The offset of the starting points of each ensemble member was also accounted for, see Figure S1 in Supporting Information S1 for an illustration.This was carried out for each integral of heat content separately, so at each depth range, in each basin, and at each latitude-longitude or latitude-depth point.Following the advice of Dittus, Hawkins, et al. (2020), we advise caution in interpreting data from before 1900 due to the impacts of any possible small shocks introduced by the abrupt change in AA forcing at 1850 (1900 has been marked with a dashed line in all time series plots).
Ocean basins were defined using the standard masks for this ocean model configuration.A proxy volcanic activity time series was generated by summing the time series of long-wave absorption due to volcanic aerosols over latitude, atmospheric model level, and wavelength.Qualitatively similar time series were produced summing any of the other volcanic forcing fields.
From Section 3.2.2on, we concentrate on two 30-year time periods : 1960-1991 and 1980-2011.We choose to focus on these in particular for two reasons: firstly, they overlap with the observational time series, and secondly, because they cover the period of increasing OHC trends (see Section 3.2.1).Both time periods have distinctive patterns of historic CMIP6 AA emissions, which are amplified or dampened by the AA scaling factor: During 1960-1991, historic emissions are dominated by positive trends in Europe, Asia, and North America.From 1980 to 2011, sharp falls in Europe and North America are contrasted by continued rises across Asia and positive contributions from South America (Dittus, Hawkins, et al., 2020 their Supporting Information S1).
The observations used in this study come from two datasets.First, from the Chinese Academy of Science Institute of Atmospheric Physics Ocean Gridded Product (Cheng et al., 2017), referred to as IAP17, we use the monthly 0-700 and 0-2,000 m time series, downloaded from http://159.226.119.60/cheng/ in January 2021.Secondly, from the World Ocean Atlas 2018 (Boyer et al., 2018), referred to as WOA18, we use the yearly 0-700 m and pentadal 0-2,000 m Global OHC time series and associated standard error, downloaded from https://www.ncei.noaa.gov/products/world-ocean-atlas in July 2021.

Simulated Ocean Heat Content
OHC was calculated from the model potential temperature fields (θ) as follows: where c p is the heat capacity of sea water (taken to be 3,850 J/(kg C)), ρ 0 is a reference density (taken to be 1.027 × 10 3 kg/m 3 ), and V is the volume considered, which can be the global ocean or a given basin, and/or a given depth range.Whilst the surface of the model ocean is allowed to vary in time, the impact this has on the OHCs presented here is at maximum a few percent, and so we use a fixed volume for simplicity.
Comparing global OHC anomalies (with respect to the 1850-1900 ensemble mean) from different ensemble members (Figure 1a), the effect of different AA emissions scale factors becomes apparent at around 1950, with a spread on the order of 5 × 10 5 J/m 3 between the 0.2 and 1.5 experiments by the end of the simulations in 2014 (Figure 1a).The 0.2 ensemble members warm rapidly, with a close to linear increase in time, as the weaker AA forcing offsets less GHG-induced warming.The 1.5 ensemble members cool until approximately 2000, as the stronger AA forcing more than offsets the GHG-induced warming.Atlantic OHC shows the most divergence between forcing ensembles by the end of 2014 (Figure 1b), whereas the Indian ocean shows considerable overlap between ensembles for much of the 20th century (Figure 1e).Additionally, the form of the time series varies by basin-regardless of AA forcing factor, Atlantic OHC exhibits clear multi-decadal variability around an overall warming trend.The Pacific and Indian OHC anomalies become negative in the latter half of the 20th century in the 1.0 and 1.5 simulations.Spatial variations in sensitivity to AA forcing magnitude are explored further in Section 3.3.1.
Considering contributions to global OHC by depth range, it is also apparent that the influence of AA forcing factor decreases with depth.By 2015, the largest spread in global OHC anomalies is at upper-depths (0-300 m, Figure 1k, on the order of 30 × 10 5 J/m 3 ), and the lowest is in the deep ocean (2 km+, Figure 1u, on the order of 1 × 10 5 J/m 3 ), where there is overlap in ensemble members for the entire simulation.

10.1029/2022JC018725
5 of 20 sensitivity, even if there is still considerable overlap between ensemble members.The depth dependence of the sensitivity of OHC to AA forcing magnitude is discussed more in Section 3.3.2.

Centennial Scale Changes: Linear Analysis
Centennial scale changes in OHC in all basins and all depth ranges (apart from global 2 km+ OHC) are significantly linearly correlated with the 2014 ERF.Table 2 gives the squared Pearson correlation coefficient (R 2 ) and slope of the linear fit for linear correlations between changes in volume-scaled OHC from 1850 to 1870 to 1995-2015 for all 25 ensemble members and the magnitude of the 2014 ERF strength given in Table 1.We use the magnitude of the 2014 ERFs rather than the signed values in order that a negative correlation indicates stronger AA forcing leads to a decrease in OHC, which is more intuitive.
As outlined in Appendix A, we expect differences in the rates of change of OHC to reflect the sensitivity of OHU to AA forcing magnitude and any feedback effects on other forcings (Equation A2).The R 2 of 0.99 for the linear dependence of centennial changes in OHC on AA forcing magnitude indicates that the impact of feedback effects and non-linear sensitivities are negligible on a globally integrated, centennial scale.
We expect differences in the rates of change of in OHC for individual basins and depth ranges to reflect the sensitivity of ocean circulation to forcing magnitude, as well as changes in OHU and feedbacks (Equation A4).Thus the slopes of the fit indicate the combined impact of OHU changes and circulation changes.The lower R 2 values for the individual basins and depth ranges indicate the relative importance of non-linear impacts and feedbacks on other forcings on OHU and transport for these sub-domains.
The strength of the relationship (indicated by the slope of the linear fit) is strongest for the Atlantic at all depths.
In the deep ocean (2 km+), the Atlantic and Southern OHC fits are similar in strength but opposite in sign-the 2 km+ Atlantic OHC is the only volume to show a positive relationship between OHC and AA 2014 ERF magnitude, indicating increased AA forcing increases deep Atlantic OHC.This is linked with changes in the AMOC, as discussed in Section 3.2.3.The Global deep ocean shows no significant linear relation with AA 2014 ERF on a centennial timescale, due to the Southern OHC, and, to a lesser extent, the Indian and Pacific OHC, which have negative linear relations with AA 2014 ERF magnitude, acting in combination to offset the Atlantic relation of the opposite sign.
While the linear fits show that the magnitude of trends in OHC on centennial scales are sensitive to aerosol forcing magnitude, there is considerable non-linear behavior in many basins and at many depth ranges at decadal scales (Figure 1).This is investigated further in the following section.

Centennial Scale Changes: Non-Linear Analysis
In order to assess the sensitivity of multi-decadal variability in OHC to AA forcing magnitude, we fit a polynomial to the full-depth global and basin-wise OHC anomaly time series (Figures 1a-1e).The degree of polynomial was Note.R 2 indicates the square of the Pearson correlation coefficient, with bold indicating statistical significance at the 99% level.Slope indicates the slope of the linear fit, in units of 10 5 J/m 3 /century/(W m −2 ).chosen by experimenting with different degrees and choosing the lowest order fit (to favor simplicity and avoid over-fit) that provided a relatively small residual (at least one order of magnitude smaller than the fit).A fourth order polynomial fits these criteria well for global OHC and basin-wise OHCs (Figures 2a-2e), with multi-decadal variability captured by the polynomial fits and residuals an order of magnitude smaller (Figures 2f-2k).
To determine the dependence of the multi-decadal variability on AA forcing magnitude, we perform a principal component analysis over the multi-ensemble dimension of the 25 polynomial time series shown in Figures 2a-2e.The form of the first and second principal components for the global and basin-wise polynomial fits are shown in Figure 2k-2o.This is similar to analysis of global surface temperature time series from works such as Schlesinger and Ramankutty (1994), who carry out singular spectrum analysis on de-trended global and regional time series.
Our analysis indicates differences in multi-decadal linear trends between ensemble members are driven almost entirely by differences in aerosol forcing magnitude: The first PCs (blue lines), which take the form of multi-decadal linear trends, explain 99% of the global variance (and at least 92% in the basins), and the weights of the first PCs are extremely significantly correlated with AA forcing magnitude for all basins and globally (R 2 > 0.95).Indeed the form of the first PC resembles the form of the ERF time series, which shows an increasing positive trend toward the end of the time period (Dittus, Hawkins, et al., 2020).The ERF time series includes all forcings, with the prominent increase in the ERF time series primarily due to GHGs.The overall shape of PC1 mainly reflects the impact of GHG forcing, while the fact that the weights are highly correlated with the AA ERF indicates that the time series is modulated by AA.We note that this analysis is indicative as a full attribution of all modes of variability has not been carried out.
Differences in multi-decadal non-linear variability, represented by the second PC (orange lines in Figure 2k-2o), are responsible for very little of the differences between the OHC time series, explaining maximum 4% of the variance, and the weights of the second PC are not significantly correlated with AA forcing magnitude (R 2 ≈ 0).This indicates that AA forcing magnitude does not drive differences in multi-decadal non-linear variability in large scale OHC, and therefore differences in multi-decadal non-linear variability are driven by other forcing factors (kept constant in these ensembles).2k) are small, indicating that sub-decadal OHC variability (as defined here by the residuals of the 4th order polynomial fit) is not primarily driven by differences in AA forcing magnitude in this model.Small/some differences on these timescales may exist (e.g., following volcanic eruptions) but are not investigated further here.
Sharp drops in residual OHC are associated with spikes in volcanic activity (Figure 2, thick green lines), consistent with Church et al. (2005) and Gleckler et al. (2006).Additionally, there is a circa 70 year periodic feature in the global OHC residual, also apparent in the Atlantic and Pacific OHCs to lesser degrees, consistent with the results of Schlesinger and Ramankutty (1994), who attribute the periodic feature to IV originating in the Atlantic.It could also be the result of a fitting artifact, or indicative of known residual modes of multi-decadal IV such as the Atlantic Multidecadal Oscillation (AMO, Deser et al., 2010) and/or the Interdecadal Pacific Oscillation (Parker et al., 2007), although it should be noted that Mann et al. (2021) argue that the AMO is entirely driven by volcanic forcing and not internally generated.The amplitude of the periodic feature is small compared with the multi-decadal variability in the polynomial fits, so we have not investigated further.
Using a LOWESS (LOcal WEighted Scatterplot Smoothing) fit instead as in Cheng et al. (2022, not shown) results in a smaller magnitude residual and a different form of PC2, but does not change the form of PC1 or subsequent interpretation.

Ocean Heat Content Trends
Time series of OHC trends were calculated from the un-scaled global OHC as defined in Equation 1 (without the factor 1/V to allow for easier comparison with observations) as follows: At each time series point, a centered 30 years linear regression was calculated.For the model OHC, we used the linregress function from the scipy python library (Virtanen et al., 2020).For the observed OHC, we used the WLS function from the statsmodels python library (Seabold & Perktold, 2010), with the weights =  1∕ 2  , where σ E is the standard error provided with the observations.Both functions provided a standard error in the linear slope, as well as the slope itself.

Ocean Heat Content Trends Versus Observations
In order to compare our modeled OHC with those from observations (WOA18 and IAP17, see Section 2 for details), we calculated OHC trends for 0-700 and 0-2,000 m.Absolute values of simulated OHC are less likely to match observations, and the uncertainty in observations increases at earlier times due to measurement sparsity, whereas trends have relatively lower uncertainty, even when taking the uncertainty at individual times into account.The standard error provided with the observational datasets do not take into account all sources of uncertainty (Carton & Santorelli, 2008;Wang et al., 2018), and so we show two different products to indicate the magnitude of additional uncertainty.
The time series of simulated OHC trends (colored lines, Figures 3a-3j) and observed trends (black solid and dashed lines, Figures 3a-3j) have standard errors one or two orders of magnitude (respectively) smaller than the OHC trends themselves, and so are not shown.
Both 0-700 and 0-2,000 m simulated OHC trends drop to a local minimum at around 1965 (representing the trend for 1950-1981), even becoming negative for the larger forcing factors, indicating ocean heat loss from these depth ranges to either the atmosphere or greater depths.This corresponds to the period with the greatest increase in aerosols (Dittus, Hawkins, et al., 2020).From 1965 on, OHC trends for both depth ranges increase for all forcing factors, peaking at the end of the simulation.This is consistent with the form of the GMST (global mean surface temperature) time series in Dittus, Hawkins, et al. (2020), which show faster than observed warming from 1990 onwards.Dittus, Hawkins, et al. (2020) suggest this could indicate a possible warm bias in the transient climate response (TCR) of the model.
The observational OHC trends vary less than the simulated OHC trends: both are relatively flat until around 1970-1990, when they begin to increase, with the IAP17 data set beginning to rise before the WOA18 data set in both depth ranges (Figures 3a-3j).The 0.4 and 0.7 scaling trends show the most similarity to one or other of the observations for many decades in both depth ranges.The 0.2 and 1.0 scaling trends show some similarity at the start or ends of the observational time series.The 1.5 scaling trends are the only to not match observations in any time period or depth range.This is summarized in Figures 3k and 3l, which show the linear trends for 1940-1970, 1960-1990, 1980-2010, and a single 70-year trend for 1955-2014.
Overall, the observations imply a scaling of 0.2-0.7 for the periods 1940-1980, when considering the smallest residuals between the trend time series for both depth ranges.Similarly, from 1980 onwards, the 0-700 m observations imply a scaling of 0.4-1.0, and the 0-2,000 m observations imply a scaling of 0.2-0.7.This is consistent with the results of Dittus, Hawkins, et al. (2020), who also find the 0.4 and 0.7 scaling simulations match observations of global mean surface temperature, even when accounting for the warm bias in the model's TCR.Anthropogenic aerosols are not evenly distributed around the planet (Stern, 2006), thus changes in the forcing factor will amplify/dampen regional differences and the resultant impacts on OHC.To further investigate the regional sensitivity of OHC to aerosol forcing magnitude, we look at spatial patterns in both latitude/longitude (Section 3.2.2) and latitude/depth (Section 3.2.3).
Figure 3. Observations of ocean heat content (OHC) trends are inconsistent with the 1.5 forcing factor for the period 1980-2010, and with both 1.0 and 1.5 for the periods 1940-1970 and 1960-1990.Panels (a-j) show OHC trends over time, by ensemble member (color) and depth range-0-700 m (a, c, e, g, i) or 0-2,000 m (b, d, f, h, j).The black solid lines are derived from OHC observations from IAP17, the black dotted lines from WOA18.Panels (k and l) summarize the upper panels by taking data points from three, 30-year periods (1940-1970, 1960-1990, 1980-2010)

Depth-Integrated Ocean Heat Content Trends
To investigate the spatial distribution of OHC trends in the SMURPHS ensemble,we calculate the depth-integrated OHC as follows: where ∫dz indicates an integral over the full depth of the model.
We calculate trends in OHC xy using linear regression as in Section 3.2.1 at every latitude and longitude point for the 30 year periods 1960-1990 (Figure 4a) and 1980-2010 (Figure 4b).Globally, the simulations show relatively low OHC xy trends in 1960-1990, and relatively high OHC xy trends in 1980-2010 (as in Figure 3).Figure 4 shows the global ensemble mean trends, regional plots of the trends in the North Atlantic and Southern Ocean for each ensemble member can be found in Figures S2 and S3 in Supporting Information S1.
The patches of warming in the North West Pacific and across the North Atlantic in 1980-2011 are consistent with observations of trends in 0-300 m OHC from ocean state estimates (Meyssignac et al., 2019).The strong OHC trends in the Southern Ocean are consistent with observations of strong abyssal warming at 2,000 m+ in the Southern Ocean (Desbruyères et al., 2016;Purkey & Johnson, 2010).
The 1.0 and 1.5 factor simulations show significant regions of negative OHC xy trends, covering large areas in 1960-1990, then mostly confined to the tropical Pacific in 1980-2010, indicating that heat is being lost and/or redistributed from these regions and time periods.
Overall, the spatial patterns of OHC xy trends can be seen to depend on time period, AA forcing factor, and on ensemble member (ensemble member standard deviation is shown by the gray contours in Figure 4).The sub-polar North Atlantic, in particular, shows extremely large variability between ensemble members, with variability peaking around the path of the Gulf Stream.During 1960-1991, different ensemble members with the same forcing factor show opposite-signed patterns of OHC xy trends (Figure S2 in Supporting Information S1).In 1980-2011 the pattern is slightly less variable, with broad warming and a cold hole, which appears most often south of Iceland, but there are still large variations in the pattern size and location.

Zonally Integrated Ocean Heat Content Trends
We calculate the zonally-integrated global OHC as follows: where ∫dx indicates a zonal integral over the global ocean.When calculating for an individual basin, we scale by the zonal integral for each grid cell in order to allow the comparison of trends to be independent of basin size: where r E is the radius of the earth.
Looking at OHC yz trends (Figure 5 The near-surface warming is stronger in the Northern Hemisphere for all forcings and both time periods.This asymmetric warming is also seen in observations of the Atlantic by Zanna et al. (2019), who attribute it to changes in ocean circulation.As with patterns of OHC xy trends, in general the ensemble member variability in OHC yz trends is largest where trends are largest.
BOLAND ET AL.

Depth-Integrated OHC Trend Sensitivity to Aerosol Forcing Factor
In order to more robustly assess the impact of changing AA forcing factor on the OHC trends calculated in Section 3.2.2,we regress the OHC xy trends for all 25 ensemble members against the AA ERF magnitude for each experiment (table 1) at each latitude-longitude point.Figure 6 shows (in color) the slopes of this regression for 1960-1990 and 1980-2010, with stippling showing where the regression is significant according to the student T-test.
In other words, the impact of increased aerosol forcing magnitude on OHC xy trends across all ensemble members is shown, independent of other forcings (which are identical across all experiments).As it is a linear regression, this will not reflect non-linear effects.Such effects might include processes or regions that are in some way saturated with respect to OHC, such that increasing the radiative forcing beyond a certain magnitude does not result in an increase in OHC trends, but decreasing radiative forcing results in a decrease in OHC trends, or the converse.Additionally, regions where the variability in OHC trends is large compared to the dependence on AA forcing magnitude within the bounds tested in these simulations will fail the significance test.There is broad negative linear sensitivity of OHC xy trends to aerosol forcing magnitude across much of the tropics and sub-tropics, indicating increasing the forcing magnitude results in a reduction in OHC xy trends in these regions, with a magnitude of 2-5 × 10 −6 [J yr −1 ]/W in both time periods, with peaks in dependence at high latitudes.There are a few patches of positive linear sensitivity, notably in the sub-polar  or polar  North Atlantic, and close to Antarctica.
The sub-polar North Atlantic and Southern Ocean show regression patterns that differ between the two time periods, and both show large regions without statistical significance, likely because these regions have strong IV (Figure 4).The strongest dependence on forcing magnitude is found in a dipole pattern centered over Iceland, which does show statistical significance and is of opposite signs in the two different time periods.The depth-dependence of the regression of OHC trends on aerosol forcing magnitude is discussed in Section 3.3.2,and possible links with the overturning circulation are discussed in Section 4.
The North Pacific is another region where the regression differs between 1960-1991 and 1980-2011.Similarly to the North Atlantic and Southern Oceans, the North Pacific is an OHC xy trend hotspot and a region of strong IV in 1980-2011 (Figure 4b), and not linearly dependent on AA forcing magnitude by our test (Figure 6b).Overall, the patterns of OHC xy trend regression on forcing magnitude in the North Pacific are very similar to the patterns of surface air temperatures (SAT) regressions from 1951 to 1980 and 1981-2012 (see Dittus et al. (2021) their Figure 8).
The SAT regression in 1981-2012 resembles a negative Pacific Decadal Oscillation, and this link is investigated in Dittus et al. (2021) in both the SMURPHS ensemble and other CMIP6 GCMs.They conclude that AA can induce an increase in North Pacific sea-level pressure (SLP) which promotes a negative PDO (Pacific decadal oscillation) in this time period.This introduces a relative cooling in the North Pacific surface air temperature in the same regions where we find significant OHC xy trend regression on forcing magnitude.They also note a high level of IV in the North Pacific SLPs in the SMURPHS ensemble and other GCMs, consistent with the high OHC trend variability in this region that we find in the SMURPHS ensemble.
Overall Pacific SATs in the SMURPHS ensemble warm in 1981-2012 in all but a few patches, due to stronger GHG forcing than AA forcing in this time period (see Dittus et al. (2021); Figure 4) which is reflected in the warming of the upper 200 m of the ocean in the same time period (Figure 1).However, both the 1.0 and 1.5 scaling ensembles shows overall negative OHC xy trends for large parts of the Pacific (Figure 4b), which is due to negative OHC trends in mid-depths (Figures 1m and 1r; Figure S3 in Supporting Information S1).The depth-dependence of the OHC trends regression on aerosol forcing magnitude is addressed in the next section.

Zonally-Integrated OHC Trend Sensitivity to Aerosol Forcing Factor
Regressing the trends in  O ĤC (defined in Equation 4) against the 2014 AA ERF, as in Section 3.3.1,indicates how the impacts of aerosol forcing strength vary with depth, Figures 7a-7d.The regressions are generally negative outside of the polar regions and above ∼1,500 m for both 1960-1991 and 1980-2011.However, the strength of the regression varies with depth, time period, and basin.As before, where there is no stippling (indicating statistical significance of the linear regression), then either IV is dominant, there is non-linear dependence on aerosol forcing magnitude, or another forcing term is dominant.
The Atlantic shows the strongest overall relative dependence on aerosol forcing magnitude (Figures 7a and 7b), with strong, significant negative dependence in the upper 1,000 m in the tropics in both periods.The dipole of positive and negative sensitivities centered on Iceland seen in the OHC xy regression (Figures 6c and 6d) can be seen to extend to depth in the OHC yz trend regressions (Figures 7a and 7b).There is no significance in the OHC yz regression in the positive patch at ∼50N in 1960-1991, whereas there is in the OHC xy regression, likely because the patch of significance is limited to the east of the sub-polar North Atlantic.
The high-latitude Atlantic shows a dipole structure at depth, where negative sensitivity switches to positive sensitivity below ∼1,500 m, indicating increased aerosol forcing magnitude is resulting in deep warming.This is also present in the Southern Ocean, but the scaled sensitivity is far weaker than the Atlantic (Figures 7e and 7f).
In both the Pacific and Indian basins,  O ĤC shows relatively strong, significant negative dependence on aerosol forcing magnitude in the upper ∼600 m.There are also subsurface patches where the regression is positive but not significant-the latitudes of these patches coincide with the latitudes of non-significant patches in OHC xy trends.
In the Southern Ocean, increased aerosol forcing leads to decreased  O ĤC trends in the upper 1,000 m at 40-50°S, most pronounced in 1960-1991 (Figures 7e and 7f).In 1980-2011, increased aerosol forcing also leads to relative cooling at the surface close to the continent and warming at depth (Figure 7f).

Multi-Decadal Variability
We find that aerosol forcing magnitude is responsible for changes in multi-decadal global OHC linear trends at global and basin-wide scales (R 2 ≥ 0.92), but that interannual to multi-decadal variability is relatively insensitive to forcing magnitude.
The reason we do not see changes in multi-decadal non-linear variability in our ensemble may be because the overall ERF time series (including all forcings) does not show a lot of multi-decadal variability, instead showing large interannual variability and long-term linear trends.It also may be that the impacts of multi-decadal variability are confined to near the surface, and don't impact on depth-integrated OHC.This would be consistent with the results of Qin et al. (2020) who find multi-decadal variations of volcanic aerosols and AA are responsible for multi-decadal variations in SSTs in all three major ocean basins, using models and observations.

Comparison With Observations
Trends in 0-700 m OHC are most consistent with observations for the 0.4-1.0scaling experiments, consistent with Dittus, Hawkins, et al. (2020), who find the 0.4 and 0.7 scaling experiments most consistent with GMST  14 of 20 observations.Trends in 0-2,000 m OHC are most consistent with observations for the 0.2-0.7 experiments.This inconsistency between the different depth ranges is likely due to a combination of factors.First, discrepancies in the model's representation of reality, both in terms of impacts of AAs on OHU and circulation changes.Secondly, as discussed in Section 3.2.1,uncertainties in observations of OHC are not fully quantified and so it may be that the true uncertainty in observations spans the same forcing factors.

Regional Changes
Regional changes in OHC due to changes in aerosol forcing are due to a combination of changes in OHU and ocean transport.Increased aerosols in the atmosphere decreases shortwave radiation which leads to a relative cooling at the surface of the ocean which can then be transported into the interior.Indeed we find over large parts of the ocean, particularly the tropics and sub-tropics, that the impact of increased aerosol forcing magnitude is a fairly uniform linear cooling of depth-integrated OHC (Figure 6).
Changes in ocean heat transport can be caused by direct changes in radiative forcing and subsequent changes in winds and fresh water forcings.Even if AAs were spatially uniformly distributed, the impact of changes would not be spatially uniform due to the impacts of ocean circulation.The impact of ocean circulation changes can be seen most clearly in the polar regions, where the vertical limbs of the global meridional overturning circulation are located.Zanna et al. (2019) demonstrate that up to 50% of the increase in ocean heat stored in the mid-latitude Atlantic is due to transport changes.The impacts of increased aerosol forcing magnitude on depth-integrated OHC in these regions contains regions of warming, and are stronger than in other regions (especially in the Atlantic), and highly variable (Figure 6), all likely due the additional impacts of heat convergence/divergence.The increased vertical transport in these regions leads to a stronger impact of aerosol forcing changes at depths (Figure 7).Thus, while we expect the impact of aerosols on OHC to be non-uniform spatially, we might expect that the impact is similar at different times.In fact, we find that the regional impacts can vary significantly with time period.We focus in this paper on two 30-year periods near the end of the historical simulation , both because they are the time periods with most observations and because the model behavior is different in both periods-1980-2011 sees a strong acceleration of global warming with subsequent impacts on cryosphere (Andrews et al., 2020;Dittus, Hawkins, et al., 2020), shows a reversal in the trend in Atlantic Multidecadal Variability (AMV) (see Andrews et al. (2020) and Figure S8 in Supporting Information S1), an increase in the trend of equivalent radiative forcing from all sources (Dittus, Hawkins, et al., 2020), although there are significant regional variations (Dittus et al., 2021) and the overall AA forcing is stable in this period.Indeed, Andrews et al. (2020) link the change in AMV (and AMOC) trend sign in the standard historical forcing simulation (our 1.0 forcing case) with regional variations in aerosol forcing.We now discuss the differences between the two forcing periods and how they compare with literature for each basin in turn, with the Indian and Pacific basins discussed together.

Atlantic
Atlantic OHC shows the strongest relative response to increased AA ERF in both time periods by all measures presented-on a centennial, basin-integrated scale (Table 2), and on a multi-decadal scale in both depth-integrated (Figure 6) and zonally-integrated (Figure 7) trends.
We find that the strength of the AMOC is dependent on aerosol ERF strength (Figure S8 in Supporting Information S1).Additionally, there is a significant relationship between the centennial trend in AMOC strength and AA ERF (not shown), consistent with the results of Cai et al. (2006), Collier et al. (2013), Irving et al. (2019), Menary et al. (2020), Robson et al. (2022), andShi et al. (2018).However, there is not a clear link on shorter timescalesthe links between the 30-year trends in AMOC and AA ERF strength are not significant for 1960-1991 and of opposite sign to 1980-2011 (see Figures S8f andS8g in Supporting Information S1).The statistical correlations are also weak and of similar magnitude to IV, implying either multi-decadal variability in AMOC strength is not strongly controlled by AAs or that the processes are more complex than a correlation can represent.We hypothesize that AA regional trends rather than absolute AA forcing strength could influence AMOC strength: we see a centennial scale strengthening in the AMOC alongside increasing AA forcing up until circa.1980, then regional decreases in AA across Europe and N America from 1980 onwards (strongest in the 1.5x scenario), alongside a sharp rise in GHGs, drive a decreasing AMOC (strongest in the 1.5x scenario).BOLAND ET AL.The strengthening of the AMOC with increased AA forcing is consistent with the significant positive sensitivity of depth-integrated OHC in the Sub-Polar North Atlantic to AA ERF magnitude in 1960-1991.A similar pattern is seen in the depth-integrated temperature-trend response in Cai et al. (2006) (their Figure 4), and the SST responses in Collier et al. (2013), Robson et al. (2022), andShi et al. (2018).This pattern is consistent with the increased convergence of heat in the region, which Robson et al. (2022) and Shi et al. (2018) find leads to increased heat loss to the atmosphere, decreasing upper ocean stratification and further strengthening the AMOC.This feedback effect may explain the link between centennial AMOC trends and AA ERF magnitude.
During the period 1980-2011, we still see relatively strong cooling in the south Atlantic in response to increased AA forcing magnitude, but we no longer see warming south of Greenland (Figure 6b).Instead, the regression in this region resembles the pattern of OHC trends in the simulations (Figure 4), with a "warming hole" south of Greenland, linked in Liu et al. (2020) to a slowing of the AMOC, consistent with the weak dependence of AMOC trend on AA forcing magnitude in this time period (Figure S7g in Supporting Information S1).We still see an increase in warming at depth in the north Atlantic (Figure 7b) but the upper OHC response to increased aerosols appears to be under the influence of multiple interacting and possible competing processes-a strong but decreasing AMOC, a strong slowing in the loss of Arctic sea ice (Dittus, Hawkins, et al., 2020), and changes in Northern Hemisphere aerosol composition with North American and European emissions dropping against a background of increasing Asian emissions (Dittus et al., 2021).

Indo-Pacific
In the Pacific, Cai et al. (2006) find the inclusion of aerosols in historic gcm simulations induce a cross-equatorial overturning circulation, with northward transport at the surface and southward transport down to circa 800 m depth, inducing warming north of the equator and cooling south of it.This resembles the pattern of OHC sensitivity to aerosol forcing magnitude in 1960-1981, both depth-integrated (Figure 6a) and zonally integrated (Figure 7c), although the latter is not statistically significant.In 1980-2011 the Pacific regression pattern is instead dominated by a PDO-like signal, linked by Dittus et al. (2021) to a surface pressure response to increased aerosol forcing, possibly triggered by a Rossby wave response to increased Asian aerosol emissions since the 1980s.Sun et al. (2022) find that a weakened AMOC leads to compensating northward transport in the Indo-Pacific basins, causing heat to be redistributed from the Atlantic to the Indo-Pacific basins via the Southern Ocean, leading to sub-surface (∼1-3 km) warming.The implied opposite effect due to the AMOC strengthening is consistent with the sub-surface negative sensitivity in both the Pacific and Indian basins here.
Whilst the Pacific doesn't stand out when comparing the relative sensitivity of OHC to AAs between basins, it is the largest basin and contributes the most in absolute terms to the sensitivity of global OHC to AAs [not shown].Thus, the Pacific's broad statistically significant negative sensitivity to AAs (see Figure 6) indicates it plays an important role in the ocean's energy budget during the historical period.

Southern Ocean
Whilst the Southern Ocean is far from the strongest aerosol forcing locations, it is clear that its dominant role in ocean heat storage is significantly dependent on aerosol forcing magnitude.Especially notable is during 1960-1991, when GHG forcing is relatively weaker than in 1980-2011, the strongest aerosol forcings (1.0 and 1.5) lead to a cooling in much of the upper 2,000 m in the Southern Ocean (see Figure 5a).
The warming at depth to the north of the Southern Ocean in both time periods in response to increased AA ERF magnitude is likely linked to the strengthening of the global meridional overturning circulation, as indicated by the changes in the AMOC discussed above.
In the Southern Ocean in 1960-1991 (Figure 7e), we see a cooling pattern that is the opposite of the GHG-induced trends in Southern ocean heat storage (Dias et al., 2020;Liu et al., 2018)-surface cooling north of 55S and a concentration of statistically significant cooling at 40-45S in the upper 1,000 m due to weakening of the overturning circulation and a reduction in heat convergence.This is consistent with a negative SAM (Southern Annular Mode)-like pattern in the regression of SLP against aerosol forcing found in this time period in Dittus et al. (2021) (although not statistically significant), implying decreasing zonal winds.However, Steptoe et al. (2016) find a similar negative response in the SAM to AA changes in CMIP5 models is not robust and model-dependent.
In 1980-2011, we instead see a concentration of surface cooling at 60S, and less relative cooling in the interior.This could be due to the contribution of warming at these latitudes in the Indian sector (Figure 6b).This may also be influenced by the positive SAM-like response to aerosol forcing in this time period (Dittus et al., 2021), which acts to increase the wind-induced overturning circulation, strengthening the GHG induced effects and relatively warming the interior, although there is still a cooling effect of aerosol forcing overall.Additionally, there is an increased loss of Antarctic sea ice extent in the model in 1981-2011 in both summer and winter (see Figure S9 in Supporting Information S1), which is slowed by increased aerosol forcing when considering the ensemble means, although there is large IV.Increased sea ice cover is consistent with relative cooling at the surface close to the continent, whilst the warming at depth below is consistent with the production of cold, dense waters slowing (Figure 7f).

Summary
Using a unique ensemble of 25 simulations of the historical climate with five different AA forcing time series, we have been able to determine the influence of aerosol forcing magnitude on OHC in a CMIP6-era gcm.We find that the 20th century volume-scaled global OHU sensitivity to AA forcing magnitude is −3.0 ± 0.1 × 10 5 (J m −3 century −1 )/(W m −2 ) for the HadGEM3-GC3.1-LLmodel.Centennial changes in the OHC of the major ocean basins and globally integrated depth ranges above 2,000 m also show significant linear dependence on AA forcing magnitude (R 2 ≥ 0.94), indicating that the impact of non-linear effects and feedbacks from other forcings are negligible.We find that aerosol forcing magnitude is responsible for changes in multi-decadal global OHC linear trends at global and basin-wide scales (R 2 ≥ 0.92), but that interannual to multi-decadal variability is relatively insensitive to forcing magnitude.
Trends in 0-700 m OHC are most consistent with observations for the 0.4-1.0scaling experiments, consistent with Dittus, Hawkins, et al. (2020), who find the 0.4 and 0.7 scaling experiments most consistent with GMST observations.Trends in 0-2,000 m OHC are most consistent with observations for the 0.2-0.7 experiments.Both results are consistent with Robson et al. (2022), who find that CMIP6 models with the strongest aerosol forcings are inconsistent with observations.We find the responses to increased AA strength is significantly dependent on region and time period.In general, the strongest responses are found in the regions where there are the strongest trends in OHC, specifically the North Atlantic and Southern Oceans.The responses in these regions are summarized in Figure 8, with proposed mechanisms.
The difference in aerosol sensitivity in the different time periods implies a strong state dependence of the aerosol impacts on OHC, such that the impact of aerosols forcing changes on OHC is different depending on a combination of some or all of: the ocean, atmosphere, and cryosphere state; the magnitude and distributions of other forcings (GHGs, volcanic, natural aerosols); the magnitude of the aerosol forcing itself (the ERF is generally higher and increasing for all forcing factor experiments in 1980-2011compared with 1960-1991, see Dittus, Hawkins, et al. (2020); Figure 1b).
Our results give, for the first time, a well-constrained estimate of the dependence of historic global OHU on aerosol forcing magnitude.Our results suggest that OHC could potentially be used to constrain the estimate of the true aerosol forcing magnitude, but that accurate and sustained measurements would be required.
We find that there is significant regional and decadal variability in the sensitivity of OHC to aerosol forcing magnitude in regions of high OHU.The uncertainty in aerosol forcing magnitude is not the dominant source of regional variability-instead the strong state dependence means that careful process based studies are required to disentangle the various mechanisms at play that determine the overall regional impact of aerosol forcing in different time periods.

Figure 1 .
Figure 1.The sensitivity of ocean heat content (OHC) to anthropogenic aerosol (AA) forcing factor varies by basin (column) and depth range (row).Each line represents the volume-scaled OHC of a single ensemble member, relative to the 1850-1900 ensemble mean, with the AA forcing factor indicated by the color (see legend).When scaled by volume, the four basin columns sum to the global column, and the four depth range rows sum to the full depth row.

Figure 2 .
Figure 2. Multi-decadal linear trends are sensitive to anthropogenic aerosol (AA) factor, but annual to decadal variability is largely independent of AA forcing magnitude.Panels (a-e) (the upper row) show 4th order polynomial fits to ocean heat content (OHC) anomalies w.r.t.1850-1900 ensemble means, globally and by basin (Figures 1a-1e), for all ensemble members.Colors indicates the AA forcing factor (see legend).Panels (f-j) (the middle row) show the residual OHC, calculated by subtracting the 4th order fits (upper row) from the same OHC anomalies, and smoothed with an 18 months low-pass Butterworth filter.The green line is a proxy for volcanic forcing, which is the same for all ensemble members, see text for details.Panels (k-o) (lower row) show the form of the first and second principal components of the 4th order fits in the upper row.The explained variance (Evr) and correlations of the PC weights with AA effective radiative forcing (R 2 ) are shown in each panel.
Figure3.Observations of ocean heat content (OHC) trends are inconsistent with the 1.5 forcing factor for the period 1980-2010, and with both 1.0 and 1.5 for the periods1940 -1970  and 1960 -1990.Panels (a-j) .Panels (a-j)  show OHC trends over time, by ensemble member (color) and depth range-0-700 m (a, c, e, g, i) or 0-2,000 m (b, d, f, h, j).The black solid lines are derived from OHC observations from IAP17, the black dotted lines from WOA18.Panels (k and l) summarize the upper panels by taking data points from three, 30-year periods(1940-1970, 1960-1990, 1980-2010), as well as for one 70-year period covering 1955-2014.Colored dots indicate ensemble members, colored crosses ensemble means, and black lines observations as before.
), we see upper 200 m warming dominates, especially in 1980-2011.Warming extends deeper in the Southern Ocean and in the Northern mid-latitudes for 1980-2011.Surface cooling is present in the 1.5 forcing scenario in places during 1960-1991, but is confined to 200 m+ in all other scenarios and in 1980-2011.Figure 5 shows ensemble mean trends, Figures S4-S7 in Supporting Information S1 show the trends by ensemble member for each basin.

Figure 4 .
Figure 4. Spatial patterns of depth-integrated ocean heat content (OHC) trends vary by ensemble member, by anthropogenic aerosol forcing factor, and time period.Colors indicate the ensemble mean depth-integrated OHC for each forcing factor for (a) 1960-1991 and (b) 1980-2011.Gray contours indicate the ensemble standard deviation, at 2, 4, and 6 × 10 7 J/m 2 /yr.

Figure 5 .
Figure 5. Patterns of zonally-integrated ocean heat content vary by ensemble member, by anthropogenic aerosol forcing factor, and time period: Colors indicate the ensemble mean OHC yz trends for each forcing factor for (a) 1960-1991 and (b) 1980-2011.Gray contours indicate the ensemble standard deviations at 2, 3.5, and 5 × 10 11 J/m 2 /yr (a) and 2, 4, and 6 × 10 11 J/m 2 /yr (b).Note that depth intervals are not constant, and color axis limits in (b) are twice those in (a).

Figure 6 .
Figure 6.Ocean heat content (OHC) trends are linearly sensitive to anthropogenic aerosol (AA) forcing magnitude in the tropics and east-basin sub-tropics: Colors indicate the regression of (a) 1960-1991 or (b) 1980-2011 OHC trends on 2014 AA effective radiative forcing magnitude for all 25 ensemble members.Stippling indicates where the trend is statistically significant using the student T-test.

Figure 7 .
Figure 7.The sensitivity of ocean heat content trends to aerosol forcing magnitude varies with basin, depth, and time period: Colors indicate the regression of (a) 1960-1991 or (b) 1980-2011  O ĤC trends on 2014 anthropogenic aerosol effective radiative forcing magnitude for all 25 ensemble members.Stippling indicates where the trend is statistically significant using the student T-test.Note that depth intervals are not constant but change at 1,000 m as indicated by the black line.
When considering absolute values of OHC, the Pacific has the largest volume and so contributes most to global OHC (not shown).Scaling the OHC by the basin volume allows for comparison of the relative changes in OHC Considering relative contributions to global OHC from each of the four main ocean basins, it is apparent that the sensitivity of OHC to AA forcing factor is not spatially uniform.

Table 2
Linear Correlations Between Changes in Volume-Scaled Ocean Heat Content 20-Year Means, From 1850-1870 to 1995-2015, and Present Day Aerosol Forcing Magnitude, Split by Depth Range and Basin 10.1029/2022JC018725 6 of 20