Shortwave Radiative Flux Variability Through the Lens of the Pacific Decadal Oscillation

The variability of the shortwave radiative fluxes at the surface and top of atmosphere (TOA) is examined in a pre‐industrial modeling setup using the Pacific Decadal Oscillation (PDO) as a possible pacemaker of atmospheric decadal‐scale variability. Within models from the Coupled Model Intercomparison Project—Phase 6, downwelling shortwave radiation at the surface, the net shortwave fluxes at the surface and TOA, as well as cloud radiative effects show remarkably similar patterns associated with the PDO. Through ensemble simulations designed with a pure PDO pattern in the North Pacific only, we show that the PDO relates to about 20%–40% of the unforced year‐to‐year variability of these shortwave fluxes over the Northern Hemispheric continents. The sea surface temperature imprint on shortwave‐flux variability over land is larger for spatially aggregated time series as compared to smaller areas, due to the blurring effect of small‐scale atmospheric noise. The surface and TOA radiative flux anomalies associated with the PDO index range of [−1.64; 1.64] are estimated to reach up to ±6 Wm−2 for North America, ∓3 Wm−2 for India and ±2 Wm−2 for Europe. We hypothesize that the redistribution of clouds in response to a North Pacific PDO anomaly can impact the South Pacific and North Atlantic SSTs.


Introduction
Internal climate variability on a range of spatial and temporal scales and its interplay with radiative processes is important for understanding Earth's energy balance and its potential response to changing forcings over different temporal and spatial scales (e.g., Dessler et al., 2018;Lehner & Deser, 2023;von der Heydt et al., 2021).The shortwave part of Earth's energy budget, specifically absorbed solar radiation (or net shortwave flux at the top of atmosphere (TOA) and its response to greenhouse gas forcing are shown to be fundamental for Earth's climate and climate change (Donohoe et al., 2014;Trenberth & Fasullo, 2009).
The variability of the global mean Earth radiative imbalance on interannual time scales has been linked to the sea surface temperature (SST) variability in the Pacific ocean, specifically the El Niño-Southern Oscillation (ENSO) (Allan et al., 2014;Ceppi & Fueglistaler, 2021;Fueglistaler, 2019;Loeb et al., 2018;Wills et al., 2021).Decomposing this relationship into shortwave and longwave components, the variability in the global mean is mostly related to changes in absorbed shortwave radiation through changes in clouds and sea ice (Loeb et al., 2018(Loeb et al., , 2021)).Investigating the peculiar difference in trends before and after 2014, (Loeb et al., 2018), show that the spatial distribution of the anomalous absorbed solar radiation matches the SST pattern of the Pacific Decadal Oscillation (PDO), which is the focus of our study.
The PDO is identified as the dominant pattern of SST variability in the North Pacific, characterized by decadalscale warming and cooling with a characteristic horseshoe pattern (Mantua & Hare, 2002).On decadal time scales, the PDO is associated with strengthening and expansion of the North Pacific subpolar gyre in response to a deepening of the Aleutian Low and increasing variability at longer time scales due to adjustment of westward propagating oceanic Rossby waves (Qiu et al., 2007;Taguchi et al., 2007;Wills et al., 2019)).The PDO has a South Pacific counterpart, referred to as South Pacific Decadal Oscillation (SPDO, Chen & Wallace, 2015;IPCC, 2021), which is thought to be influenced by a variety of processes including extratropical modes, ENSO teleconnections and ocean dynamics (Shakun & Shaman, 2009;Zhang et al., 2018).The combined phenomenon encompassing the PDO, ENSO, and SPDO, though believed to be physically distinct modes (IPCC, 2021;Newman et al., 2016), is the Inter-decadal Pacific Oscillation (IPO, Folland, 2002;Henley et al., 2015;Power et al., 1999).The PDO is found to be related to atmospheric dynamics, specifically the Pacific North American pattern, and is strongly correlated with temperature and precipitation patterns over North America (Liu et al., 2017;Rohli et al., 2022).There are further indications of a possible link between the Pacific North American Pattern and North Atlantic Oscillation and Atlantic storm tracks through baroclinic waves (Pinto et al., 2010).
Links between the PDO and Earth's radiative balance have also been established in the literature.Investigating observed trends in absorbed solar radiation between 2005 and 2020, Loeb et al. (2018Loeb et al. ( , 2021) ) show that the strong increase in global mean shortwave anomalies after the shift from a negative to positive PDO in 2014 is mainly due to a decreased low cloud amount in the warm sector off the coast of North America.Over the continental U.S., the PDO has also been proposed as a major driver of decadal-scale surface downwelling shortwave radiation anomalies with negative radiation anomalies during a positive phase and vice versa (Augustine & Capotondi, 2022).Within climate models, the PDO has been linked to decadal cloud feedbacks and the temporal variability of the climate feedback parameter (Meyssignac et al., 2023;Zhou et al., 2016).
Records from surface observations dating back to the 1930s reveal decadal variability of the downwelling shortwave radiation at the surface (dimming and brightening), which has been linked to the temporal variability of clouds and aerosols (e.g., Stanhill & Moreshet, 1992;Sanchez-Lorenzo et al., 2015;Wild, 2009;Wild et al., 2013Wild et al., , 2015)).Taking a climate model perspective, Chtirkova et al. (2023) establish a connection between the PDO and surface shortwave radiation anomalies both above ocean and land, proposing it as an important driver for continental dimming and brightening not only for the U.S. (Augustine & Capotondi, 2022) but also for Europe and Western Asia.
In this study, instead of inferring the possible causes of radiative anomalies, we take the PDO as a starting point and ask the question what are the possible radiative anomalies that are associated with this single mode of variability.As our interest is with the PDO as a potential pacemaker of part of the internal variability of shortwave fluxes on decadal scales, we use the PDO evolution as a reference and examine the shortwave fluxes "through the lens of the PDO."We do so by analyzing both coupled climate model simulations from the Coupled Model Intercomparison Project-Phase 6 (CMIP6, Eyring et al., 2016) and our own simulations with the atmosphereonly part of the icosahedral nonhydrostatic Weather and Climate Model (ICON), where we constrain the PDO spatial anomaly to a specific value and assess the response of the shortwave radiative fluxes.We focus on various aspects of the shortwave part of Earth's energy budget, including the downwelling shortwave radiation at the surface-F ↓ S , net shortwave radiation at the top of the atmosphere (TOA, positive downward)-F ↓↑ T , net shortwave radiation at Earth's surface (positive downward)-F ↓↑ S , shortwave atmospheric absorption-A atm , and the shortwave cloud radiative effects at the surface and TOA-CRE S , CRE T .The paper is structured as follows: in Section 2 we describe the data and experimental setup.The results are presented in three parts: in Section 3.1, we identify patterns of decadal trends in the shortwave fluxes; in Section 3.2, we estimate the fraction of total variability attributable to SSTs and the PDO and in Section 3.3 we quantify the anomalies in the radiative fluxes as a function of the PDO index value.In Section 4, we discuss potential differences that arise from the lack of coupling or from using specifically ICON-A; we also compare the model-based PDO anomalies used in this study to observations.We conclude in Section 5.

Data and Methods
The study consists of two parts-we use coupled model data to pinpoint the qualitative relationship between the shortwave fluxes and the PDO and we give a quantification of the fraction of variability attributable to the PDO and the flux anomalies in Wm 2 based on our own numerical simulations.

CMIP6 Data
The first part of the study uses the unforced control simulations (piControl) of CMIP6 (Eyring et al., 2016) to investigate the unforced variability of the shortwave fluxes at the surface and TOA and relate them to the PDO.The PDO index for the CMIP6 analysis is computed using the Climate Variability Diagnostics Package-CVDP, version 5.1.1 (Phillips et al., 2014), developed by NCAR's Climate Analysis Section.We include 52 simulations from 44 coupled models.The multiple simulations for some models come from slightly modified model versions in terms of physical parameterizations used or as different realizations.The full list is given in Table 1.Our analysis uses annual mean data, and is performed on a per grid box level.The radiative fluxes are interpolated to a 1°grid (same as the GFDL-ESM4 grid) using 2nd order conservative remapping (Jones, 1999).This part is an expansion of Chtirkova et al. (2023) to all shortwave fluxes, and the analysis is done in the following order.
1.For each of the 52 simulations, we compute the distribution of all possible trends (linear regressions) in the annual mean PDO index.2. We take the periods with trend magnitudes below the 10th and above the 90th percentile of the PDO index.3. We combine them into composite (mean) trend maps of the shortwave fluxes-one for the increasing and one for the decreasing phases of the PDO. 4. The individual model maps from steps 1-3 are subsequently combined into a multi-model median map per shortwave flux component.
In addition, to provide a region-specific discussion, we select four regions by visual inspection of the shortwave flux trend maps (that are derived from the PDO index trends), focusing on regions that do show trends upon strong changes of the PDO.We choose regions with approximately equal area (between 3.7 and 3.9 10 6 km 2 ).The regions we come up with are the following and we refer to them with the names of the nearest (but not the whole) geographical features they cover: North America (30-50°N, 250-270°E), Europe (40-55°N, 0-30°E), India and Indochinese Peninsula (10-18°N, 71-110°E) and Australia (15-30°S, 128-150°E).A visual depiction of the regions is provided on the last subplot of Figure 2.

ICON-A Simulations
We use the atmospheric part of the ICON Earth System Model (ESM): ICON-A (Jungclaus et al., 2022).The atmospheric general circulation model ICON-A (Crueger et al., 2018;Giorgetta et al., 2018) is in a configuration using the Max Planck Institute climate physics package on the rotated R2B4 grid (grid id 0013) with an approximate grid increment of 160 km.The grid is the same as in the ICON ESM-Low Resolution (ICON-ESM-LR), which is the ICON participation as a coupled model in CMIP6.The tuning configuration of our ICON-A instance is also the same as in the coupled version's (ICON-ESM-LR) piControl simulation.As boundary conditions, we use monthly mean SST data from the first 100 (out of 500) years coupled simulation within CMIP6.
The SST climatological mean fields for the simulations as well as sea ice climatologies are computed over all 500 years of ICON-ESM-LR data.All other boundary conditions are kept at pre-industrial climatological values, including sea ice data, stratospheric ozone and aerosol optical properties.Well-mixed greenhouse gasses have a constant vertical mixing ratio and the solar constant is kept at 1,360.744Wm 2 in accordance with CMIP6 protocols.Our motivation for using SST data from the coupled CMIP6 run instead of observed data is two-fold: (a) we have 500 years of statistics that are not contaminated by global warming and (b) the SST field is in part characteristic of the atmospheric circulation.By using SSTs from the coupled counterpart of the same model, we avoid a possible mismatch between the SST data and the specifics of the atmospheric circulation of the ICON model and its R2B4 grid.
We run experiments with prescribed SSTs, where we retain only the low-frequency component of the PDO.This is done by first decomposing the SST field (500-year of ICON-ESM-LR) into eigenvectors (empirical orthogonal functions, EOFs), eigenvalues and their principal components (PCs).Following the recommendations of NCAR Climate Data Guide (Schneider et al., 2013), we deseasonalize and decompose the SST field in a rectangular box in the North Pacific.The box boundaries are taken from the CVDP source code (which we use for the CMIP6 analysis).The empirical orthogonal functions (EOF) decomposition is done using the Python package: pyEOF (Zheng, 2021), which uses the Numpy singular value decomposition methods (Harris et al., 2020).The PDO is defined as the first principal component (PC) (PC 1 , the one with the largest eigenvalue) of the decomposed field-Figure 1 shows the eigenvector expressed as correlation (a), its corresponding time series (b) and power spectra (c).We further low-pass filter the PC time series using a logarithmic function in frequency space with a time window of 10 years.To obtain a complete SST field as model boundary conditions, we re-project the resulting low-frequency time series (lfPDO) on the PDO eigenvector and add the per-grid box seasonal cycle on top.We do not sum up the remaining eigenvectors and PCs, giving us an SST field, where the only inter-annual variability comes from the PDO pattern.All grid boxes outside the PDO box are kept at their climatological values, apart from a smoothing distance around the box which is included for numerical considerations and detailed in Appendix A.
To better carve out variability related to SST variability in general and SST variability related to the PDO in particular, we run two main sets of simulations with 20 ensemble members each: a set of 100-year simulations where the SST field is taken as is from the first 100-year of the piControl run ("allSST") and another set of 100year simulations where the only variability in the SST field comes from the low frequency component of the PDO obtained as described above via low-pass filtering and projecting the PC 1 time series ("lfPDO").Ensemble members share the initialization of the atmosphere but differ in the calendar day and associated insulation at simulation start (e.g., January 1, January 5 etc.).We discard the first 1-5 years after initialization, considering this as model spin up, and analyze only subsequent data.This data corresponds to years 4006-4099 (year numbers as they appear in the CMIP6 data archive) of the ICON-ESM-LR piControl run.
To further quantify the influence of a concrete PDO phase (index value) on the shortwave fluxes, we run seven 100-year simulations with a "perpetual" PDO phase with different amplitudes of the SST anomaly and analyze the last 95-year to reduce the impact of the initial condition.To obtain the field, we take the 99th, 95th and 68th percentiles of the distribution of monthly PDO index values from the 500-year ICON-ESM-LR CMIP6 piControl simulation, and project constant time series with these index values onto the PDO eigenvector, seasonality is added on top.The seven PDO index values we consider based on the above mentioned percentiles and their negative counterparts are 2.58, 1.64, 0.44, 0, 0.44, 1.64 and 2.58.Sea surface temperature anomalies that correspond to these values are further detailed in Appendix A.

Results
We examine the impact of the PDO on the various shortwave flux components.Each subsection takes a different perspective in terms of time scales (decadal scale trends and year-to-year variability) and with regard to the details of the PDO.Section 3.1 looks at the PDO arising in fully coupled CMIP6 piControl simulations, where it combines with other modes of variability, and asks about the imprint of strong PDO phase changes on decadal    ).Panel (c) displays the power spectral density estimated using Welch's method (Welch, 1967) for the entire 500year simulation.

Patterns Related to PDO
Our first quest is to investigate the patterns of the shortwave radiative flux components related to PDO phase shifts within CMIP6 models.The analysis is done for the distribution of all possible 20-year trends (linear  regressions) of the PDO index, selecting the periods with strongest transitioning of the index and computing the trends of the energy balance components for those periods (same as in Chtirkova et al. (2023) for F ↓ S ).The strongest transitions of the PDO index in our case correspond to 20-year trends above 90th and below 10th percentiles for PDO↑ and PDO↓ transitions respectively.The results are presented in Figures 2 and 3 for all-sky and clear-sky fluxes, respectively.For simplicity, we only discuss results related to PDO↑ in the text, that is, the PDO changing phase from negative to positive, bearing in mind that the same is valid with opposite sign for the opposite transition.
For the all-sky shortwave fluxes, we find that the patterns are remarkably similar both in spatial structure and in magnitude for F ↓ S and the net shortwave fluxes at the surface and TOA (F ↓↑ S , F ↓↑ T ).This suggests that the reduction of F ↓ S (in N. America and Europe during PDO↑) is also associated with a reduction in F ↓↑ T and in F ↓↑ S , that is, less shortwave radiation absorbed over these regions by the whole system and in particular-at Earth's surface.A significantly different pattern is displayed by shortwave atmospheric absorption (A atm ).Regions bearing a statistically significant imprint (no hatching) of the PDO phase shift in A atm may not show up in the other shortwave fluxes and vice versa (e.g., parts of Africa or the Americas).Regions featuring a statistically significant impact from PDO phase changes may show opposite signs for A atm than for the other shortwave fluxes.This is notably the case in wide parts of the Pacific and neighboring regions.For example, the increase in F ↓↑ T over Australia combines with a reduction in A atm to result in an even stronger increase in F ↓ S downward.
On the clear-sky side (Figure 3), patterns are overall least pronounced for F ↓↑ T,cs .The much more pronounced and geographically more extended patterns in F ↓ S,cs and A atm,cs follow each other closely but with opposite sign, suggesting that changes in F ↓ S,cs are primarily mediated by changes in A atm,cs .The absorbed clear-sky shortwave radiation at the surface (F ↓↑ S,cs also resembles the F ↓ S,cs and A atm,cs for most ocean areas except for the parts where sea ice plays a role.The clear-sky F ↓↑ T,cs in the Arctic regions increases with PDO↑, which is associated with less reflectivity in the system, likely due to the decreased sea ice area during a warm PDO (Simon et al., 2022).This surface albedo effect is also evident in the all-sky patterns.The difference in the all-sky A atm and clear-sky A atm,cs implies that increased absorption in cloudy conditions due to multiple scattering is significant for regions like the North Pacific.
CREs at the surface and TOA, associated with the different PDO phases, also yield identical patterns and similar magnitudes to F ↓ S , F ↓↑ T and F ↓↑ S , with the exception of the Arctic.The increase of F ↓ S in the warm North Pacific pool and the decrease in the cold pool further serve to enhance the SST anomalies associated with the PDO.An interesting feature is the shift from enhanced F ↓ S , F ↓↑ T and F ↓↑ S in the warm region off the west coast of North America during a PDO↑ to trends of the opposite sign over North America, indicating that a PDO↑ leads to dimming over the continent.
We assess the uncertainty of the presented results on a grid box level in two ways: we test whether the samples assembled upon results from individual models for PDO↑ (first sample) and PDO↓ (second sample) are part of the same distribution using the Mann-Whitney-U test (Mann & Whitney, 1947); we also compute the spread among models (difference between the 90th percentile and 10th percentile values) for each grid box.Depending on the variable, using the Mann-Whitney-U test, the null hypothesis that the two samples are part of the same distribution cannot be rejected at the 0.05 level in 27%-39% of grid boxes (34% for F ↓ S ), that is, for these grid boxes there is no statistically significant difference between periods of strongest negative or positive changes in the PDO index.These regions are hatched in Figures 2 and 3.The CMIP6 multi-model spread across grid boxes is on average 0.12 Wm 2 year 1 for the all-sky fluxes and CRE (apart from A atm ) and on average 0.02 Wm 2 year 1 for the clear-sky fluxes and all-sky A atm (given are median values across grid boxes representative for all variables).Model spread also depends on the region and is highest in regions where the trend values are the highest.In relative units, the spread in clear-sky fluxes across models is higher in clear-sky as compared to all-sky, which was also highlighted in Chtirkova et al. (2022).
The differences among models with regard to the PDO and its relation to shortwave fluxes may be illustrated by inspecting the model specific correlations between the PDO index annual mean time series and regionally averaged annual mean shortwave flux time series, even though they cannot be directly related to the decadal scales discussed above.Corresponding Pearson correlation coefficients are given in Table 1.We choose the four regions (see last map in Figure 2): North America and Europe, where our trend maps show a negative correlation for the all-sky fluxes (apart from A atm ) and India and Indochinese Peninsula and Australia where the maps show a positive relationship.For F ↓ S , F ↓↑ T and F ↓↑ S , almost all models agree on the sign of the relationship in these regions.Exceptions include EC-Earth3, INM-CM4-8, INM-CM5-0 and NorCPM1.Exceptions for A atm , occur more often but trend patterns are also different for A atm as compared to F ↓ S and the choice of regions is not optimal for A atm .We also highlight the models which show the strongest correlation coefficients between the PDO index and shortwave fluxes in North America, India and Australia: CESM2-FV2, CESM2-WACCM-FV2, CMCC-CM2-SR5, CMCC-ESM2, GFDL-ESM4, GISS-E2-1-G (p5), ICON-ESM-LR and MIROC6 (given are models showing a correlation coefficient above 0.5 for at least one region and variable).The weaker correlations for Europe are expected due to its larger distance from the Pacific Ocean and other causes for variability there but models still agree on the sign of the relationship.
Going back to trend magnitudes and Figures 2 and 3, we mark as important that even though trend magnitudes are given in absolute units, they are strongly dependent on model and percentile of the PDO index trends over which we select periods.Trend magnitudes are on average 23%-28% stronger (depending on variable) when increasing the percentile range from 90-10 (as shown in Figures 2 and 3) to 95-5.This happens for a few reasons: (a) PDO trends across models differ and (b) PDO trends in the same model above the 90th percentile differ, (c) the atmospheric response to the same PDO trend across models might be different.Taking a smaller subset of PDO indices does not necessarily provide better statistics because of too much non-PDO related variability that is not averaged across many data points.To further quantify what fraction of the variability is attributable to the PDO and what are the specific radiative anomalies related to its phases, we rely on results from simplified numerical experiments described in the following sections.

PDO Fraction in Variability
In our simplified ICON-A numerical setup, we would like to give an estimate on what fraction of total shortwave flux variability can be attributed to the SST field and what fraction of that is attributable to the PDO pattern.Our whole analysis is based on annual mean data and, as the results suggest, a large fraction of it is still due to variability intrinsic to the atmosphere.A large influence of atmosphere-only on inter-annual variability of temperature and sea level pressure is also shown in (Lehner and Deser (2023)).We attempt to separate variability into sources intrinsic to the atmosphere and sources related to the SST field by running ensemble simulations.Firstly, we separate atmosphere-only variability from SST-related atmosphere-ocean variability and as a second step we carve out only variability related to the time-varying PDO pattern within the SST field.
For the first part, we run 20 ensemble members with SSTs directly taken from the coupled run ("allSST"): one ensemble member of this experiment (our control run) should possess both atmosphere and coupled variability.The atmosphere-only variability should be averaged out in the ensemble mean time series (average time series across ensemble members for each grid box) and the only remaining variability in the 95-year ensemble mean should be related to the prescribed time-varying SSTs, as they are the only time-evolving boundary condition in our "allSST" experimental setup.
Taking this idea further, annual mean time series averaged over the 20 ensemble members of the "lfPDO" experiment should retain only variability associated with the low frequency PDO variability.The standard deviation of the latter (computed upon the 95-year ensemble mean time series per grid box) is shown in Figure 4, first column.We also show results for the surface temperature T S (last row in Figure 4) because they best illustrate the direct effect of the experimental setup.Testing for a different number of ensemble members, we find the resulting reduction of variability converges for around 10 ensemble members.
To gain an impression of the relative importance of SST variability in general and PDO variability in particular, we compute the differences relative to the control run of the standard deviations of the time series per grid box as: σ allSST,ensmean /σ allSST,ens01 , where σ allSST,ensmean is the per grid box standard deviation of the ensemble mean time series (atmosphere-only variability is averaged out), and σ allSST,ens01 is the per grid box standard deviation of one of the ensemble members (our control run; includes atmospheric variability; results are not sensitive to the choice of specific ensemble member).The ratio is shown in Figure 4, second column.We see that SST-related variability is responsible for almost all of total (atmosphere and ocean) variability in the tropical regions.The standard deviation for the ensemble mean is reduced in roughly half for the shortwave flux variables for continental regions like North America and Europe (on a grid box level).Looking into spatially averaged (aggregated) time series instead of individual grid boxes (i.e., taking the spatially averaged time series and computing its standard deviation instead of , ratio between the standard deviation of the 20-ensemble-member mean in the "allSST" experiment and the standard deviation for the control run (second column) and ratio between the standard deviation of the 20-ensemble-member mean in the "lfPDO" experiment and the control run (third column).The control run is one ensemble member of the "allSST" experiment.
averaging the standard deviations of the time series per grid box), the resulting ratio σ allSST,ensmean /σ allSST,ens01 is 10%-30% more (not shown), that is, SSTs explain a larger fraction of variability for spatially averaged time series across regions than for individual grid boxes.This is not surprising because the influence of the SSTs over land usually occurs thorough the large-scale circulation, the effect of which is reduced at the grid-box level by smallscale noise and this small-scale noise is already averaged out in the aggregated time series, yielding a smaller reduction in variability.For the aggregated time series in our regions of interest, we find that the ratio is 0.78 for F ↓ S and 0.69 for F ↓↑ T for North America, 0.53 for F ↓ S and 0.54 for F ↓↑ T for Europe and 0.74 for both F ↓ S and F ↓↑ T for India and Indochinese peninsula.The fraction of variability of CRE that can be related to the time evolving SSTs is 0.62 of the variability in the Northern Hemisphere (0.57 for North America, 0.45 for Europe and 0.67 for India and Indochinese peninsula).
The second step is to estimate how much of the total variability is related to the low-frequency PDO pattern.We do so by computing the ratio between the 20-ensemble-member mean in the "lfPDO" experiment and the control run: σ lfpdo,ensmean /σ allSST,ens01 .The result (third column in Figure 4) is that the low frequency component of the PDO in our experimental setup is related to almost none of the variability in the tropics for the variables we investigate, it is not, or little, related to the SST-related variability in Australia and South America.For the Northern Hemisphere aggregated time series, the PDO is related to a fraction of 0.37-0.44(depending on the variable) of SST-related variability.For the shortwave radiative fluxes σ lfpdo,ensmean /σ allSST,ens01 this is 0.24-0.25 (depending on the variable) on average for individual grid boxes in the Northern Hemisphere.For spatially aggregated time series in the Northern Hemisphere, the fraction of total variability attributable to the PDO is 0.33-0.35,for North America-0.39-0.42,for Europe-0.27-0.31,for India and Indochinese peninsula-0.22-0.27.The fractions of variability presented in this section are indicative for the variability of annual mean time series and become larger for longer time scales, at which higher frequency atmospheric noise plays a smaller role (not shown).

Quantification of Shortwave Radiative Flux Anomalies Associated With PDO Phases
We remove the temporal element of the PDO evolution and use the "perpetual" PDO set of simulations to quantify the increase in each of the radiative fluxes as a function of the PDO index value which is embroidered in the climatological SST field as anomalies in space (exact SST anomalies are shown in Figure A3).Each simulation is conducted with a constant in time PDO index value; the index values are chosen as different percentile of the distribution of all PDO index values obtained from the 500-year coupled simulation.The remaining oceans (outside the North Pacific) are kept at their climatological monthly means.By computing the difference between the mean time series of each experiment and the control run (neutral PDO), we obtain radiative flux anomalies per grid box which correspond to the different simulations (PDO indices).

Spatial Extent
The spatial distribution of the radiative flux and surface temperature anomalies are shown in Figure 5. Comparing the spatial patterns to the trends distribution obtained from CMIP6 (Figure 2) for F ↓ S , F ↓↑ T and F ↓↑ S , we observe similarities in the Northern Hemisphere, including the North Pacific region, North America, Europe and Southern Asia.A prominent difference observed in the Southern hemisphere is Australia.In the fully-coupled CMIP6 simulations, where the PDO occurs in combination with other modes of ocean variability, Australia evolves in phase with South Asia.By contrast, in our simulations where the PDO exists in isolation, this connection between Australia and South Asia is broken.The differences between our ICON-A simulations and CMIP6 will be further addressed in the discussion.
Comparing the simulations that differ in the strength and sign of the PDO SST anomaly (the different columns in Figure 5), we observe prominent and statistically significant differences in the simulations that correspond to higher percentiles (95, 99).In the p68, p68 experiments, where the spatial SST anomaly is not as strong, the patterns of the radiative flux response cover smaller areas that appear scattered, but the resulting anomaly is of the same sign as in the simulations with stronger anomalies.We assess the statistical significance per grid box by performing Student t-test (Student, 1908) between the time series in the control run and the one in the corresponding experiments.Our null hypothesis that the two samples come from the same population is rejected in 23%-29% (depending on the radiative flux variable) of the boxes in the Northern Hemisphere for the p95 simulation and in 40%-42% for the p99 simulation.These include most of the North Pacific Ocean, North America, South India, the Indochinese Peninsula and East Africa, more pronounced in the Ethiopian Highlands.The null hypothesis cannot be rejected at the 0.05 level for most of Europe (depending on the variable and experiment), indicating that the atmosphere-only variability can easily mask the PDO signal in a statistical test.
Based on the anomalies that correspond to p95 and p95 (second and fifth columns in Figure 5), the surface and TOA radiative flux anomalies reach up to ±6 Wm 2 for parts of North America, ∓3 Wm 2 for India and ±2 Wm 2 for parts of Europe.These numbers represent per grid box anomalies associated with the PDO range of [ 1.64; 1.64] in the 95-year mean where atmospheric noise is averaged out.
To draw a more general picture and exploit the connections and similarities between the radiative fluxes, we decompose the observed changes to changes in the optical properties of the atmospheric column (transmittance, reflection and absorption) and changes in surface reflection (albedo) in relation to the PDO anomaly (our different simulations).To do so, we closely follow the decomposition described in Stephens et al. (2015) and Loeb et al. (2019), who use the upwelling and downwelling shortwave fluxes at TOA and at the surface: Exactly following Stephens et al. (2015), where further details may be found, we define the effective surface albedo α, the system transmittance T and the system reflectance R as: Figure 5. All-sky shortwave flux, cloud radiative effect and surface temperature anomalies computed for each "perpetual PDO" experiment ( PC 1 (t) ∈ [ 2.58( p99), 1.64( p95), 0.44( p68), 0.44( p68), 1.64( p95), 2.58( p99)]) relative to the control run ( PC 1 (t) = 0).Hatched regions represent areas where the Student-t test is passed at the 0.05 significance level, that is, the null hypothesis of equal population means of the two samples (time series in experiment and control in individual grid boxes) cannot be rejected.
where t and r are the atmospheric transmissivity and reflection for a single beam, given the assumption that the atmosphere reflects and transmits equally in both upward and downward directions.The system reflectance R includes the reflected energy from the atmosphere rS plus the multiple scattering between the surface and atmosphere, it corresponds to Earth's planetary albedo (Stephens et al., 2015).
We can compute t and r from the T and R (and the radiative fluxes respectively) by: to better carve the energy partitioning, we express the atmospheric transmissivity t as: where a = 1 r t is the fraction of the beam that is absorbed by the atmosphere.This gives a more clear separation on whether the excess energy that does not reach the surface goes into heating the atmosphere.
The spatial patterns of t, r, a and α are shown in Figure 6.In addition, we explore several physical and microphysical variables: cloud cover, cloud liquid water and ice, and water vapor reveal similar patterns shown in Figure B1.According to our simulations, the changes in the radiative fluxes are due to both changes in cloud cover fraction as well as cloud water and ice content.It is evident that r closely resembles the patterns of changes in CRE (Figure 5), cloud cover, cloud liquid water and cloud ice content (Figure B1). a is similar to A atm , as is to be expected, but their pattern does not closely follow the column integrated water vapor content (Figure B1), as clouds are also essential for atmospheric absorption (see Section 3.1).The surface albedo α has its own unique pattern, which is enhanced when looking only at winter months (not shown), indicating that the changes are associated mainly with changes in snow cover.There is an increase in α in North America during a positive PDO, in line with increased r in the region.Over Eurasia, the change in α is not symmetric for the positive and negative PDO phases: the increase in α during a positive PDO seems to extend South to around 55°N, while during a negative phase, there is an increase in North Siberia (South to 60°N), and a decrease further South.

Linearity of the Response
To draw whether the response of the large-scale radiative fluxes to the SST anomaly related to the PDO is linear or not, we take the anomalies of the spatially averaged time series in the regions discussed in the previous sections-North America, Europe and India and Indochinese peninsula and plot the mean anomalies for each simulations as a function of the "perpetual" PDO index driving the simulation.The anomalies are computed by aggregating the time series in the region, taking the time-mean across the 95-year run and subtracting the control run aggregated temporal mean.Resulting dependencies are shown in Figure 7 with the anomaly for India and Indochinese peninsula given with a negative sign to match the direction of the other curves.In the range of PDO indices we work in, the response of the radiative fluxes is mostly linear, even though there is asymmetry around zero for some regions and variables.For the aggregated time series in North America and India and Indochinese peninsula, the anomaly is about ±4 Wm 2 for PDO index values of 2.58, 2.58 ( p99, p99 in the ICON-ESM-LR piControl distribution) and about ±1 Wm 2 for Europe.These values depend on the exact choice of the regions, which is only done for illustrative purposes and is arbitrary.
Results for Europe are indistinguishable from the control run, given the 95% confidence interval of the uncertainty of the mean, which is also evident from the statistical tests performed per grid box (Figure 5).However, the confidence intervals do not overlap for the highest and lowest PDO values, showing a distinct and statistically significant difference between the positive and negative values.The increase of both a and r from negative to positive PDO values, shows that they both contribute to the reduction F ↓↑ S with the contribution of r being much larger.This is also evident when comparing the net radiative anomalies at the surface and TOA with F ↓↑ T being always slightly smaller than F ↓↑ S , implying that the increased atmospheric absorption slightly offsets the reduction in F ↓↑ S for a positive PDO phase.Unlike the radiative fluxes, the surface temperature, T S , dependence is different in North America and Europe-the anomaly is negative in North America and small but positive in Europe, which is also evident in the patterns in Figure 5.

S
Lastly, we assess the relative contributions of r, a and α to changes in F ↓ S for the different PDO anomalies (simulations).We do so by expressing F ↓ S as: Since S in our case does not have any inter-annual variability and by also omitting covariance terms between δr, δa and δα, we can approximate the change in F ↓ S as: (1 r a)( rδα αδr) The linear equation above allows us to attribute the changes in downwelling solar radiation δF ↓ S to changes in the atmospheric and surface properties: δr, δa and δα.Similar expressions can also be obtained for the net fluxes at the surface and TOA but for this we limit ourselves only to downwelling surface radiation.
Taking the relative contributions as δaS/δF ↓ S , δrS/δF ↓ S and δαS/δF ↓ S , we find that for the individual regions they are approximately constant for the different PDO index values (i.e., relative contributions of the different terms do not change with the PDO value), especially the higher ones, where δF ↓ S is more distinct.The first thing we highlight is that the relative contributions depend on the region.For the ocean regions, where α is essentially constant, we note the following: in the Tropical extension of the PDO, almost all of the change in F ↓ S comes from changes in the reflectivity of the atmosphere δr; the same is true for the West part of the PDO region (cold region during a positive PDO phase) with only 2%-3% of the changes in F ↓ S being attributable to δa; in the East part of the PDO (warm region during a positive phase), 86% of the changes in F ↓ S are attributable to δr and around 14% to δa.Over India and Indochinese peninsula, 91% of the changes are attributable to δr and around 9% to δa.In the North American region, 93% of the changes can be attributed to δr, 13%-to δa, and 6% are offset by δα, which dampens the reduction in F ↓ S by increased snow cover during the positive PDO phase and the other way around for the negative phase.For Europe (excluding Scandinavia), the increase in r and a, and a reduction in α, all contribute to a reduction in F ↓ S by 80%, 17% and 3% respectively.All percentages given are average contributions of the terms for the p99, p95, p95 and p99 simulations.The residual contributions of the covariance terms are less than 0.1% for all regions and experiments investigated.

Discussion
The discussion section is structured as follows: in the first part we compare results from our "perpetual" PDO ICON-A simulations to the coupled ICON-ESM-LR and CMIP6 multi-model median to carve out differences that are due to coupling and differences between ICON and the median across CMIP6 models.In the second part, we compare the PDO index amplitudes within ICON and observations, namely ERSST V5 (Huang et al., 2017).are evident in both CMIP6 multi-model median (Figure 2) and ICON-ESM-LR (Figure C1), implying that they result from the coupling or a large-scale process related to the PDO that is not included in our experimental setup rather than differences between the ICON model and CMIP6.This suggests that the mechanisms responsible for Australian and South American radiative flux anomalies do not depend only on the PDO as a horseshoe phenomenon in the North Pacific, but on its basin-wide manifestation-IPO, which does not exist in that form in our PDO-only simulations.Furthermore, the radiative fluxes in the South Pacific appear to be affected in a similar (but weaker) way as in CMIP6, which in a coupled simulation would have an effect on the SSTs and the dynamics thereafter.This change in the shortwave radiative fluxes (including unforced dimming and brightening) might be one of the mechanisms through which Northern extra tropics communicate with Southern hemisphere and the PDO influences the SPDO and they combine into IPO.
The dimming above the Atlantic during PDO↑ that is evident in CMIP6 and ICON-ESM-LR (Figure C1) but not in our ICON-A (Figure 5) simulations implies that it either arises from the interaction between the atmosphere and ocean or there are other processes, independent from the North Pacific horseshoe pattern, that are related to the PDO.A mechanism that we propose is that the North Atlantic dimming in the coupled models might arise from the ocean transport (by surface currents like the Gulf stream) of cold waters from the East coast of North America into the ocean interior.The East coast is influenced by PDO-related dimming that is also evident in our atmosphere-only ICON-A simulations.This gives a hypothetical example of a communication mechanism between the North Pacific and North Atlantic oceans through surface shortwave fluxes.
The difference in the region of the Mediterranean Sea, Arabian peninsula and Sahara: dimming in CMIP6 multimodel median and no-specific response or weak brightening in our ICON-A simulations (Figure 5) are also evident between the CMIP6 multi-model median (Figure 2) and ICON-ESM-LR (Figure C1), which implies that they are due to differences in the specific model response of ICON.
We also highlight that the correlation coefficients between ICON-ESM-LR and the shortwave fluxes are higher as compared to the majority of CMIP6 models, but such high correlation coefficients are also evident in model families such as CESM2, CMCC, GFDL and MIROC (Table 1).This stronger relationship is also evident in the stronger trend magnitudes in ICON-ESM-LR (Figure C1) as compared to CMIP6 multi-model median (Figure 2).It is also related to the large inter-model spread in the trend magnitudes estimated in Section 3.1.Whether the differences among models are due to a weaker SST anomaly or a different atmospheric response requires further investigation.
We next compare the PDO that we derived from ICON-ESM-LR to the PDO time series computed with the same method upon the ERSST-v5 reconstructed SST data (Huang et al., 2017) that covers the period 1920-2015 on a 2°grid with a monthly resolution.We decrease the imprint of global warming on the observational SST field by subtracting the global annual mean SST anomaly from each grid box.Extracting the first eigenvector (that corresponds to the PDO horseshoe) from the deseasonalized monthly anomalies in the observations, yields the historical PDO time series.The variance fraction explained by this PC is 17%, in comparison for ICON-ESM-LR, the variance fraction is 25% (based on the ratio between the first eigenvalue divided by the sum of all eigenvalues).To account for differences in the variance fraction, the PC time series are usually normalized by their eigenvalues, which relates the magnitude of the spatial pattern EOF and the temporal variability PC.We then compare the 95th percentiles of the resulting distributions of PDO index values, that are 1.64 for ICON-ESM-LR and 1.58 for ERSST-v5, yielding a slightly stronger PDO amplitudes in ICON-ESM-LR as compared to ERSST-v5.The differences in variance fraction also imply that the observational SST field contains much more "noise" on top of the idealized PDO evolution that we analyze in this study.The pattern correlation for PDO between the ICON-ESM-LR and ERSST-v5 is 0.85 (computed on the ICON grid), typical for CMIP model ranges of 0.8-0.9(Newman et al., 2016).
One might pose the question on why we use SSTs from the coupled piControl run instead of directly taking them from observations.The reasons for this are the following: (a) the piControl simulation is not contaminated by global warming and climate change; (b) we can build better statistics and derive a cleaner eigenvector from 500 years of data as compared to the shorter observational period; (c) the atmospheric dynamics within the model is aligned with the SSTs, which might be of relevance because SST patterns are a combination of atmospheric and oceanic processes and having a numerical mismatch between the atmosphere and SSTs might compromise a study that focuses solely on internal variability.Given the comparison to observed SSTs above, we do not have indications of significant biases in the PDO pattern selected for our experiments.

Summary
The present study targets the role of one specific mode of ocean variability-the PDO, in the variability of the shortwave flux components within Earth's energy balance: F ↓ S , F ↓↑ T , F ↓↑ S , A atm , CRE S , CRE T .First, we look at them through the lens of the PDO in coupled simulations, that is, the PDO as it is diagnosed via the PDO index in CMIP6 piControl simulations, where the PDO is naturally co-occurring with other modes of variability, notably ENSO and the SPDO.Seen through this PDO lens, the all-sky components of F ↓ S , F ↓↑ T , F ↓↑ S , CRE S , CRE T show remarkably similar patterns associated with the PDO, highlighting the redistribution of clouds in sync with the PDO in the coupled system.On the other hand, changes of all-sky A atm and the clear-sky components in relationship to the PDO show different patterns that reflect the redistribution of shortwave absorbers in the atmosphere (mainly water vapor) and changes in effective surface albedo in response to the PDO.
The shortwave flux redistribution as a result of the PDO phase shifts described here has several prominent features which find support in observations.The first one is the increase in absorbed solar radiation off the coast of the western U.S. and the decrease above the continent.Above ocean, Loeb et al. (2018Loeb et al. ( , 2022) ) link the PDOassociated increase in absorbed solar radiation to decreased low cloud amounts in the Eastern Pacific, while above the continent, Augustine and Capotondi (2022) link the PDO to a decrease in surface solar radiation, which both agree with our model-based results.Furthermore, an observational record from Nandi-Fiji in the West Pacific reveals brightening during 1970s-1990s (PDO↑) and dimming afterward 1990s-2010s (PDO↓) (Correa et al., 2024)-opposite to the observed trends in Europe and the U.S. during the same periods (e.g., Augustine and Hodges (2021);Sanchez-Lorenzo et al. (2015)) and in agreement with our model-based results.
Adjusting the PDO lens to focus, by experiment design, exclusively on the North Pacific horseshoe and its yearto-year variability, we ask the question what is the fraction of total year-to-year variability that can be related to the PDO.To do so, we run ensemble atmosphere-only simulations with ICON-A (each with length 100 years) that allow us to average out atmosphere-only variability for an all SST ("allSST") case and a low frequency PDO case ("lfPDO").The "allSST" experiment contains the SSTs directly taken from the coupled piControl run of ICON-ESM-LR, while the "lfPDO" experiment contains only the low-frequency component of the PDO time series projected on the SST field.We show that SST-related variability (including ENSO) is related to almost all variability of shortwave fluxes in the tropics and about half of the variability in the extratropics.The PDO alone (constrained to the North Pacific) is related to almost none of the variability in the tropics and around 20%-40% of variability of the shortwave fluxes over the Northern Hemispheric continents.We also show that over land, SSTrelated variability is more evident in spatially averaged (aggregated) time series as compared to individual grid boxes.This is because the SSTs impact on large-scale dynamics is prone to be blurred by small-scale noise.
Keeping the idealized spatial form of the PDO lens but now idealizing also the temporal manifestation, we carry out a set of "perpetual" PDO simulations, each corresponding to a different PDO anomaly in space.When comparing the geographical patterns of shortwave flux anomalies we obtain this way to geographical patterns of trends within the CMIP6 data set, we find overall agreement in the Northern Hemisphere.We show that the mean anomalies on a grid box level for North America reach up to ±6 Wm 2 ; for India and Indochinese peninsula-up to ∓3 Wm 2 ; and for Europe-up to ±2 Wm 2 .Numbers represent anomalies associated with the PDO range of [ 1.64; 1.64] ( p95 and p95) taken from the temporal mean of each 100-year simulation, where atmospheric variability is averaged out.

Focusing on changes in F ↓
S , induced by the PDO, we find that the relative contributions of changes in fractional atmospheric reflection r, absorption a, and effective surface albedo α depend on the region.r, which is related to changes in both cloud cover and cloud water content has the largest contribution in all regions with α (primarily related to snow cover) slightly offsetting the effect for North America and a, which is related to both water vapor and clouds, playing a larger role for PDO-related F ↓ S anomalies in Europe and a smaller role for India and Indochinese peninsula.Magnitudes in F ↓↑ T are in general less than those in F ↓↑ S , which is explained by the increase of both r and a in those regions-increased atmospheric absorption slightly offsets the decrease in absorption at the surface.
The combined view through these different PDO lenses, ranging from a "realistic" PDO embedded in its natural context of other modes of variability in CMIP6 to a highly idealized, eternally phase-locked PDO within our ICON-A simulations, inspires some hypotheses and speculations.By assessing the similarities and differences in these patterns, we mark regions that show a different response in coupled and atmosphere-only simulations.The For a better illustration of how the SST field is being modified for each experiment, we include the anomaly time series (before adding the seasonality) for different grid boxes in Figure A2.The locations of the points are shown in Figure A1a and we choose them as: the green point represents a location with a strong positive correlation with the PC 1 (t) time series; the orange point-strong negative correlation; the blue point is within the North Pacific  A1 for 100 years.The green point as chosen with a strong positive correlation with PC 1 (ρ(PC 1 )), the orange point-with a strong negative, the blue point-weak correlation.The yellow point is strongly correlated to PC 1 but lies outside the analysis box and is affected by the distance factor.Black curve shows the original SST time series from the piControl run (also the control simulation), orange curve shows time series from the SST field after projecting the low-pass filtered PC 1 .Gray lines depict the climatological simulation (all values are zero), blue lines depict locked negative phase Pacific Decadal Oscillation (PDO) simulations and red lines-locked positive phase PDO simulations.Seasonality is not shown in the plots but included in the SST boundary conditions.box but is not correlated with the PC 1 (t) time series; and the yellow point is strongly correlated with the PC 1 (t) time series but lies outside the North Pacific box and is therefore again weakly affected by our modification (projected anomalies are smaller in magnitude as compared to the green and orange points).The spatial representation for each of the PC 1 (t) = const projections is given in Figure A3, where we also show the mean SST field and the climatology of sea ice that is kept constant across simulations, that is, does not exhibit any inter-annual variability.

Note.
The shortwave components are taken as spatially aggregated annual-mean time series for four different regions (drawn in Figure2).Bold font indicates absolute value of the correlation coefficient larger than or equal to 0.5.All correlation coefficients with absolute values above 0.1 are statistically significant at the 0.05 level.Correlation coefficients are calculated over the entire piControl simulations, which yields a different number of years for each model.scale trends of the various shortwave fluxes.Sections 3.2 and 3.3, by contrast, look at a PDO the pattern of which is highly idealized by design, comprising only the North Pacific region.Associated atmosphere-only simulations with the ICON-A model are used to examine the imprint of this idealized PDO on shortwave fluxes: how the yearto-year variability of the PDO mirrors in shortwave flux time series (3.2) and how shortwave fluxes differ between ever-lasting positive and negative PDO phases, respectively (3.3).The latter experiments offer an idealized view on strong PDO phase changes.

Figure 1 .
Figure 1.First eigenvector and principal component (PC) time series derived from empirical orthogonal function (EOF) analysis of the sea surface temperature (SST) field in ICON-ESM-LR's piControl run.The first EOF (eigenvector), shown in panel (a), is presented as the correlation with the monthly SST time series for each grid box.Blue box shows the region in which the EOF decomposition is done.In (b), the normalized value of the original first PC (Pacific Decadal Oscillation index) time series is displayed in blue and its low-pass filtered counterpart in orange.Shown are only the first 100 years.Low pass filtering is done by smoothing out frequencies of larger than 1/10 years with a sharp logarithmic function (described in Appendix A).Panel (c) displays the power spectral density estimated using Welch's method(Welch, 1967) for the entire 500year simulation.

Figure 3 .
Figure 3. Same as Figure 2 but for the clear-sky fluxes.

Figure 2 .
Figure 2. All-sky shortwave flux and cloud radiative effect trend patterns during 20-year periods with a strong negative-to-positive phase transition of the Pacific Decadal Oscillation (PDO) index (PDO↑) and periods with a strong positive-to-negative transition (PDO↓) obtained from Coupled Model Intercomparison Project-Phase 6 (CMIP6) piControl simulations.Each map is composed of radiative flux trends which correspond to the strongest 10% trends in the PDO index, averaged for each CMIP6 model.Multi-model median is computed among the maps obtained from 52 CMIP6 piControl simulations.Hatched regions represent areas where the differences between PDO↑ and PDO↓ are not statistically significant at the 0.05 level as assessed with the Mann-Whitney-U test.Black boxes on the last map are those, over which we compute aggregated time series.

Figure 4 .
Figure 4. Standard deviation of the annual mean time series per grid box of the 20-ensemble-member mean "lfPDO" experiment for different variables (first column), ratio between the standard deviation of the 20-ensemble-member mean in the "allSST" experiment and the standard deviation for the control run (second column) and ratio between the standard deviation of the 20-ensemble-member mean in the "lfPDO" experiment and the control run (third column).The control run is one ensemble member of the "allSST" experiment.

Figure 6 .
Figure 6.Same as Figure 5 but for anomalies in atmospheric transmittance, reflection, absorption and surface albedo.

Figure 7 .
Figure 7. Scatter plots for different variables as a function of the Pacific Decadal Oscillation index value for regions denoted in Figure 2. The regional means for India and Indochinese peninsula are multiplied by 1. Uncertainty ranges show the 95% confidence interval determined via t-statistics.Values on y-axis depend on the exact choice of region boundaries, which is arbitrary.

Figure A2 .
Figure A2.Time series of the first principal component and deseasonalized monthly sea surface temperature (SST) anomalies at different locations, depicted in FigureA1for 100 years.The green point as chosen with a strong positive correlation with PC 1 (ρ(PC 1 )), the orange point-with a strong negative, the blue point-weak correlation.The yellow point is strongly correlated to PC 1 but lies outside the analysis box and is affected by the distance factor.Black curve shows the original SST time series from the piControl run (also the control simulation), orange curve shows time series from the SST field after projecting the low-pass filtered PC 1 .Gray lines depict the climatological simulation (all values are zero), blue lines depict locked negative phase Pacific Decadal Oscillation (PDO) simulations and red lines-locked positive phase PDO simulations.Seasonality is not shown in the plots but included in the SST boundary conditions.

Figure A3 .
Figure A3.Sea surface temperature (SST) anomaly associated with different percentiles ( p) of the monthly Pacific Decadal Oscillation index distribution that we uses for the corresponding experiments.Middle plot shows the mean SST field on top of which we add the anomalies and the three shades of turquoise represent the maximum, mean and minimum extent of sea ice within 1 year.

Table 1
Pearson Correlation Coefficients Between the Annual Mean PDO Index and All-Sky Shortwave Components of the Energy Balance: Downwelling at Surface (F ↓