Subsurface Ocean Temperature Responses to the Anthropogenic Aerosol Forcing in the North Pacific

Separating the climate response to external forcing from internal climate variability is a key challenge. While most previous studies have focused on surface responses, here we examine zonal‐mean patterns of North Pacific subsurface temperature responses. In particular, the changes since 1950 driven by anthropogenic aerosol emissions are found by using a pattern recognition method. Based on the single‐forcing large‐ensemble simulations from two models, we show that aerosol forcing caused a nonmonotonic temporal response and a characteristic zonal‐mean pattern within North Pacific, which is distinct from the pattern associated with internal variability. The aerosol‐forced pattern with the nonmonotonic temporal feature shows a substantial temperature change in subpolar regions and a reversed change on the southern flank of the subtropical gyre. A similar characteristic pattern and nonmonotonic time evolution are extracted from the subsurface observations, which likely reflect the subsurface responses to the aerosol forcing, although differences exist with the simulated responses.

• A pattern recognition method is applied to extract aerosol-induced subsurface temperature pattern from internal variability in North Pacific • Based on two large-ensemble simulations, aerosols govern a nonmonotonic feature in temperature change with a characteristic pattern • This spatiotemporal feature can be obtained from observations since 1950, although there are some differences with the modeled results

Supporting Information:
Supporting Information may be found in the online version of this article.
non-uniform across the globe, which can give rise to characteristic forced spatial patterns at the surface (Deser et al., 2020;Shi et al., 2022).An intensification of Pacific trade winds, Walker circulation, and wind-driven ocean circulation is found to be partially attributed to AA forcing since the 1990s (Allen et al., 2014;England et al., 2014;Takahashi & Watanabe, 2016).In contrast, Oudar et al. (2018) found that the observed global cooling in the early 21st century is attributable to internal variability, and the intensified trade winds are not robust in AA-forced simulations over 1998-2012.Therefore, the role of AAs in historical climate change remains unclear.
The North Pacific has attracted much attention recently in terms of the effect of AAs.AA forcing might give rise to a weakening of the Aleutian low (Dow et al., 2021;Smith et al., 2016) and a negative phase of the Pacific Decadal Oscillation (PDO) from 1981/1998 to 2012 (Dittus et al., 2021;Smith et al., 2016).Moreover, it has been found that the AA forcing from extratropics can lead to substantial changes in SST in the North Pacific extratropical regions (Diao et al., 2021;Kang et al., 2021;Luongo et al., 2022;Shi et al., 2022).In particular, Shi et al. (2022) have shown the AA-induced patterns at the sea surface based on single-forcing simulations, but the responses in the interior ocean remain unclear.Hence, using the pattern recognition method, we here focus on separating the subsurface patterns of forced ocean temperature responses in the North Pacific from its strong internal variability.The zonal-mean pattern is investigated here to remove the internal variability to some extent.

Large Ensemble Simulations
The two climate models with large ensemble simulations (more than 10 members) used in our analysis are the Community Earth System Model version 1 (CESM1) and the Canadian Earth System Model version 5 (CanESM5), in which the AA effect can be isolated from their single-forcing ensembles.
For CESM1, the all-forcing large ensemble (LENS, 40 members) and all-but-one-forcing large ensemble (XAER, 20 members) are used (Deser et al., 2020;Kay et al., 2015).Then, we use the method from Deser et al. (2020) to obtain the AA and GHG single-forcing ensemble (20 members each), respectively (Equation 1).We call these derived single-forcing ensembles the CESM1-AER and CESM1-GHG.Each ensemble member is driven by the identical single external forcing but exhibits different internal variability.
AERmember = (XAERmember − XAERensmean) + (LENSensmean − XAERensmean) (1) CanESM5 has the largest ensemble size in the CMIP6 single-forcing experiments compared with other CMIP6 models.We use 15 members from its historical aerosol-only simulation (hist-aer) and 15 members from the GHG-only simulation (hist-GHG) in which the subsurface variables are available (Swart et al., 2019).Please note that the different experimental designs used in two models, that is, single-forcing versus all-but-one-forcing, can give rise to potential differences in aerosol effects if there are substantial non-linearities in the system.In this study, we show that the results from the two models are quite consistent.The effect of the differences in experimental design merits further careful investigation.
For each model, we also use the preindustrial control simulation (1,800 years in CESM1, 1,000 years in CanESM5) to calculate the PDO index, which is defined as the leading principal component of the area-weighted SST anomalies in the North Pacific (20°N-70°N, 110°E−100°W; Deser et al., 2010).Subsurface temperature is then regressed on the PDO index to obtain the pattern associated with internal variability.
Here, we focus on the upper 1,000 m temperature change in the North Pacific, 0° to 60°N, from 1950 to 2014.We interpolate all the outputs from the two models to a regular 1° × 1° latitude-longitude grid.

Subsurface Temperature Observations
It is challenging to attribute observed climate variability to AAs on the timescale of one or two decades (Schmidt et al., 2014).Here we use the optimal interpolated EN4.2.2 potential temperature product (EN4-g10 and EN4-c14) from the Met Office Hadley Center and also the recent product from the Institute of Atmospheric Physics (IAP; Cheng et al., 2017) since they cover long periods in the subsurface observations.There are different versions with different methods to correct the bias from expendable bathythermograph (XBT) and mechanical bathythermograph measurements (Cheng et al., 2014;Gouretski & Cheng, 2020;Gouretski & Reseghetti, 2010).
The data collection, bias correction methods, and gap-infilling methods will create some differences in the objective analyses from these observationally-based datasets.Nevertheless, we focus on their common features extracted from the S/N-maximizing pattern filtering methods (see below).The data are on a 1° × 1° grid from the surface to 1,000 m depth from 1950 to 2020.

Signal-to-Noise-Maximizing Pattern Filtering
The goal of S/N-maximizing pattern filtering is to extract the anomaly patterns, for which different ensemble members agree on the temporal evolution (Allen & Smith, 1997;Delsole et al., 2011;Déqué, 1988;Schneider & Griffies, 1999;Schneider & Held, 2001;Ting et al., 2009;Venzke et al., 1999;Wills et al., 2018Wills et al., , 2020)).This method has been shown to give a clearer forced response than the simple ensemble mean given the limited ensemble member (Wills et al., 2020).In this method, the extracted patterns are associated with the maximization of the ratio of signal (e.g., variance of ensemble mean) to the total variance (from all ensemble members).Full details of this pattern-based method are given in Supporting Information S1 and Shi et al., 2022.

AA-Forced Patterns From Model Simulations
The first two S/N-maximizing patterns of CESM1-AER North Pacific zonal mean temperature anomalies and their time series are shown in Figure 1.Only for these two modes, S/Ns are much greater than 1.For the rest of the modes (not shown), the S/Ns are less than one, and the explained variances are less than 1% without showing a clear long-term temporal evolution.Therefore, we regard the first two patterns as forced patterns (FPs) driven by AA forcing.
FP1 shows a broad subsurface cooling (Figure 1a) with a quasi-monotonic time evolution (Figure 1c).The increasing rate slows down since 2000, which can be explained by the emission regulation.On the other hand, FP2 shows a broad cooling in the upper layer with a much stronger signal in the extratropical region (Figure 1b).The cooling signal mainly penetrates along the steep isothermals in the mid-latitudes.The ventilated thermocline facilitates The locations of the negative and positive anomalies relative to the mean pattern indicate that this pattern is associated with the displacement of isotherms due to AA forcing (see Section 3.2).The time series for FP2 exhibits an obvious nonmonotonic evolution, evolving in opposite directions before and after the 1980s.To examine the robustness of the results between models, we calculate the two FPs from CanESM5 AER single-forcing runs (Figure S1 in Supporting Information S1).Its second mode captures a common pattern with a clear nonmonotonic temporal feature.Please note there are some differences between the two models, such as cold anomalies in FP1 are stronger and poleward intensification in FP2 is much more muted in CanESM5, which deserve further investigation.
The regressed pattern of aerosol optical depth (AOD) on FP1 time series (Figure 1e) picks up the geographical transition effect of AA forcing (Shi et al., 2022).This AOD regression pattern shows the transition of the main emissions sites of AA from North America/Europe to East/South Asia, and also the dominant role of the emissions in Asia in the net aerosol increase during this period.The AOD regression pattern on the time evolution of FP2 shows the primary sources of AA forcing located in eastern North America and Europe, regions that underwent early industrial development and then a decline in emission-associated regulations implemented in the mid-1980s.This extratropical AA forcing in FP2 is responsible for the substantial oceanic temperature change in the North Pacific to the north of 30°N (Figure 1b).
Although FP1 explains more variance than FP2, the monotonic increase of GHG forcing can largely offset and even overwhelm the effect of FP1, making the separation of FP1 and GHG signals difficult.Hence, the effect of the nonmonotonically evolving FP2 may be more detectable in the observed historical record.In the following section, we will focus on FP2 and discuss the possible physical mechanisms to explain its characteristic zonal mean pattern.

Physical Mechanisms Associated With Extratropical Aerosol Forcing
We regress the wind stress anomalies and their curl on the time series of FP2 (Figure 2).The trade winds to the north of the equator and westerly winds are both enhanced relative to the mean wind fields (Figure 2a).This zonal mean wind stress increase gives rise to anomalous Ekman downwelling within 10°-35°N (Figure 2b), driving broadly distributed downwelling wind stress curls in the mid-latitudes (Figure 2d).Consistent with previous studies, a clockwise meridional atmospheric overturning circulation response over the tropics is evident, which is a strong signal of a southward shift of the Hadley Cell (Figure 2f), in response to the NH cooling driven by aerosol forcing.This shift is consistent with the enhanced trade winds to the north of the equator (Figure 2a) and also gives rise to a southward shift of the wind-driven ocean gyres.The southward shift of tropical surface winds and associated enhanced Ekman downwelling together contribute to the subsurface warming on the equatorward flank of the subtropical gyre shown in Figure 1b.It is worth noting that these regression patterns reflect the nonmonotonic evolution of the second mode.Thus, the responses are in the opposite direction after the mid-1980s.
The regression of surface heat flux on the FP2 time series shows a widespread heat loss in the North Pacific (Figure 2c).The heat loss is substantial in the subpolar regions, at similar latitudes to the aerosol sources (Figure 1f), which contributes to the poleward enhanced cooling in the North Pacific (Figure 1b).This poleward enhanced cooling increases the meridional SST gradient, which further affects the tropospheric responses (e.g., Xu & Xie, 2015).The regression of the zonal mean of zonal atmospheric circulation over the North Pacific shows an enhancement of mid-latitude jet streams and a clear poleward shift in the upper troposphere.These atmospheric responses are quite consistent between the two models (Figure S2 in Supporting Information S1).Thus, a poleward shift of the westerly jet can be driven by extratropical aerosol forcing in North America and Europe, likely through the strengthened meridional SST gradient.In this nonmonotonic mode, the responses of the Hadley Cell and N. Pacific jet in the upper troposphere to AA forcing show opposite directions, which is consistent with the result from Diao et al. (2021).They concluded that the shift of the jet stream is attributed to aerosol forcing from South/East Asia, rather than from North America and Europe.In this study, there is a tremendous temperature drop around 40°-45°N (Figure S1b in Supporting Information S1), which is to the north of the jet core (∼30°N).Therefore, based on the thermal wind relation, the enhanced SST meridional gradient in North Pacific can produce a strengthening of zonal winds whose maxima are located to the north of the jet core.

Results From Observations
As the mid-1980s are the turning point of AA-induced FP2, we calculate the observed trends of North Pacific zonal mean temperature for 1950-1985 and 1985-2020, respectively (Figure 3).The former 36-year trend shows a significant cooling spot in the mid-latitudes, which penetrates mainly along the isotherms (Figures 3a, 3c,  and 3e).In addition, the subsurface warm anomalies are found on the equatorward flank of the subtropical gyre, similar to the AA-induced FP2 from the model simulations.For the latter period, the strong warming in the upper layer dominates the overall zonal mean pattern with weaker subsurface cooling on the equatorward flank of the subtropical gyre (Figures 3b,3d,and 3f).The contrasting trend patterns also appear in the AA single-forcing runs, while the deviations across the ensemble members reflect the impact of background noise (Figures S3 and S4 in Supporting Information S1).
To further show the time evolutions of the observed subsurface temperature, we choose two boxes to calculate the averaged temperature anomalies (Figure 3a).The temperature difference between the two boxes shows a clear nonmonotonic evolution in all three datasets (Figure 3g).The -1xPDO index also exhibits a similar long-term variability but is not able to fully capture the subsurface temporal evolution, especially for the first and last 10 years.We would like to note that this observed PDO can include both the external-forced and internal components.
The results from CESM1 simulations for the same boxes are shown in Figure 3h.The GHG effect is relatively  monotonic (red curve in Figure 3h).Importantly, we find that the nonmonotonic feature is predominantly driven by AER2 (reconstructed temperature from FP2 based on Equation S6 in Supporting Information S1; cyan curve in Figure 3h).Hence, the existence of AER2 stemmed from the AA forcing in NH extratropical region can to some extent affect the nonmonotonic trend pattern discovered in observations (Figure 3g).
Although the nonmonotonic feature can be detected, the magnitudes of piecewise trends shown in Figure 3 do not match the AA-induced FP2 well, especially over 1950-1985.To increase the signal-to-noise, we also apply the pattern recognition method to the three observational datasets to extract the common features from them.In this case, it is gridding and bias corrections differences and not internal variability, which is treated as noise.In addition, we expect the observations to reflect both internal variability and all external forcings together.
As we focus on the nonmonotonic feature here, the zonal-mean temperature fields are detrended during 1950-2020 before the pattern recognition analysis.The first mode (EP1) shows strong warming in the upper layer with a poleward amplification and a relatively weak cooling on the southern flank of the subtropical gyre (Figure 4a).The time series of EP1 shows piecewise, long-term trends before and after the 1980s (Figure 4c), which resembles the time series based on the AER2 in Figure 3h.The second extracted mode (EP2) shows a clear tripole pattern (Figure 4b), and the corresponding time series consists of a recognizable multidecadal variability.The S/N ratios in these two modes are much larger than 1.For the remaining modes (not shown), no clear long-term signal is found.The pattern correlation between FP2 and EP1 (Figure 1b vs. Figure 4a) is relatively high (r = 0.80 for CESM1, r = 0.64 for CanESM5) compared with the correlations between FP2 and the other extracted observed patterns (Table S1 in Supporting Information S1).Hence, EP1 potentially reflects the aerosol fingerprint.In addition, we also conduct the low-frequency component analysis based on the method from Wills et al. (2018), and find the results are quite consistent with the SST responses from their work (Figure S5 in Supporting Information S1).While it is remarkable that the observed zonal-mean temperature trends in the North Pacific bear similarity to the AA-induced FP2 determined from the model simulations, it is unclear how much the observed multidecadal trends reflect internal variability, especially associated with the PDO.To show the subsurface pattern associated with internal variability alone, we regress the zonal mean temperature from the pre-industrial control run on its PDO index (Figure S6 in Supporting Information S1).The regression shows a tripole pattern that is broadly similar to EP2 from observations.Responses from the control run are more intense in extratropical regions as the PDO index is obtained from 20°N-70°N (Deser et al., 2010).The pattern correlation between the PDO pattern from models and EP2 is relatively high (r = 0.64 for CESM1, r = 0.63 for CanESM5; see Table S2 in Supporting Information S1).Hence, EP2 is likely to reflect the unforced PDO variability.

Summary and Discussion
In this study, we have extracted the influence of anthropogenic forcings on the North Pacific subsurface ocean temperature change since 1950 based on the single-forcing large ensembles of simulations.Despite the confounding effects of internal variability and GHG-forced signals, the AA-forced signals can be extracted using an S/N-maximizing pattern filtering method.Our analysis indicates that aerosol emissions from North America and Europe can cause a clear nonmonotonic time evolution and a characteristic zonal mean pattern of the subsurface ocean temperature within the North Pacific.This response is found in two different climate models with large ensemble simulations, and can be used to understand the AA-induced fingerprints in observations.Time series analyses, such as that of global/regional mean SSTs, are limited in being able to separate the forced signals from internal noise.The situation is even worse when the observed record is sparse and internal variability is strong and complex.In this study, the ocean subsurface temperature pattern of forced response provides additional information to better characterize the difference between internal modes and changes driven by anthropogenic forcing.For instance, FP2 from the AA simulations is distinct from the PDO-related pattern from the preindustrial run.There are substantial surface heat flux changes in the North Pacific subpolar regions and also tropospheric circulation changes, which leave "fingerprints" in the interior ocean: The increasing surface heat loss in the subpolar regions (prior to ∼1985) leads to strong ocean surface cooling which then subducts into the thermocline.The surface wind change, associated with the equatorward shift of the Hadley cell, strongly affects the temperature change on the equatorward flank of the subtropical gyre.Importantly, after ∼1985, these changes reverse.
Disentangling the roles of internal variability and external forcing in historical climate change remains challenging (e.g., Hua et al., 2018;Qin et al., 2020).We have focused on the nonmonotonic temporal evolution of AAs and their induced zonal mean pattern of subsurface ocean temperature.Here we propose a way to use the obtained FPs from models to compare with the observed changes, although the observed patterns are still open to interpretation.As a caveat, there are some differences between the patterns obtained from models and observations on the regional scale, such as the weaker observed subsurface response at around 500 m, which deserves further investigation.Another useful extension of our work should include the pattern-based fingerprint analysis (e.g., Santer et al., 2018) to understand potential interactions between GHG and AA forcing and how these show up in observations versus models.For the atmospheric circulations, the Hadley cell and jet over the North Pacific shift in opposite directions.The shift of the Hadley cell is attributable to the interhemispheric energy imbalance, and the shift and enhancement of the jet stream are found to be associated with the tremendous temperature change in the North Pacific subpolar region to the north of the jet core.These results are consistent with Diao et al. (2021) but are found to be induced by aerosol emissions from North America and Europe.Rather than the monotonic effect of increasing aerosol emissions widely discussed in previous studies, we show that these atmospheric circulation responses have a clear nonmonotonic feature that further affects ocean responses.

Figures S1 to S6
17 Tables S1 and S2 18 and form the ensemble data matrix X, in which the data from each ensemble member are 25 concatenated in the time dimension.We then conduct the EOF analysis on the matrix X.The 26 matrix is weighted by the square root of the grid cell area.We define linear combinations of 27 eigenvectors (ak) from the EOF analysis: 28 ) where uk is the pattern from the linear combination of N leading EOF eigenvectors. ! is the 30 square root of the i-th EOF eigenvalue, and ek is a coefficient vector that is determined as that 31 which maximizes the S/N ratio as described below.In this analysis, we pick N = 50 retaining 32 around 99% of the total variance.

Figure 1 .
Figure 1.The first two forced patterns (FPs) and time series from CESM1 AER simulations.(a, b) FPs (FP1 and FP2) of North Pacific zonal mean temperature from 1950 to 2014.The mean temperature is shown as contours.(c, d) Standardized time evolution for FP1 and FP2, respectively.(e, f) Regressions of aerosol optical depth on each time evolution.

Figure 2 .
Figure 2. Regressions of various fields on the FP2 time evolution in Figure 1d.(a, c) Regressions of zonal wind stress, wind stress curl, and surface heat flux anomalies on FP2 time evolution from CESM1-AER (black curves) and CanESM5-AER (red curves).(a, c) show the zonal mean of regression and panel b shows the zonal sum of regression of wind stress curl within the North Pacific (0°-60°N, 100°E−100°W).Negative values in panel c denote heat losses from the ocean.The mean fields are shown as dashed curves.(d) Map of regression of surface wind stress (arrows) and wind stress curl (shadings) from CESM1-AER.(e, f) Regression of North Pacific atmospheric zonal mean velocity and global meridional overturning streamfunction, respectively, on FP2 time evolution from CESM1-AER.Positive values in panel f denote clockwise circulations.The corresponding mean fields are shown as contours.

Figure 3 .
Figure 3. Changes in North Pacific temperature from observations and CESM1 simulations.(a-f) 1950-1985 and 1985-2020 trends of zonal mean temperature from observations.Climatological isotherms are shown as contours.Stippling indicates regions exceeding 95% statistical significance from the two-tailed t-test.Panel (g) The 5-year running mean of temperature difference between the average temperature in 25°N-45°N, 50-300 m (green box in panel a and 12°N-25°N, 200-600 m (magenta box) from observations.The 5-year running mean of -1xPDO index from the WMO Regional Climate Center is shown as the blue curve.(h) Temperature difference between two boxes in panel a from CESM1-AER (AER1 and AER2 from FP1 and FP2, respectively) and CESM1-GHG (GHG1 from its first mode).Dashed curves in panel h are the results of a simple ensemble mean of subsurface temperature, <X>.

Figure 4 .
Figure 4. First two extracted patterns and time evolutions from three observational datasets.(a, b) Extracted patterns of North Pacific zonal mean temperature from 1950 to 2020.Climatological isothermals are shown as contours.(c, d) Time evolutions for the extracted patterns.The thick curve denotes the 10-year running mean.

First, we calculate
the zonal mean of the annual mean subsurface temperature in the 20 North Pacific for each ensemble member.The temperature anomalies are calculated relative to 21 the mean of 1950-2014.The obtained temperature anomalies are divided by the temporal 22 standard deviation at each depth for the respective ensemble member to give equal weight to the 23 temporal variability vertically in the following Empirical Orthogonal Function (EOF) analysis, 24

Table S2 .
Uncentered pattern correlations (absolute values) between the unforced PDO regression 104 pattern from the preindustrial control simulations (PI) and the extracted patterns (EPs) from