Surface and Tropospheric Response of North Atlantic Summer Climate from Paleoclimate Simulations of the Past Millennium

: We investigate the effects of solar forcing on the North Atlantic (NA) summer climate, in climate simulations with Earth System Models (ESMs), over the preindustrial past millennium (AD 850–1849). We use one simulation and a four-member ensemble performed with the MPI-ESM-P and CESM-LME models, respectively, forced only by low-scaling variations in Total Solar Irradiance (TSI). We apply linear methods (correlation and regression) and composite analysis to estimate the NA surface and tropospheric climatic responses to decadal solar variability. Linear methods in the CESM ensemble indicate a weak summer response in sea-level pressure (SLP) and 500-hPa geopotential height to TSI, with decreased values over Greenland and increased values over the NA subtropics. Composite analysis indicates that, during high-TSI periods, SLP decreases over eastern Canada and the geopotential height at 500-hPa increases over the subtropical NA. The possible summer response of SSTs is overlapped by model internal variability. Therefore, for low-scaling TSI changes, state-of-the-art ESMs disagree on the NA surface climatic effect of solar forcing indicated by proxy-based studies during the preindustrial millennium. The analysis of control simulations indicates that, in all climatic variables studied, spurious patterns of apparent solar response may arise from the analysis of single model simulations.


Introduction
Changes in the climate system may be a result of internal climate variability or external influences that may be either anthropogenic or natural (i.e., orbital, solar, or volcanic forcing). The human influence on global temperature since the start of the industrial era (after AD 1850) has very likely exceeded the impact of natural forcings [1]. However, during the pre-industrial era, externally driven climate change at decadal-to-multidecadal time scales is considered to be forced primarily by variation in solar output and by volcanic eruptions [2][3][4]. In contrast to volcanic forcing, solar forcing may cause long-term (decadal mean) continental and regional climate changes that are greater than unforced (internal) variability [5].
The Sun's climatic impact arises from changes in its total solar irradiance (TSI), spectral solar irradiance (SSI) and energetic particle precipitation (EPP) [6][7][8][9][10]. Recently, more efforts were placed on the representation of variations in SSI and EPP in climate model simulations, with new recommendations regarding the magnitude of SSI variations used to force state of the art global climate models. The SSI uncertainty and the possible impacts of the higher SSI variations are discussed extensively in Matthes et al. [11].
In the current study, we are interested in the climatic impact of TSI variations during the preindustrial last millennium (AD 850-1849), as TSI is used as proxy for identifying the effect of solar variability in paleoclimate proxy data. The effect of TSI changes on the climate of the past millennium was identified by early works in the context of the historical Maunder Minimum (AD 1675-1715) [12] and confirmed by more recent studies investigating climate changes during periods of solar maxima and minima [13][14][15][16]. In the early 2000s it was generally accepted that a long-term change in solar activity between the mean solar output during the Maunder Minimum (decline in solar activity) and the present-day climate, a period with presumably higher solar output, was about 0.3% [17]. However, the amplitude of these TSI changes was subsequently revised and reduced to a value of approximately 0.1% [18,19]. This new scaling was applied to the models participating in the 5th phase of the Coupled Model Intercomparison Project-CMIP5 [20,21]. Even though this change in the TSI scaling has implications for our understanding of the effects of changes in solar output on climate, a quantitative assessment of this effect has not yet been made.
We use simulations forced with the latest estimations of the solar amplitude and investigate how these TSI variations are reflected in the North Atlantic (NA) climate of model simulations in the preindustrial last millennium. We focus on the NA, where oceanatmosphere interactions crucially affect the climate of North America and Europe [22][23][24], and on the summer season, as it is more relevant for comparison to studies using proxy data. Biological proxies, such as trees and bivalve shells, tend to reflect the growing season more strongly. The variability in their growth increment widths is related to environmental conditions and their growth is biased towards summer [25,26].
Investigating the Sun-climate relationship prior to industrial times is important for the comparison between climate models and proxy records and for disentangling anthropogenic climate change from natural variations [27]. Modeling and proxy-based studies have so far disagreed on the importance of TSI forcing on driving variations in global annual mean surface temperatures [28][29][30], but they agree on larger regional effects of TSI forcing [31][32][33]. Nevertheless, these effects might be too large to be attributed to direct solar heating, and positive feedbacks have been invoked [34], particularly in regions that usually are cloud covered [35][36][37]. Regarding the NA region, some authors have claimed the detection of influence of solar forcing in marine proxies over the past millennium for Sea Surface Temperatures (SSTs) in centennial and multi-decadal time scales [33,38,39].
Jiang et al. [38] used regression analysis and paleoclimate proxy data, and identified a positive correlation between SSTs and solar forcing for the region of the North Icelandic shelf. Further, these authors used composite analysis and found that surface temperatures north of Iceland increase for stronger solar forcing. In order to investigate whether such signals occur in model simulations, we make use of linear regression and composite analysis, and investigate the spatial climatic changes that are solely induced by changes in solar activity, during the summer seasons of AD 850-1849, over the North Atlantic. For this goal, we use solar-only forced CMIP5-type simulations available for the Max Planck Institute Earth System Model (MPI-ESM) and the Community Earth System Model (CESM). These simulations employ reduced long-term solar amplitudes, smaller than that used in simulations carried out one or two decades ago, and implement solar activity changes for different wavelength bands.
The use of only-solar forced simulations is important for having a clearer solar forcing signal, as solar forcing effects might be concealed by the larger effect of volcanic eruptions, for example, in periods when TSI minima coincided with volcanic eruptions during the last millennium. For this reason, the use of linear regression might not be appropriate for detecting the signal of solar forcing [40]. Moreover, the CESM output includes an ensemble of solar-only forced experiments. That is important in order to quantify the magnitude of model internal variability and the magnitude of the TSI effects on the modeled climate.
Other studies also used simulations with current estimations of long-term solar amplitudes and investigated the sun-climate relationship, but they investigated the winter season in the recent historical period and used fully forced experiments [41][42][43][44]. Our analysis focuses on summer season and the preindustrial period, therefore being more relevant to model-proxy comparisons. Moreover, the use of only-solar forced experiments during the preindustrial era automatically excludes possible influences from volcanic eruptions and anthropogenic forcing. Therefore, the investigation of the pure impact of changes in solar activity on the simulated climate is more robust. Another aspect of the paper is that our analysis includes different techniques for signal identification. Therefore, our results might help to disentangle potential solar-only induced changes from those of the fully forced runs, focusing spatially on the effect of model internal variability.  [45]. This configuration is identical to MPI-ESM-LR with two exceptions: the dynamical vegetation module is not active, so that land use changes can be prescribed, and the orbital parameters are prescribed by a table providing annual values for eccentricity, obliquity, and perihelion [46].

Materials and
The atmosphere has a T63/1.9 • horizontal resolution and 47 hybrid sigma pressure levels (ECHAM; [47]). The ocean grid is curvilinear bipolar with an effective 1.5 • resolution (near the equator) and 40 z-levels (MPIOM; [48]). The atmospheric grid extends to 0.01 hPa in the vertical domain, where the center of the topmost level in the upper mesosphere is set. This vertical extension of the atmospheric grid includes for the first time a better resolved stratosphere in CMIP simulations for MPI-M [49]. The solar spectrum is split into 14 spectral bands and temporal variations of solar irradiance depend on the wavelength. Here the relative variations among spectral bands are the result of a physically reasonable estimation, since the information provided by the ice-core records cannot differentiate between spectral bands [50]. The monthly average ozone concentrations for the period AD 850-1849 are also prescribed and were calculated using the AD 1850-1860 monthly climatology of ozone concentrations from the AC&C/SPARC Ozone Database as a basis. The ozone dependency on solar irradiance is represented through regression coefficients between historical ozone concentrations and the annual 180.5 nm solar flux [45].
We use a single simulation of the MPI-ESM-P model, the solar-only forced simulation MPI-R1. For solar radiation, the Vieira et al. [18] total solar irradiance reconstruction over the Holocene is employed with a TSI scaling of 0.1% from the 17th century Maunder Minimum to present [51]. Well mixed-greenhouse gases (CO 2 , CH 4 , N 2 O) and orbital parameters are annually varying [46].

CESM-LME
The second ESM used is the Community Earth System Model-Last Millennium Ensemble (CESM-LME; [52]). The CESM-LME employs the CESM-Community Atmosphere Model version 5 (CESM1-CAM5; [53]) with a difference in the resolution of its individual model components. The ocean component of CESM-LME has a 1º horizontal resolution, while the atmospheric component at 2 • horizontal resolution extends from the surface to the middle stratosphere (~40 km) [54].
Radiative fluxes and heating rates in CAM5 are calculated using the Rapid Radiative Transfer Method for GCMs (RRTMG) [55], which separates the short-wave spectrum into 14 bands extending from 0.2 µm to 12.2593 µm. The TSI value is scaled in each spectral band through a specification of a time-varying solar spectral irradiance [56]. To provide fluxes at the top of the atmosphere, RRTMG uses an additional layer above the CAM5 model top in both the longwave and shortwave. This extra layer is specified by replicating the composition of the highest CAM5 layer into a layer that extends from the top of the model to 10−4 hPa. The radiation parameterization requires monthly mean ozone volume mixing ratios to be specified as a function of the latitude grid, vertical pressure levels, and time. The ozone path lengths are evaluated from the mixing-ratio data derived by Chervin [57] from analysis of Dütsch [58]. The CESM-LME solar-only experiments are forced with the same annually varying well mixed-greenhouse gases and orbital parameters as the MPI-ESM-P simulation [46], following the PMIP3 protocol. The CESM-LME uses the [18] total solar irradiance reconstruction over the Holocene to produce a four member solar-only forced ensemble over the period AD 850-2005. The only difference among its four ensemble members is the application of small random differences (order of 10-14 • C) in the air temperature field at the start of each ensemble member. Both the CESM-LME and the MPI-ESM-P simulations are analyzed for the preindustrial period AD 850-1849. For better comparison, the original model output was post-processed and re-gridded to a regular grid (1 • × 1 • degree horizontal resolution). In Table 1, details on the models and experiments used in this study are given.  [18], Orbital forcing parameters and well mixed-greenhouse gases [46] 4/E1, E2, E3, E4 Solar forcing [18], Orbital forcing parameters and well mixed-greenhouse gases [46]

Methods
To identify the effects of TSI forcing on the NA climate, we used three methods. The first method, Pearson Correlation, and the second method, Linear Regression, are linear parametric techniques that identify the linear dependence between variables. Both of these methods take into account the full preindustrial period, AD 850-1849. This analysis is performed for decadal low-pass filtered series, motivated by the multi-decadal surface response found in proxy data [38,39]. The decadal filtered series were obtained by filtering the climatic parameters and the TSI series in the spectral domain (cut off frequency = 1/11) using a low-pass filtering technique [59]. The low-pass filter applied isolates the signals with a frequency lower than the selected cutoff frequency. Therefore, variability is suppressed for time scales shorter than 11 years. Moreover, for the correlation analysis we also consider possible temporal lags in the solar-climate relationship.
The third method used is Composite Analysis, focusing on the effect of prominent TSI periods with high and low TSI values. Periods with high TSI are defined as those with positive TSI anomalies exceeding the +2.2 standard deviation, whereas periods with low TSI are defined as those with negative TSI anomalies lower than the −1.4 standard deviation ( Figure 1a). These asymmetric requirements were set in order to have approximately the same amount of selected cases with positive (25 cases) and negative (20 cases) extreme TSI values, as the TSI distribution resembles a positively skewed non-Gaussian distribution ( Figure 1b). It must be noted that such analysis does not account for lagged responses in the climate due to TSI changes.
We focus on the spatial climatic response pattern of the North Atlantic Ocean due to TSI changes during summer (June-August), including changes in atmospheric circulation. The geographical domain is between 80 • W-30 • E and 21 • N-75 • N. For all methods used and variables, the model data have been detrended prior to the analysis. Moreover, the TSI data have been detrended prior to the analysis. The reason is to avoid the identification of spurious relationships caused by possibly existing, but physically unrelated, trends in the data [60]. We focus on the spatial climatic response pattern of the North Atlantic Ocean due to TSI changes during summer (June-August), including changes in atmospheric circulation. The geographical domain is between 80° W-30° E and 21° N-75° N. For all methods used and variables, the model data have been detrended prior to the analysis. Moreover, the TSI data have been detrended prior to the analysis. The reason is to avoid the identification of spurious relationships caused by possibly existing, but physically unrelated, trends in the data [60].
The statistical significance of the results obtained by these methods is tested with two approaches. One approach is the standard two-sided Student's t-test. This test is applied in the linear methods to estimate the significance of the regression and correlation parameters. To account for the high degree of autocorrelation in the low-pass filtered series, we use a random-phase bootstrap method. This method preserves the spectrum and autocorrelation of the original time series, but the surrogate time series have a random time evolution [61,62]. A t-test is also applied in the results obtained by Composite Analysis to test the mean difference between the high-solar and low-solar samples. Because of the small sample size, the confidence intervals are estimated using a bootstrapping technique with replacement. We select two random subsets from the original time series with lengths equal to the original composite subsamples. A distribution of the differences is then constructed by repeating this procedure 1000 times.
Our null hypothesis H0 is that there is no effect of solar forcing on climate, i.e., that the regression and correlation parameters are zero and that there is no mean difference between the populations selected from years with high and low solar activity. If the null hypothesis is rejected at the significance alpha level α = 0.05, then we can reject the null hypothesis, and interpret that there is an effect of solar output on climate, with a 5% risk. However, statistical significance does not necessarily imply physical relevance. Especially for the linear methods, we use a large pool of samples (n = 1000) and the sample size is important for the estimation of the nominal level of statistical significance. As consequence, very large sample sizes will ultimately lead to the phenomenon that even very small physical quantities reach the statistical significance threshold for the alpha level under consideration.
The second approach we use to test the statistical significance of our results, regarding the possibly emerging signals of solar impact on climate, is the analysis of Preindustrial Control (PiControl) simulations conducted with the same ESMs. These control runs are conducted with the same models, but with constant external forcing (with the exception of a constant annual cycle of solar insolation). Any signal of solar impact on climate The statistical significance of the results obtained by these methods is tested with two approaches. One approach is the standard two-sided Student's t-test. This test is applied in the linear methods to estimate the significance of the regression and correlation parameters. To account for the high degree of autocorrelation in the low-pass filtered series, we use a random-phase bootstrap method. This method preserves the spectrum and autocorrelation of the original time series, but the surrogate time series have a random time evolution [61,62]. A t-test is also applied in the results obtained by Composite Analysis to test the mean difference between the high-solar and low-solar samples. Because of the small sample size, the confidence intervals are estimated using a bootstrapping technique with replacement. We select two random subsets from the original time series with lengths equal to the original composite subsamples. A distribution of the differences is then constructed by repeating this procedure 1000 times.
Our null hypothesis H0 is that there is no effect of solar forcing on climate, i.e., that the regression and correlation parameters are zero and that there is no mean difference between the populations selected from years with high and low solar activity. If the null hypothesis is rejected at the significance alpha level α = 0.05, then we can reject the null hypothesis, and interpret that there is an effect of solar output on climate, with a 5% risk. However, statistical significance does not necessarily imply physical relevance. Especially for the linear methods, we use a large pool of samples (n = 1000) and the sample size is important for the estimation of the nominal level of statistical significance. As consequence, very large sample sizes will ultimately lead to the phenomenon that even very small physical quantities reach the statistical significance threshold for the alpha level under consideration.
The second approach we use to test the statistical significance of our results, regarding the possibly emerging signals of solar impact on climate, is the analysis of Preindustrial Control (PiControl) simulations conducted with the same ESMs. These control runs are conducted with the same models, but with constant external forcing (with the exception of a constant annual cycle of solar insolation). Any signal of solar impact on climate emerging in the control simulations, using the very same methods as for forced simulations, is due to random sampling and to the autocorrelation structure of the climate variables. Its amplitude yields a benchmark that needs to be surpassed by the solar-climate signals that are statistically identified in the forced simulations.

Results
Using the methodologies described before, we investigate the surface and tropospheric response to solar forcing over the North Atlantic during summer. Previous studies, mostly investigating the winter season, hypothesized that the simulated changes in solar forcing are integrated into the lower troposphere and therefore might affect both directly and indirectly the surface climate. The analysis is performed for decadal time scales using correlation and regression analysis. As suggested by [44], the use of the composite method avoids data filtering. Therefore, the results of the composite method do not take into account decadal time scales.

Linear Methods
The TSI time series is correlated and regressed to the oceanic and atmospheric variables for the period AD 850-1849. The Pearson correlation coefficients of TSI with sea surface temperature (SST), sea level pressure (SLP) and geopotential height at 500 hPa (geo500) are shown in Figure 2a-f, respectively. For the regression analysis, the TSI anomalies are standardized and subsequently regressed over the North Atlantic (NA) atmospheric and oceanic variables ( Figure 3). Therefore, the regression results represent changes of the climatic variables that correspond to changes of one standard deviation (SD) of TSI. The results are given in the first row for the MPI-ESM single simulation and in the second row for the CESM ensemble mean.
On decadal time scales, the TSI signal is generally weakly correlated (r < ±0.2) with SST and SLP in both models. Therefore, during the summer season, only a small fraction of the total variance over the surface and lower troposphere can be statistically linked to changes in TSI. Higher and statistically significant correlation coefficient values are shown over the subtropical regions (r < ±0.4) for the geo500 in both models, indicating a higher sensitivity to changes in TSI in the mid-troposphere.
The TSI regression patterns show similarities to the patterns shown by the correlation maps for all variables and both models. Although statistically significant at times, the response is generally weak over most of the North Atlantic, and its spatial pattern varies across models. This highlights the advantage of investigating two profoundly different models to rule out physically and externally forced robust patterns versus model-dependent ones. However, the comparison between the single MPI-ESM realization and the ensemble mean of CESM is complicated, because averaging within the CESM ensemble filters out changes in internal variability and hence should result in a higher signal-to-noise ratio. However, the signal-to-noise ratio (SNR) in the CESM simulations is low for all climatic variables studied (SNR < 0.4) and indicates that internal variability could mask the climatic signal of TSI forcing (for the definition of SNR see Supplementary Materials, Figure S1. A similar conclusion can be drawn by the correlation patterns of TSI to SST, SLP and geo500 for the individual CESM ensemble members (Supplementary Materials, Figure S2).      The spatial climatic response to TSI forcing is generally weak in the CESM ensemble and in the MPI-R1 simulation. Therefore, in order to test whether this statistical relationship between TSI and the variables under consideration could arise by chance, we performed a correlation analysis using the output from control simulations and the TSI series from the forced simulations (Supplementary Materials, Figure S3). The rationale here is that, due to natural decadal variability, spurious correlations to solar forcing can emerge between TSI and the climatic variables because of similar spectral characteristics included in the respective time series. As the control simulations are forced with constant solar forcing, potentially emerging patterns showing statistical significance indicate that the autocorrelation structure of the climatic variables may give rise to spurious correlation patterns. The control simulations of both models show nominally statistically significant correlation coefficient values of the same magnitude as their respective forced simulations. Focusing on the MPI-ESM model, we see that spurious patterns of apparent solar response may arise from the analysis of single simulations when using these statistical methods. The TSI-SST CESM control indicates a correlation pattern similar to the forced one, suggesting that there is no robust TSI response for the summer SSTs of the NA region. However, the correlation patterns emerging for the CESM SLP and geo500 are different between the control and forced simulations, with a large number of nominally statistically significant areas, indicating a robust forced response.

Lagged Climatic Response
In order to investigate whether there is a summer lagged climatic response to solar forcing, we performed a correlation analysis between the climatic variables and the TSI time series at different lead times. In Figure 4, we show the results of the correlation for the SSTs (column 1), SLP (column 2) and geo500 (column 3) for the realization MPI-R1. In the first row, the correlation with no lead-time is provided as reference. As the row number increases, the TSI lead time increases by one year. The respective results are given for the CESM ensemble mean in Figure 5.    In Figure 4, we see that within the MPI-R1 realization there are weak changes at all lags, with no intensification of SST, SLP, and geo500 correlation patterns resulting from the lag-0 analysis or emergence of new robust patterns. In the CESM ensemble mean (Figure 5), as the lead-time increases, there are more nominally statistically significant regions over the mid-latitudinal NA and Nordic Seas indicating a possible SST response to solar forcing. However, as previously discussed, there is a similar pattern emerging in the control simulation, indicating that the SST signal in the forced simulation is not robust. A slightly stronger signal over the subtropics is identified for the geo500 lagging two years, but no lagged response is found in SLP.
TSI forcing is expected to promote a latitudinal heating response at the stratosphere due to the ultraviolet absorption by the ozone layer. This should affect the zonal (U) wind component, which would communicate the signal to the lower atmosphere and to the surface. Another way that the zonal wind can be influenced by solar forcing is due to latent heat release and consequently though the upper troposphere. Therefore, we investigated the response of the zonal wind at the pressure level of 200 hPa (U200). Comparing the correlation maps TSI-wind obtained from the forced runs (Supplementary Materials, Figure  S4) to the correlation maps obtained from the control runs (Supplementary Materials, Figure  S5), we see that there are statistically significant changes in the U200 for both the MPI-R1 and the CESM ensemble. The strength of the U200 climatic response identified over the subtropics for the MPI-R1 weakens as the lag increases following the weakening in the elevation of the subtropical geo500 ( Figure 4). The CESM forced U200 response increases with TSI over the south of Greenland, agreeing with the negative correlation between TSI and SLP or geo500 over the same region ( Figure 5). It suggests that an increase in TSI reinforces the low-pressure system to the west, strengthening the jet above. As the lag increases, the reinforcement of the low-pressure system decreases, agreeing with the decreased values of the jet above.

Composite Analysis
The composite response of each atmospheric and oceanic variable analyzed in this study is estimated by subtracting the mean of the atmospheric and oceanic variables in the low TSI period from the respective atmospheric and oceanic variables within high TSI periods. In contrast to the low TSI sub-period, the high TSI sub-period is comprised mostly by single years and it does not contain periods longer than three years ( Table 2). As seen in Figure 1a, most of the years during high TSI are sampled from the first half of the study period (25 cases), while the low TSI years are sampled from the second half of the study period (20 cases). However, as all variables have been detrended, the results do not capture differences in the climatic variables' trends between the first and the second half of the preindustrial millennium. The composite difference maps are given in Figure 6, for the SSTs (a,d), SLP (b,e) and geo500 (c,f). As a benchmark, we also calculate the composite different maps derived from the control runs, by sampling the nominally same years as in the forced simulations (Supplementary Materials, Figure S6). The forced composite map of the CESM SSTs indicates that no signal can be attributed to TSI by the composite method during the summer season. The MPI-R1 shows an SST cooling off the European coast that is absent in the control simulation, but as this result comes from a single simulation it cannot be considered a robust forced response.

Common TSI Signals in the CESM Ensemble
As seen so far, the low TSI signal could not be robustly identified in the modele summer surface climate of the NA basin. A possible reason is that the SST-solar signal i concealed by the larger model internal variability. Even in the case of the four-membe CESM ensemble mean it might be that a larger number of the ensemble members i needed to clearly distinguish the solar signal from internal variability [63].
We have also explored whether the solar signal can be identified by analyzing th correlations between the individual simulations of the CESM ensemble. The rationale i that the correlations between climatic variables of two individual simulations can onl arise by common external forcing. If these correlation patterns are stable among all poss ble pairs in the ensemble, the signal will likely have a physical origin. We show in th Supplementary Materials, Figure S7 the correlation between all possible ensemble mem ber pairs of the CESM model, for the decadal filtered variables SST, SLP, and geo500. Ou null-hypothesis is that the CESM simulations do not show any response to changes i solar activity. The null-hypothesis is tested through point-wise correlations between th individual CESM ensemble members and can eventually be rejected for those region showing correlations exceeding the 5% significance level.
Given a resulting large number of statistically significant regions identified in th pairwise correlation maps between the same variables of the individual ensemble mem bers, one could conclude on a commonly forced signal exerted by changes in TSI as this i the only external information shared by all simulations. In order to identify which region show a common forced response, we then take into account from all the individual paire correlations only the regions of statistical significance. We find that with this approach we cannot identify regions with a commonly forced signal for any of the variables. Inter nal variability dominates the TSI signal for such small forcing and therefore the ensembl mean approach is critical for deciphering the forced response.

Sea Surface Temperature TSI Response
Two mechanisms act together in order to amplify solar forcing at the surface [64 These two mechanisms are a stratospheric response of ozone to fluctuations of shortwav solar forcing (top-down) [65,66], and the coupled ocean-atmosphere surface respons (bottom-up) [67][68][69]. As ozone is a radiative active gas susceptible to shortwave radiation The CESM ensemble mean indicates a weak tropospheric forced response to higher TSI values, which exhibits a different pattern from the one shown by the control simulation. This response is expressed with decreased SLP values on the southwest of Greenland and elevated geo500 values over some subtropical regions.
The forced patterns resulting from the composite method and from linear methods are expected to be different. The linear methods relate the evolution of the climatic variables to the TSI series, while the composite method depicts the mean difference of the climatic variables within extremely high and low TSI years. However, a similar forced response of the tropospheric variables was identified in the CESM ensemble mean by both methods, agreeing on a weak response with decreased SLP south of Greenland and increased geo500 over the subtropics. Moreover, all methods suggest that no robust forced response can be identified in the summer SSTs of the CESM ensemble mean in the NA basin.

Common TSI Signals in the CESM Ensemble
As seen so far, the low TSI signal could not be robustly identified in the modeled summer surface climate of the NA basin. A possible reason is that the SST-solar signal is concealed by the larger model internal variability. Even in the case of the four-member CESM ensemble mean it might be that a larger number of the ensemble members is needed to clearly distinguish the solar signal from internal variability [63].
We have also explored whether the solar signal can be identified by analyzing the correlations between the individual simulations of the CESM ensemble. The rationale is that the correlations between climatic variables of two individual simulations can only arise by common external forcing. If these correlation patterns are stable among all possible pairs in the ensemble, the signal will likely have a physical origin. We show in the Supplementary Materials, Figure S7 the correlation between all possible ensemble member pairs of the CESM model, for the decadal filtered variables SST, SLP, and geo500. Our null-hypothesis is that the CESM simulations do not show any response to changes in solar activity. The null-hypothesis is tested through point-wise correlations between the individual CESM ensemble members and can eventually be rejected for those regions showing correlations exceeding the 5% significance level.
Given a resulting large number of statistically significant regions identified in the pairwise correlation maps between the same variables of the individual ensemble members, one could conclude on a commonly forced signal exerted by changes in TSI as this is the only external information shared by all simulations. In order to identify which regions show a common forced response, we then take into account from all the individual paired correlations only the regions of statistical significance. We find that with this approach, we cannot identify regions with a commonly forced signal for any of the variables. Internal variability dominates the TSI signal for such small forcing and therefore the ensemble mean approach is critical for deciphering the forced response.

Sea Surface Temperature TSI Response
Two mechanisms act together in order to amplify solar forcing at the surface [64]. These two mechanisms are a stratospheric response of ozone to fluctuations of shortwave solar forcing (top-down) [65,66], and the coupled ocean-atmosphere surface response (bottom-up) [67][68][69]. As ozone is a radiative active gas susceptible to shortwave radiation, variations of the incoming solar UV radiation (120-350 nm) lead to changes in stratospheric ozone concentration and to heating that amplify the effect of the initial UV radiation changes in the Earth's atmosphere [65]. In contrast, the visible and infrared bands, which are much more weakly absorbed in the atmosphere, heat directly the Earth's surface and in turn the lower atmosphere.
For these reasons, both ESMs used in our study are expected to reproduce the North Atlantic surface response through direct solar heating, especially in the summer. The MPI-ESM-P is a high-top model, resolving the atmosphere as high as an altitude of 0.01 hPa. Even though it has a prescribed and not-interactive calculation of ozone, the ozone concentration scales with the strength of the TSI variations. It is therefore, expected to represent the top-down mechanism as well, which would act to enhance the surface response to TSI. The CESM-LME is a low-top model with no prognostic ozone chemistry. Therefore, it does not include the "top down" effect of solar variability, so an amplification of the surface response due to this mechanism can a-priori be excluded. Nevertheless, as shown by Mitchell et al. [43], the NA surface response to solar forcing is reproduced in both highand low-top models, but is more prevalent in the former. Moreover, both of the model configurations used in our study were shown by the study of Misios et al. [42] to be able to reproduce global responses to the 11-year solar cycle. Finally, a study using the CESM-LME solar-only forced ensemble found a response of the northern Asian winter monsoon to TSI forcing [70].
The results of other modelling studies investigating the TSI response of surface variables during the preindustrial period have indicated that the magnitude of surface response to TSI forcing is small but robust. One reason for the discrepancy of our results to other studies might be the considerably smaller TSI scaling that our ESMs use and the fact that we investigate the summer season. Regarding the NA basin, Ammann et al. [71] used a TSI scaling corresponding to a 0.25% change between the late Maunder Minimum and present. These authors found a statistically significant annual surface temperature response over a small region off the south coast of Greenland equal to 0.15 • C per W·m −2 . The response was calculated during the period AD 850-1849 using linear regression on yearly mean surface temperature. Swingedouw et al. [29] regressed 13 year low-pass grid-point annual SSTs on the TSI time series and found a general NA SST warming during AD 1001-1860, in simulations also forced with a 0.24% TSI scaling. This warming corresponds to approximately 1.75 • C per W·m −2 and is considerably larger than the one calculated in the current study. We used the CESM ensemble mean and decadal filtered SSTs and found a statistically significant, albeit not robust, SST warming of +0.1 K per one TSI SD or +0.3 K per W·m −2 .
Finally, a more recent study using the weaker 0.1% TSI scaling factor investigated which forcings were responsible for the near-surface air temperature (SAT) variability of AD 850-1850 period on sub-decadal time scales in all-forcing CMIP5 simulations [41]. These authors identified a relationship between TSI and SAT over the tropics and subtropics but no seasonal or annual relationship was found over most of the NA basin. Even though these results refer to SAT, they support our findings regarding the lack of a robust summer SST response over the NA basin.
Studies that explore the response of climate to TSI forcing during the industrial period (after AD 1850) have focused on winter circulation or on the annual effects of TSI forcing. These studies identify 2-4 year lagged surface temperature solar response signals in observations and in fully forced model simulations with low TSI scaling [44,60,[72][73][74]. Apart from the different period and the different season used by our study for the investigation of the solar signal, the dissimilarity of our study's results to the aforementioned modeling studies might come from the fact that these studies use single model simulations. Misios et al. [42] used an ensemble of fully forced CMIP5 simulations, and lead/lag multiple linear regression analysis, to investigate the observed and simulated surface response to solar forcing. Even though these authors identified observed solar signals on the surface and in the troposphere over several regions, there was no statistically significant observed annual signal identified in the NA basin. Moreover, over the NA region there was a high disagreement between the CMIP5 models for the annual, winter, and summer surface temperature solar response.
Not being able to identify the weak solar signal in the summer NA modeled climate during the last millennium might have implications on our ability to identify weak solar signals in proxy data collected in the NA basin. In this context, several studies based on proxy data report a surface climatic response to solar variability. Moffa-Sánchez et al. [39] used a marine sediment core from the basin of Iceland to reconstruct the temperature of the North Atlantic Current (NAC). These authors found a multi-decadal TSI response of the NAC temperature with correlation coefficient values of around +0.5. Moreover, the study of Sejrup et al. [33] demonstrates co-variability between a SST proxy and proxies for solar irradiance for the last millennium. Contradicting these studies, Turner et al. [75] examined peatland archives, which are widely used for the investigation of the sun-climate relationship [76][77][78]. These authors found that periodicities similar to the ones reported by the proxy archives as solar-type signals can be the product of random variations alone. This underpins the fact that more critical approaches might be required for the interpretation of proxy data.
Palaeoclimate assessment has demonstrated that changes in large-scale features of climate show consistent responses to external forcing, and that these responses are reproduced by paleoclimate models [79]. However, both numerical and interpretative uncertainties perplex the comparison between climate reconstructions and climate models at a regional scale [80,81]. This comparison is also hampered by the different background climate in the reconstructions and simulations, for example due to the different initial conditions at the start of the simulations compared to real world. Nonetheless, it can be hypothesized that sufficiently large changes in external forcings will have an imprint on both the proxy-based and the simulated climate, emerging from the internal background climate variations.
Even though we cannot directly compare the model world to the real world, it might be that internal variability is, in both cases, large enough to mask the response to solar activity changes. Still, an alternative explanation for the differences between our study and studies using proxy data relates to the small long-term scaling of the solar forcing being decisive for the very weak response in the modelled solar-climate relationship. Especially when investigating single model simulations, the model's own internal variability or model sensitivity is much higher for weak TSI forcing than for strong forcing [82]. Moreover, amplifying mechanisms of solar signals in the upper atmosphere, for instance related to the presence of interactive ozone might also be a factor that needs to be taken account for explaining the weak modelled response.

Tropospheric TSI Response and Dynamical Considerations
Studies investigating winter atmospheric circulation using ESMs [39], proxy data [83,84] and observational data [85,86] indicate the existence of blocking highs, as a response to low TSI forcing. Moreover, higher blocking persistence was found during low solar activity by a model-based study [87], while another modeling study supports that the winter SLP response to solar minima resembles a North Atlantic Oscillation-like pattern (NAO) [88]. Our results drawn by the CESM ensemble mean suggest that for the summer season the modeled TSI response in SLP and geo500 does not resemble a prominent pattern of tropospheric variability. Internal model variability is large, as depicted by the composite patterns of the individual CESM ensemble members (Supplementary Materials, Figure S8), with one member resembling a negative NAO (E2), and another a blocking high (E4) (for summer weather patterns classification see Figure 6.1 in [89]).
Except from the difference in season, the aforementioned modelling studies differ from our study, as they are not based on an ensemble mean and they use a combination of climate forcings. It is known that solar forcing and volcanic forcing appear correlated in time [90], and that there is a confounding effect of volcanic forcing when using simulations driven by all natural forcings. Dynamically, volcanic eruptions are known to produce a boreal-winter warming in some areas of the NH extra tropics due to the intensification of the winter Arctic Oscillation [5]. However, this winter warming depends on the strength of the volcanic eruption and also the ability of the models in realistically representing stratospheric processes in the context of stratospheric warming events due to volcanic aerosols [91]. The ocean warming in winter caused by volcanoes may partly persist into the summer season and thus confound the effect of changes in TSI in those regions. Moreover, uncertainties exist also in the studies using observational data to investigate the TSI-NAO connections, as the solar signal is investigated only in single SLP reconstructions [60].
Our analysis based on the CESM ensemble mean indicates a summer solar forced response of tropospheric variables such as SLP, geo500 and U200, but no SST response to TSI forcing is found. We showed that changes in the mean atmospheric flow are connected to changes in the strength of the tropospheric jets, and represented by the CESM-LME, albeit being a low-top model. Changes in the position of the subtropical tropospheric jets have been associated to anomalous heating of the upper tropical troposphere due to latent heat release in solar maximum [92]. Such changes in the jets position do not necessarily require a warming in the stratosphere and a consequent top-down response [93]. Moreover, our results for a weak tropospheric response to TSI are consistent with the study of Ortega et al. [94] supporting a weak role of TSI forcing in NAO reconstructions over the past millennium.

Conclusions
We investigated the response of the summer NA modeled climate to solar forcing in the preindustrial last millennium. We used solar-only forced simulations with the latest recommendation for a weak scaling of the TSI amplitude in CMIP5. These simulations were performed with the models MPI-ESM-P (single simulation) and CESM-LME (ensemble of four simulations). Our study included ensemble runs in order to identify the effect of model internal variability on the models climatic response to TSI forcing. This analysis is especially relevant for the comparison of model output to the results of empirical studies, as we investigated the summer season and use methodologies commonly applied in studies looking into the effect of external forcing on past climate. Our results show that:

•
The statistical methods used indicate similar results for the CESM tropospheric forced response, agreeing on a weak response with decreased SLP over Greenland (−0.1 hPa per TSI SD) and increased geo500 over the subtropics (+0.8 gpm per TSI SD). This response is weak at all lags studied, with the values of the correlation coefficient being r = 0.1 for SLP and r = 0.3 for geo500.

•
The solar signal was not robustly identified during the preindustrial period in the summer SSTs of the CESM ensemble mean over the NA basin. This result was indicated by linear methods and composite analysis, and further justified from the analysis of the PiControl simulations.

•
The investigation of the CESM individual ensemble members signified that, for the low TSI scaling, model internal variability is larger than the changes induced by TSI forcing in single model simulations regarding surface and tropospheric variables. This result is further supported by the SNR values being lower than 0.4 for the climatic variables studied.

•
The control simulations of the two ESMs used in our study indicate that spurious correlation patterns between climatic variables and TSI may arise with linear regression, Pearson correlation, and composite analysis when using single model simulations.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/atmos12050568/s1, Figure S1: Signal-to-noise ratio (SNR), calculated for the CESM forced simulations as the ratio of the ensemble spread to the ensemble mean standard deviation. The calculation follows the method described in Jerez et al. [95]. The SNR is shown for (a) SST, (b) SLP, and (c) geo500., Figure S2: Correlation patterns between TSI and the SST (column 1), the SLP (column 2) and the geo500 (column 3). The results are given in the first, second, third and fourth row for the realizations E1, E2, E3 and E4 of the CESM forced simulations. The hatched areas indicate statistical significance at the 5% level. Figure S3: Correlation patterns between TSI and the SST (column 1), the SLP (column 2) and the geo500 (column 3) for the MPI-ESM-P (row 1) and CESM (row 2) control simulations. The hatched areas indicate statistical significance at the 5% level. Figure S4: Correlation patterns between TSI and U200 for the MPI-ESM-P (column 1) and CESM (column 2) forced simulations. The hatched areas indicate statistical significance at the 5% level. The correlation is given in the first row for no lead-time. As the row number increases, the TSI lead time increases by one year., Figure S5: Correlation patterns between TSI and U200 for the MPI-ESM-P (row 1) and CESM (row 2) control simulations., Figure S6: Composite patterns for high TSI minus low TSI periods for the SST (column 1), the SLP (column 2) and the geo500 (column 3) for the MPI-ESM-P (row 1) and CESM (row2) control simulations. The hatched areas indicate statistical significance at the 5% level. Figure S7: Correlation patterns between the SST (column 1), the SLP (column 2) and geo500 (column 3) of the individual CESM ensemble members pairs, for the period 850-1849 AD. The parameters are 11-year low-pass filtered and the hatched areas indicate statistical significance at the 5% level. Figure S8: Composite patterns for high TSI minus low TSI periods for the SST (column 1), the SLP (column 2) and the geo500 (column 3). The results are given in the first, second, third and fourth row for the realizations E1, E2, E3 and E4 of the CESM forced simulations. The hatched areas indicate statistical significance at the 5% level.