Interhemispheric Contrasts of Ocean Heat Content Change Reveals Distinct Fingerprints of Anthropogenic Climate Forcings

During recent decades, both greenhouse gases (GHGs) and anthropogenic aerosols (AAs) drove major changes in the Earth's energy imbalance. However, their respective fingerprints in changes to ocean heat content (OHC) have been difficult to isolate and detect when global or hemispheric averages are used. Based on a pattern recognition analysis, we show that AAs drive an interhemispheric asymmetry within the 20°‐35° latitude band in historical OHC change due to the southward shift of the atmospheric and ocean circulation system. This forced pattern is distinct from the GHG‐induced pattern, which dominates the asymmetry in higher latitudes. Moreover, it is found that this significant aerosol‐forced OHC trend pattern can only be captured in analyzed periods of 20 years or longer and including 1975–1990. Using these distinct spatiotemporal characteristics, we show that the fingerprint of aerosol climate forcing in ocean observations can be distinguished from both the stronger GHG‐induced signals and internal variability.

• A pattern recognition method is applied to extract the aerosol-induced pattern of interhemispheric ocean heat content (OHC) change • Anthropogenic aerosols dominate the historical interhemispheric OHC asymmetry within the 20°-35° latitude band over the 1975-1990 period • Observed OHC changes have an asymmetry that is consistent with the aerosol-induced OHC fingerprint

Supporting Information:
Supporting Information may be found in the online version of this article. , which demands a cross-equatorial atmospheric heat flux to the "colder" hemisphere. In model simulations, an interhemispheric asymmetry of the radiative flux at the top of the atmosphere and ocean heat uptake are predominately driven by AAs (Irving et al., 2019;Kang et al., 2021;Lembo et al., 2019;Luongo et al., 2022;Shi et al., 2022). The global observations from the Argo Programme show that the upper 2,000 m of the ocean has continued to warm since 2005, but that the SH ocean has stored more excess heat than the NH ocean (Roemmich et al., 2015;Wijffels et al., 2016). However, Rathore et al. (2020) find, using 11-year OHC trends (2005, that the observed hemispheric asymmetric ocean warming cannot be distinguished from internal climate variability. It thus remains unclear whether, where, and when aerosol-driven signals appear in the ocean interior. Much past work on this topic has focused on hemispheric or global averages, which reduces the system response to a set of 1-dimensional time series. We confirm below this approach is not effective in signal separation. In contrast, we contend that the meridional structure across different hemispheres provides additional information. Besides the spatially inhomogeneous nature, AAs also undergo a distinct non-monotonic time evolution Diao et al., 2021;Shi et al., 2022;Wang & Wen, 2022) due to the history of industrial development and emissions regulation. We show that this temporal information can also be exploited to track the AAs' impacts in the ocean interior.

CMIP6 Simulations
This study uses simulations from the Coupled Model Intercomparison Project Phase 6 (CMIP6) (Eyring et al., 2016). We choose 11 CMIP6 models that performed all of the following 4 experiments: historical all-forcing simulations (HIST), AA single-forcing simulations (AER), greenhouse gas single-forcing simulations (GHG), and natural-only simulations (NAT) (Gillett et al., 2016). The sizes of ensemble for each model/experiment are listed in Table S1 in Supporting Information S1.
The zonally and vertically integrated OHC is obtained as: where is the annually averaged potential temperature at time t, depth z, latitude y, and longitude x. The reference density is 0 = 1,025 kg/m 3 and the specific heat of water is = 4,000 J/(kg K). Except for the full-depth integration for the heat budget analysis (see below), the ocean depth h is set to 2,000 m.
Coupled climate models are prone to a long-term drift, which can contaminate the detection of long-term trends driven by external forcing. We estimate a model drift of OHC by fitting a cubic polynomial to the corresponding preindustrial control run over its full length in time for each model. Using the branch time information (the time when forced experiments branch out) provided in the file metadata, the correct segment of the cubic polynomial from piControl is subtracted from the forced simulation.

Ocean Heat Budget Analysis
The zonally integrated heat budget for the full-depth ocean is where ∂OHC/∂t is the trend of full-depth OHC during a certain period. Qnet is the zonally integrated net surface heat flux into the ocean, and OHT is the zonally integrated meridional ocean heat transport. The prime in these terms represents the anomaly after drift removal. As only a few CMIP6 models used in this work provide the OHT, a heat budget is formed for 6 CMIP6 models: CanESM5, CNRM-CM6-1, HadGEM3-GC31-LL, IPSL-CM6A-LR, MRI-ESM2-0, NorESM2-LM.

Signal-To-Noise-Maximizing Pattern Filtering Method
In this study, we use signal-to-noise-maximizing pattern filtering to extract the anomaly patterns, for which different ensemble members agree on the temporal evolution (R. C. Wills et al., 2018; R. C. J. Wills et al., 2020). In this method, the extracted patterns are associated with the maximization of the ratio of signal (e.g., variance of multi-model mean [MMM]) to the total variance (from all individual 54 members). One advantage of this method is that it makes the most of the spatiotemporal information and the exacted pattern reflects the best agreement among all the ensemble members. Full details of how this pattern-based method is used in this work are given in Supporting Information S1.

Estimated OHC From Observations
To obtain multi-decadal coverage of ocean observations we use EN4.2.2 and the Institute of Atmospheric Physics (IAP) gridded datasets to compare with the results from model simulations since 1950. The EN4.2.2 potential temperature products are available on a 1° × 1° grid with 42 vertical levels from the surface to about 5,500 m during 1900-2021. Here we analyze two different EN4.2.2 products, EN422_g10 (Gouretski & Cheng, 2020;Gouretski & Reseghetti, 2010) and EN422_c14 (Cheng et al., 2014;Gouretski & Cheng, 2020), which use different treatments of bathythermograph bias corrections. We also use the IAP temperature product, which has 41 vertical levels to the 2,000 m depth from 1940 to 2021 on a 1° × 1° grid. The use of different products attempts to capture the additional uncertainty due to space/time gap-filling assumptions. We find that differences in masked regions around the marginal seas have little effect on the results.

Zonal Mean Ocean Temperature Responses
The global zonal mean temperature trends from CMIP6 historical experiments are shown in Figure 1. Here, we just show latitudes within 35°S-35°N to focus on the asymmetric and symmetric patterns in the tropical/ subtropical regions driven by AER and GHG runs, respectively. The AA-induced asymmetry of cooling between two hemispheres is distinct with a stronger surface cooling in the NH and a subsurface maximum of cooling in the SH (Figure 1c). In addition, the cooling around 30°N is linked with subduction of cooling signals along the ventilated thermocline (e.g., Shi et al., 2023). The weaker AA cooling between 0 and 20°N in the lower thermocline is suggestive of the effect of a southward shift of the bowl of the subtropical gyre, as is the deep cooling around 20°S. In contrast, the temperature responses to GHG forcing are relatively hemispherically symmetric, and surface intensified everywhere ( Figure 1d). The weaker warming along the inside edges of the GHG response is suggesting of a poleward migration of the gyres in both hemispheres. In HIST runs, which experience both anthropogenic and natural forcings, GHGs dominate in driving surface warming (Figures 1a and 1b), but with a subsurface cooling in the vicinity of the equatorward edges of the subtropical gyre, with the SH cooling being much stronger than that in the NH. Interestingly, the cooling/warming signals in the deep tropics in the GHG and AA runs appear to largely cancel out in the HIST run. We revisit the physical processes associated with these changes in Section 4.2, and first focus on identifying a detectable fingerprint for application to the observation record.

Interhemispheric Latitudinal Difference Patterns of Forced OHC Responses
Interhemispheric asymmetry in energetics including OHC is receiving much recent attention because it is an important diagnostic for cross-equatorial energy transports by ocean and atmospheric motions and meridional displacements of the ITCZ (Kang et al., 2021;Luongo et al., 2022). Inspired by the contrast in meridional temperature responses (GHG are meridionally symmetric and AA are asymmetric), we form a new index that remains latitude as a dimension (vs. a simple difference of hemispheric averages). We use the de-drifted OHC to calculate the interhemispheric pattern of zonally integrated OHC change, which is defined as This new index is a function of latitude (y) from 0° to 90°, and denotes how much OHC is changed in NH latitudes relative to those at the same SH latitudes. In this work, we focus on the trend segments in 1950-2014. Following previous studies that regarded the equator as the "cutting line" and diagnosed the interhemispheric contrast (Irving et al., 2019), we use 0° latitude as the center to calculate OHC| ( ) . The annual mean ITCZ at 5°N could also be chosen for inter-hemispheric differencing but we find that the choice of the boundary latitudes (close to the equator) does not affect the results shown in this study. To most efficiently extract the forced variability from internal variability we apply the signal-to-noise-maximizing analysis to our new index, OHC| ( ) .
As we aim to extract externally forced signals by separating them from internal climate variability in this study, we need to define a threshold to evaluate when the forced response is significant. We deem the forced trend to be significant if the trends from most of 54 simulations have the same sign as the MMM trend. We use a 10% (e, f) Zonal-mean zonal wind stress from AER and GHG runs, respectively. Black curve is the trend between 1950 and 2014. Solid blue curve is the mean zonal wind. The dashed curve combines the change and mean of zonal wind stress, mean + × trend . As the trend is much smaller than the mean value, we use a = 100 to clearly show shifts in wind peaks.
threshold, so that the forced trend is deemed significant if more than 49 out of 54 trends exhibit the same sign as the MMM trend.

Dominance of AA and GHG Effects in Different Latitudes
The leading signal-to-noise-maximizing patterns of interhemispheric zonally integrated OHC (OHC| ( ) in Equation 3) anomalies relative to 1950-2014 mean and the corresponding time series from HIST, AER, and GHG ensembles are shown in Figure 2. With a signal-to-noise ratio larger than unity, HIST shows a clear dipole pattern (Figure 2a): positive responses within the 20°-35° latitude band, illustrating greater OHC increase in the NH latitudes relative to SH latitudes, and a negative response in the 42°-57° band. This pattern averages to zero in a hemispheric integration, and this explains why the HIST hemispheric averages are unrevealing (Figure S1 and discussions in Supporting Information S1). The temporal development of this dipole pattern shows a positive trend emerging around 1970 (Figure 2d).
AAs and GHGs drive distinct interhemispheric patterns of OHC change. For AER the OHC| ( ) responses are mainly positive (Figure 2b)-that is less heat content decrease in the NH compared to the SH. All of the ensemble means of CMIP6 models show consistent results ( Figure S2b in Supporting Information S1). While aerosol forcing is mainly distributed in the NH, it gives rise to a greater net OHC response in the SH due to the subsurface cooling pattern (Figure 1). The time evolution from AAs also features a relatively flat tail after around the year 2000 (Figure 2e), which may reflect reduced forcing from aerosol emission regulation . For GHG, differences within tropical and subtropical regions are negligible, while the hemispheric asymmetry mainly occurs in the extratropical regions poleward of 35° (Figure 2c). The negative peak at around 45° latitude is due to the large gain of heat in the Southern Ocean (Frölicher et al., 2015;Shi et al., 2021) (Figure S2c in Supporting Information S1). The GHG time series shows a rapid increase during recent decades, as expected given the atmospheric concentration history of GHGs (Figure 2f). For completeness, in the NAT runs ( Figure S2d in Supporting Information S1), there isn't any significant OHC change forced pattern.
Comparison of the results from AER and GHG with HIST indicates that the interhemispheric difference in the latitudinal structure reveals the effect of each anthropogenic forcing. For 0°-20° latitude, the poor additivity of single-forcing results indicates the nonlinear interactions of the aerosol and GHG effects on the tropical dynamics, which remains unclear. For latitudes poleward of 35°, GHG signals become more apparent associated with heat gain in the Southern Ocean. Hence, we focus on 20°-35° band, which exhibits positive signals from the HIST runs (Figure 2a), enabling a clear detection of the aerosol fingerprint.
Inspired by the above latitudinal response patterns, we examine when these signals emerge via the sensitivity to the choice of the period for calculating trends for different latitude bands. For the 20°-35° band, AER shows strong trends in OHC| ( ) after the mid-1970s (Figure 3a). In contrast, GHG only exhibits small and insignificant amplitudes for the trends with periods less than 20 years (Figure 3b). Strikingly, the results from HIST (Figure 3c) are governed by the aerosol effect (Figure 3a) in terms of these temporal features. The dominant effect of AAs during the historical period can be also found in the time series of OHC| ( ) anomalies within the 20°-35° band ( Figure S3 in Supporting Information S1) and on "end year versus start year" ( Figure S4 in Supporting Information S1). The domain in Figure 3 where trends are significant for HIST is slightly broader than the one for AER, which could be due to the additional changes driven by natural forcing. However, in the NAT alone, all trends are overwhelmed by internal variability (Figure 3d). Although insignificant, these positive trends can to some extent hinder the clear detection of aerosol signals. Again, these results indicate that more than 20 years are required to detect anthropogenic signals from background noise, and the period 1975-1990 is very important when looking to detect the effect of AAs on the OHC change.
For the 42°-57° band, HIST responses are dominated by GHG ( Figure S5 in Supporting Information S1). However, the significant domain in the start year/period space is limited to the trends of more than 50 years. It seems to be affected by the opposing but weaker trends driven by AAs within the same latitudes (Figure 2b).

Mechanism for the Forced Interhemispheric OHC Trend Patterns
In AER, the zonally integrated OHC decrease is a smooth function of latitude, larger in the SH than in NH subtropics (Figure 4 blue curve). A key question is what drives the larger net OHC decrease in the SH from AER forcing (Figures 1c and 2b)? In a heat budget analysis on the AER runs, aerosol forcing leads to substantial heat losses from the ocean to the atmosphere in the NH (as expected) but also surface heat gain in the SH ocean (Figure 4 orange curve). Most of this surface heat uptake or loss is balanced by lateral ocean heat flux divergence or convergence (Figure 4 cyan curve), and implies a strong northward cross-equatorial ocean heat flux in the subtropics and tropics, as found by Irving et al. (2019) for the longer historical period 1861-2005. Between 20° and 35°N, heat loss through the surface is only slightly larger than the anomalous OHT convergence, which results in weak negative OHC trends (Figure 4 and Figure S6 in Supporting Information S1). Within the corresponding latitudes in the SH, the anomalous OHT divergence overwhelms the surface heat gain, suggesting ocean circulation adjustments to the aerosol forcing play a critical role in driving the OHC decrease in the SH. These results are consistent with the subsurface maximum of cooling in the SH (Figure 1c).
It has been found that AA forcing drives a southward displacement of the Hadley cell and ITCZ to the warmer SH . In the models used here, we do find a zonal wind stress southward shift of the zonal wind patterns, with peaks in both hemispheres being displaced southwards (Figure 1e). This shift results in the "bowls" of subtropical wind-driven gyres also moving southwards, which produces a clear subsurface cooling along the equatorward edge of SH subtropical gyre and a warming (less cooling) for the NH gyre (Figure 1a). These features resemble the meridional gradient of mean temperature ( Figure S7 in Supporting Information S1), confirming the important role of migration. For GHG forcing, the Hadley cells are broadening in two hemispheres, moving wind stress peaks poleward, leading to a meridionally symmetric response (Figure 1f).
Ocean volume will also affect the nature of interhemispheric OHC patterns found here, since a larger volume act as a greater storage for the anthropogenic signal ( Figure S8 in Supporting Information S1). In other words, larger SH ocean volume amplifies its OHC decrease associated with the dynamic cooling. When the volume effect is removed, we find the interhemispheric asymmetry in temperature response within 20°-35° is largely reduced in AER ( Figure S9 in Supporting Information S1). Therefore, the usage of total OHC change, in which volume effect is considered, is helpful to extract the AA-driven interhemispheric patterns, which could be potentially detectable in the historical observations.

Are These Interhemispheric OHC Patterns of AA and GHG-Driven Fingerprints Detectable in Observations?
We examine observed trend estimates of 0-2,000 m OHC| ( ) in the period/start year space. The EN4_g10 and EN4_c14 products show similar trends in OHC, indicating the different methods for bathythermograph bias corrections have a small effect here (Abraham et al., 2013). The positive and significant trends in the observed OHC| ( ) (Figures 3e and 3f) match the simulated ones remarkably well. We contend that these trends are a clear and unique signature of AA-driven OHC patterns as found in the model simulations. The discrepancy between models and observations mainly occurred in trends for shorter periods (i.e., 10-20 years), which is primarily attributable to internal variability. The trends from the IAP product are smaller in magnitude compared with those from EN4 data (Figure 3g) but still qualitatively similar to that from CMIP6 runs. Importantly, we only find similar and statistically significant trends longer than 30 years starting during the 1970-1990 period to obtain signals. Positive trends vanish with starting years after the 1990s, which also agrees with the simulated results.
The observed results calculated within the higher latitude band (42°-57°) are quite different between the EN4 and IAP data sets in terms of magnitude ( Figure S10 in Supporting Information S1). The results from EN4 match well with simulations ( Figure S5c in Supporting Information S1), and trends become more negative during the latest decade, whereas the results from IAP data show much stronger negative trends during 1970-1980. The variations of temporal features between the different datasets are associated with observation bias correction methods, data collection, and gap-infilling methods (Lyman & Johnson, 2008). The large uncertainty at the high southern latitudes very likely stems from the sparseness of historical subpolar observations. Overall, the observed results are consistent with the ones from simulations in terms of the direction of these changes.  Figure S6 in Supporting Information S1). Blue shadings indicate the 35°-20°S and 20°-35°N latitude bands. 9 of 11

Implications and Discussion
Many previous studies have addressed the impact of aerosol forcing on hemispheric asymmetry of the surface temperature or OHC changes (Dallafior et al., 2015;Friedman et al., 2013;Lembo et al., 2019;Rathore et al., 2020;Rotstayn et al., 2015;Shindell et al., 2015). We show that a simple difference of entire hemispheric OHC change is not effective for detecting the aerosol effect in HIST model runs nor the observations. Here, we employ a novel approach for isolating AA and GHG effects on OHC changes during the historical period by adding back one dimension (latitude) in NH versus SH OHC differences. This reveals distinct AA and GHG effects in different latitude bands, especially for trends longer than 20-30 years.
The lack of separation in hemispheric integrations between GHG and AAs is also explained because the GHGand AA-driven OHC responses are of opposite signs in different latitude bands (Figure 2a), which then can cancel out in a hemispheric average. By differencing hemispheric responses by latitude band, we find a remarkable contrast in responses between AAs and GHGs, which we believe reflects the fact that AAs drive a strong northward OHT while the GHGs do not. It is consistent with Irving et al. (2019) and others. Associated with this AA-driven cross-equatorial ocean heat flux is a southward shift of the surface wind patterns. The latter appears to explain the asymmetry in AA ocean temperature responses, including the striking subsurface cooling around the equatorward edge of the SH subtropical gyre. The larger SH ocean volume further amplifies this AA signal. Further investigation is needed to understand this coupled response of the ocean and atmosphere to the asymmetric AA forcing, whereby the net system must supply heat (mostly from the ocean) from the SH to the NH to balance the AA impact at the top of atmosphere.
Our results also highlight that the choice of analyzed time period matters greatly to the detectability of the AA-driven signals above the natural noise and inter-model spread . We show that in the historical record, trends over 20 years or more are required to detect the human-induced responses detectable over the internal variability, and the inclusion of the 1975-1990 period is necessary to detect the AA-induced signals. The latter is consistent with previous studies focusing on shorter periods (e.g., [2005][2006][2007][2008][2009][2010][2011][2012][2013][2014][2015] being unable to detect a significant AA signal (Rathore et al., 2020).
While GHG-driven warming dominates the globally integrated OHC increase as has been widely discussed in the literature, latitudinal patterns of AA-induced OHC change we detect in model simulations are also evident in the observations. One caveat of the observed change is that although multiple data sets show the OHC change in the same direction, there are some differences in the magnitude among the datasets which likely stems from the sparseness of observations during the early period. Further studies are necessary to understand more about the differences among datasets and between observations and models.
In future projections, with AA emissions declining (Persad et al., 2023;Shi et al., 2018;Westervelt et al., 2015), the negative NH-SH OHC change within the 20°-35° latitude band and associated thermosteric sea-level change are expected. Furthermore, the geographic shifts of AA loading from North America and Europe to Asia also offer new dimensions for detecting AA-induced changes. Therefore, understanding the aerosol signals in the ocean could aid in a better representation of the OHC pattern and improvements in future projections.