Interactive comment on “ Internal and external variability in regional simulations of the Iberian Peninsula climate over the last millennium ”

Thank you for the kind comments and the generally positive view of the paper. Of course, all typos and suggestions regarding the redaction of the paper will be taken into account in the revised version of the paper. The bibliography references suggested by the reviewer will be checked out, and taken into account to enrich the discussion of the results. Regarding the main point of the reviewer, the discussion of the large-scale field in summer precipitation, it will be improved in the next version of the paper. However this


Introduction
The climate system fluctuates naturally over a large frequency range, from days to millions of years (Huybers and Curry, 2006).This variability is the combination of an undetermined level of internal variability superimposed on the net effect of a number of external forcings.Some of these forcings have a natural origin, such as changes in the so-lar irradiance or the radiative effect of big volcano events, whereas others have an anthropogenic cause, like land use changes or the emission of greenhouse gases and aerosols to the atmosphere during the Industrial Era (Crowley, 2000;Stott, 2003).
In the context of anthropogenic climate change, it is important to have available reliable estimations of the amplitude of natural variability on multidecadal timescales and at regional spatial scales, since this variability may hinder the attribution of trends observed to the anthropogenic forcing.In this respect, recent detection and attribution studies (Hegerl et al., 2011) have shown the fingerprint of external forcings in the temperature evolution of climate at continental scale during the last millennium.These kind of detection and attribution exercises require long historical records of climate variability.However, observations are too short to reliably assess multidecadal or even centennial climate variations, and therefore the analysis of past climate is a valuable tool in the estimation of the amplitude of climate variability and its mechanisms.
The efforts to assess the role of natural variability in climate evolution belong to two categories: climate reconstructions based on proxy indicators and climate model simulations.Comparing both approaches is important to identify systematic errors in simulations, as well as drawbacks in the methodologies used in reconstructions (González-Rouco et al., 2009).In particular, the validation of climate models in a past climate context may increase the confidence placed in climate change projections.Nevertheless, several limitations arise when comparing results from models and reconstructions.
Published by Copernicus Publications on behalf of the European Geosciences Union.
One drawback is the scale gap between both approaches.Although the use of comprehensive atmosphere-ocean global circulation models (AOGCMs) has become possible due to the impressive increase in computational power (Zorita et al., 2005;Tett et al., 2007;Ammann et al., 2007;Swingedouw et al., 2010;Jungclaus et al., 2010, among others), their spatial resolution is still too coarse to take into account the regional climate features caused by fine orographic details.Climate reconstructions are thought very sensitive to these details, which are implicit in the information extracted from the proxies.This scale gap may be bridged using Regional Circulation Models (RCMs), which simulate the climate system for a limited domain (Kittel et al., 1998;Jerez et al., 2010, among others).This downscaling approach is commonly used in climate change projections (Jacob et al., 2007;Déqué et al., 2007;Gómez-Navarro et al., 2010, among many others), but fewer studies exist nowadays focusing on the applications of RCMs in a paleoclimate context, partly due to the large computational cost involved (Diffenbaugh and Sloan, 2004;Zorita et al., 2010;Strandberg et al., 2011).In particular, most regional simulations available to date are too short to address the role of internal variability, or if long, they do not consider different runs (Gómez-Navarro et al., 2011).
Another important limitation in the model-proxy comparison is the inherent internal variability of climate models.Just like the actual climate, the models are affected by a strong chaotic internal variability over a broad band of time scales.This implies that a complete agreement at interannual timescales should not be expected when comparing the temporal evolution of model simulations and reconstructions, even if both are perfect (Yoshimori et al., 2005).Several studies have previously analyzed the role of internal variability in climate simulations.Goosse et al. (2005) used an intermediate complexity model to perform 25 simulations of the last millennium.These simulations shared the same external forcings, only differing in their initial condition.This ensemble allowed them to detect the fingerprint of climate forcings in the long-term evolution of temperature at hemispheric and continental scale.However, they found that the role of internal variability becomes greater at regional and interannual timescales.Similarly, Servonnat et al. (2010) used an AOGCM to analyse the role of external forcings against internal variability, focusing on the evolution of surface temperature.They also pointed out the greater role of internal variability at regional scales, which may blur the effect of external forcings at sub-continental scale.However, due to the large computational cost involved, most studies so far do not consider the use of downscaling techniques to analyse explicitly the role of internal variability at regional scales, despite that climate proxies contain important information at these scales.The fingerprint of the forcings in precipitation is other aspect which is seldomly considered in the bibliography.
Thus, in this study we present a comparison of two simulations performed with a climate version of the mesoscale model MM5 driven by the AOGCM ECHO-G over the last millennium  for a domain encompassing the Iberian Peninsula (IP).The model configuration and the external forcings are the same in both simulations.The only difference lies in the initial condition used to run the two simulations in the global model, and thus these experiments allow us to investigate the role of external forcing in the evolution of several climate variables, not only temperature, compared to the magnitude of the internal variability of the model at regional scale.In particular, we focus on the evolution of near-surface air temperature (SAT) and precipitation (PRE) in winter (mean of December-January-February) and summer (mean of June-July-August).

Description of the simulations
The global model ECHO-G driving the RCM consists of the spectral atmospheric model ECHAM4 coupled to the ocean model HOPE-G (Legutke and Voss, 1999).The model ECHAM4 was used with a horizontal resolution T30 (∼3.75 • × 3.75 • ) and 19 vertical levels.The horizontal resolution of the ocean model is approximately 2.8 • × 2.8 • , with a grid refinement in the tropical regions and 20 vertical levels.Two simulations were performed with this model configuration, both driven by the same reconstructions of several external forcings, which are depicted in the first two panels of Fig. 1.Cyan, pink and grey lines represent the evolution of atmospheric carbon dioxide, nitrous oxide and methane, respectively.The orange line represents the reconstruction employed for the variability of the total solar irradiance (TSI).Black lines show the estimated equivalent reduction in solar irradiance at the top of the atmosphere caused by volcanic eruptions.The sum of both lines is the effective solar irradiance, which is implemented in the model to take into account both sources of short-wave external forcing.The two simulations (hereafter referred as ERIK1 and ERIK2) cover the last millennium almost entirely, only differing in their initial condition (ERIK2 starts from a colder climate state).A full description of these simulations can be found in González-Rouco et al. (2003); Zorita et al. (2005); and references therein.
The RCM used to downscale the two AOGCM simulations is a climate version of the Fifth-generation Pennsylvania-State University-National Center for Atmospheric Research Mesoscale Model (Dudhia, 1993;Grell et al., 1994;Montávez et al., 2006;Gómez-Navarro et al., 2010).A double-nesting scheme was implemented, with a lowerresolution (90 km) outer domain, covering Western Europe, and a higher-resolution (30 km) inner domain covering the IP.Both domains are two-way coupled between them, whereas the RCM is one-way coupled to the driving GCM.These domains, as well as the chosen physical parametrisation, are the same as those previously described by Gómez-Navarro et al. (2011).The regional simulations have been driven with exactly the same external forcing reconstructions as in the global model simulations to avoid physical inconsistencies.The two downscaled simulations, which cover the period 1001-1990, will be referred hereafter as MM5-ERIK1 and MM5-ERIK2, respectively.
Here we analyse the role of internal variability at regional scales trough in a model-model comparison.Thus, the evaluation of the model skill in reproducing observations or climate reconstructions is beyond the scope of the present study.For a full description of the capabilities of this model configuration to reproduce a realistic climate over the IP, the reader is referred to Gómez-Navarro et al. (2011).

Correlation as a measure of internal variability
Figure 1 depicts the evolution of the anomalies in SAT and PRE in winter and summer averaged over the IP in the two RCM simulations.The eight series have been smoothed with a 31 yr running mean in order to filter out the high frequency signal and highlight the low-frequency variability.In the SAT series we can easily identify three main periods: a warm initial condition up to roughly 1400, followed by a long cold period that finishes around 1850, reverting to a stronger trend towards a warmer climate.These periods match the already known characteristic periods in the last millennium such as the Medieval Climate Anomaly, the Little Ice Age and the warmth of the Industrial Era.In particular, there are some noticeable cold periods such as the Spörer Minimum (around 1450), the Maunder Minimum (around 1700) and the Dalton Minimum (around 1810) which can be identified in the SAT series and are associated with a simultaneous reduction in the TSI and volcanic activity.The final trend can be linked to an increase in the TSI, but also with the large increase in GHG concentrations.In general terms, the impact of the external forcing in these cold periods seems to be more relevant in summer in both simulations, which show similar characteristics in terms of variability and temporal evolution.
The fingerprint of the forcings in the evolution of precipitation is, however, less apparent.Other than a slight decrease in precipitation, more noticeable in summer, since 1800 to the end of the simulation, it is not easy to identify the impact of external forcings on precipitation.Further, although the variability depicted by both simulations is similar in both seasons, their temporal agreement is lower than in the case of SAT, suggesting a stronger independence of precipitation of the external forcings.The variability of winter precipitation is stronger, as it corresponds to wetter conditions in this season in the IP, a feature correctly reproduced by both simulations.
Careful comparison of these two millennial simulation allows the study of the internal variability by quantifying the signal-to-noise ratio.This can be assessed through the calculation of the correlation of the temporal series associated to several variables in the two experiments.A simple conceptual model of the evolution of a climate variable that is partially driven by the external forcing may illustrate this point.The variable can be considered as a combination of the external forcing plus a contribution of noise due to the inherent chaotic nature of the simulation: where T is the variable of interest, α is a proportionality constant, f is the direct effect of the external forcing in the variable and W is a normalised random variable representing internal variability, which is uncorrelated with the forcing.If we now perform several identical simulations, only differing in the initial conditions, we have several variables T i .The point here is that random noise prevents correlation between these variables from being perfect.Instead, the correlation between these variables can be depicted as: where the variance of the variable is assumed to be the same in all the experiments (i.e.Var(T i ) = Var(T j ), ∀i,j ) and W and f are uncorrelated.Hence, according to the last equation, if forcing plays a strong role in the evolution of T , α can be considered negligible and the correlation is close to one.
On the other hand, if the evolution of the variable depends strongly on the internal variability (α is large), the right term in Eq. ( 2) becomes small and the evolution of the variable is not correlated between the different experiments.Thus, the correlation gives a quantitative measure of the relative role of internal variability in the evolution the different variables of the simulation, a parameter that will be used in the next section.
It is important to note, nonetheless, that the influence of the external forcing detected in this way may be dependent on the time scale, season and area.This is due to the different amplitude of internal variability, from daily to interdecadal, and to the amplitude of the external forcing at different time scales.In this study we have performed an analysis using different running means with increasing time intervals in an attempt to identify the temporal scales at which the signalto-noise ratio is stronger.Low-pass Fourier filters have also been employed to retain the low-frequency variability, which is more clearly governed by external forcings.
Another important aspect to consider is what correlation limit can be safely considered as statistically significant.Indeed, the running mean (or the Fourier filters) we apply to smooth the series introduces an artificial autocorrelation which could affect the calculations of correlation and overestimate its significance.Thus, the estimation of the confidence interval has to take into account this fact.We use a statistical test which calculates the correlation threshold numerically (Ebisuzaki, 1997), taking into account explicitly the artificial autocorrelation structure introduced in the smoothing process we apply to the series.

Forced vs. unforced evolution of temperature and precipitation
Figure 2 depicts the spatial distribution of the correlation between the SAT series simulated in the MM5-ERIK1 and MM5-ERIK2 experiments, both smoothed using three different running means.The evolution of this variable in the two experiments is highly correlated in both seasons (black circle denotes grid points where the correlation is significant at the 95 % confidence level according to a bootstrap method (Ebisuzaki, 1997)), and it is quite steady over the domain.
Correlation in summer is slightly stronger, and in both seasons tends to be somewhat greater when a stronger smoothing is applied.It is interesting to note that the spatial structure of the correlation is different in winter and summer, with stronger correlations in the southern part of the IP in summer and stronger correlations in the north in winter.There are hardly any changes when different running means are applied.The main conclusion we can derive from this figure is that the evolution of SAT seems to be dependent on the evolution of the external forcings, and this result is valid everywhere in the IP.This result is quite reasonable, since nearsurface temperature should be physically strongly modulated by the external forcing.The magnitude of the correlation is essentially independent of the degree of smoothing, which may be interpreted as the fact that the strongest influence of the external forcing is already attained at the smallest time scale probed here.Similarly, Fig. 3 depicts the same information for precipitation.In this case the correlation is in general lower, and, in some cases, even negative.There are, however, some well defined areas where this variable still exhibits a high and statistically significant correlation between experiments (denoted with a black circles, as in Fig. 2 ).As in Fig. 2, the correlation structure is different in winter and summer, and given a season, it hardly changes when a stronger smoothing is applied.In winter, the areas which most clearly respond to Correlation map of precipitation series in winter (December-February, left column) and summer (June-August, right column) between the MM5-ERIK1 and MM5-ERIK2 experiments.The three rows represent the correlations performed with the series smoothed through a running mean of 31 (top), 61 (middle) and 91 (bottom) years, respectively.Black circles denote grid points where the correlation is significant at the 95 % confidence level according to a bootstrap method (Ebisuzaki, 1997).
the forcings are the northeast of the IP and the main mountain systems such as the Pyrenees, and the Iberian and Betic systems.Conversely, the sensitive areas in summer are located in the west and north of the IP, and show no clear influence of the orographic features of the domain.We discuss a possible physical explanation for these patterns in the next subsection.
Given the length of the simulations, we have tried another test to evaluate the significance of the correlations.We divided the 990-yr long series in two subseries of 495 yr.We calculated the same correlations as in Figs. 2 and 3 in the two periods to check whether they are period-dependent.We found that the correlation patterns hardly change in different periods (not shown), which further suggests that the correlations found have a physical meaning rather than just being a statistical artifact.
When we further analysed the significance of the correlations calculated in the former section and their relationship with the smoothing applied to the series in Fig. 4, the mean correlation for SAT is above a 95 % confidence interval in both seasons and for all the tested running means, although in winter it tends to be slightly lower.The resemblance between series increases with longer smoothing, but it saturates around 50 yr.The limit for the confidence interval increases monotonically with longer smoothing, and is greater for summer when correlations are also higher.Apart from the nonsmoothed series of winter temperature, nearly all grid points exhibit a correlation above the 95 % confidence level, supporting our previous interpretation on the importance of the external forcings in the evolution of SAT in the IP during the last millennium.
From a regional modelling perspective, the results for precipitation are especially interesting.Figure 4 illustrates how the correlations for this variable are, on average, below the confidence level.Nevertheless, as mentioned above, there are areas where the correlation is still high.The dashed bar represents the percentage of grid points in the domain where the correlation is significant at the 95 % confidence level, and Fig. 3 allows their identification in different seasons.It is in these areas where the forcings play an important role in the evolution of precipitation.They can only be identified through the use of a high resolution model, since the average spatial process dilutes the statistical confidence, as Fig. 4 clearly illustrates.Another important aspect of these calculations is that they demonstrate that, although the correlation between experiments tends to be higher when longer smoothing is applied, the threshold for statistical significance is also larger, so that the number of grid points that show a significant correlation does not increase monotonically.

Physical meaning of the correlations
Although statistical significance is a necessary condition, it is not sufficient to assert that there is a causal relationship between the long-term evolution of forcings and SAT and PRE.The physical link between forcing and temperature is straightforward: the stronger the external radiative forcing, the higher the temperature.It depends neither on regional features nor on the season.The correlation map at global scale (not shown) between the external forcings and SAT is homogeneously positive over most of the globe, independent of the season, and it is similar in magnitude to that shown in Fig. 2 for the IP.
However, the link between forcing and precipitation is less obvious.In principle, higher temperatures tend to increase the evaporation, and hence the moisture content of the atmosphere.However, higher temperatures would tend to reduce the relative humidity for a given level of moisture, which tends to diminish the cloud cover and precipitation.The net result may depend on the regional features, the large-scale circulation or the season.In fact, we found that the largescale correlation between forcings and precipitation shows no clear homogeneous signal over the global (not shown), as in the case of SAT.Thus, the net effect of higher/lower forcings is dependent on other indirect factors such as modifications in the local circulation or the interaction with the orography.In addition, this relationship strongly depends on the season, as illustrated by Fig. 3.For this reason, we have investigated the physical mechanism linking the evolution of precipitation and forcing separately for winter and summer, and focusing only in the IP, since in other areas it could be different.In the remaining part of this section we focus on the low-frequency variations, since they show a larger signal-to-noise ratio.To do so, we use a simple lowpass Fourier filter, which eliminates the variability timescales shorter than 50 yr.We analyse the low-frequency variations of SAT and PRE separately through an Empirical Orthogonal Function (EOF) analysis (Hannachi et al., 2007), a methodology that reduces the high dimensionality of complex phenomena, such as climate, and has been used in other studies regarding long regional climate simulations (Gómez-Navarro et al., 2010).Finally, since we are interested in the variations relative to the mean state, and the precipitation over the IP is strongly heterogeneous (Serrano et al., 1999), we have used standardised precipitation series to avoid an overrepresentation of the wettest areas in the northwest of the IP.Hence, all EOF maps shown here are dimensionless.
The first EOF of the normalised low-frequency variations of SAT and PRE in summer in the experiment MM5-ERIK2 are shown in Fig. 5a and b and explain 89 % and 44 % of the variance, respectively.The corresponding figures for MM5-ERIK1 are very similar and are not shown here.The associated Principal Components (PCs) for the two variables and experiments are shown in Fig. 5d, together with the external forcings.The close relationship between the evolution of SAT and PRE is apparent in the two simulations.The correlation between the PCs of SAT and PRE is 0.82 and 0.79 for the MM5-ERIK1 and MM5-ERIK2, respectively.Figure 5c shows the correlation map between the low-frequency evolution of the two variables for MM5-ERIK2 (the corresponding map for MM5-ERIK1 is similar and has also been omitted), and independently illustrates the close link between the two variables.The resemblance between the maps in Fig. 5b and c can be better understood by looking at the PCs.The low-frequency variations of SAT are dominated by a spatially homogeneous EOF, whereas the variations of precipitation display several spatial characteristics.The strong correlation between the PCs, as well as the high percentage of variance that the first EOF for each variable explains, drive the clear correlation between these variables and their spatial structure.
On the other hand, the correlations between the PCs of the two experiments are 0.81 and 0.34 for SAT and PRE, respectively.This is in good agreement with our previous finding of a stronger influence of the external forcings in the evolution of SAT than in the case of PRE.In addition, the similarity between the map in Fig. 5c (or Fig. 5b) and correlation maps in Fig. 3 for summer is clear.Again, the explanation for the structure and intensity of these correlation patterns is better sought in the EOF analysis.The homogeneous first EOF, together with the large amount of variability it explains and the large correlation between the associated PCs, force a high and homogeneous correlation between the SAT in the two experiments.On the contrary, the lower correlation between PCs associated with the evolution of precipitation in the two experiments precludes strong coupling between them.Despite this, the shape of the main variability mode is similar in the two simulations, and the correlation between both PCs is not negligible (0.35, statistically significant), which explains why the areas most affected by this pattern stand out in the correlation maps of Fig. 3.
Having identified that the response to the forcing in summer precipitation over the IP is due to the main variability mode, which has the same spatial structure and is clearly correlated in the two simulations, we sought the physical mechanism behind this link.First, we separately considered largescale and convective precipitation.The correlation maps equivalent to those shown in Fig. 3 for convective precipitation alone, depict very low values and no spatial structure, whereas the maps corresponding to large-scale precipitation look very similar to those in Fig. 3 (not shown).This suggests that the response to forcing is in the large-scale field, in particular in the response of regional circulation to external forcings.To confirm this, Fig. 6a shows the correlation All calculations were performed with the smoothed low-frequency series.Black circles denote grid points where the correlation is significant at the 95 % confidence level according to a bootstrap method (Ebisuzaki, 1997).
between the low-frequency filtered series of external forcings and SLP for summer in the ERIK1 experiment (the map corresponding to ERIK2, not shown, exhibits the same pattern and supports the same physical explanation).We show the calculations in the AOGCM fields since this area lies outside the domain simulated by the RCM, and, in any case, the differences between the RCM and the AOGCM in the SLP field are small.This figure shows how the local circulation is affected by the driving forcings.In particular, the strengthening of the Azores high reduces the large-scale precipitation over the northwest, whereas the low over Morocco has the opposite effect and is responsible for the simultaneous increase of precipitation in the southeast of the IP and over the Mediterranean Sea, which is precisely the shape of the first EOF shown in Fig. 5b.
We analyzed the correlation pattern depicted in Fig. 6a.It is closely related to the leading variability mode of simulated summer SLP, extracted through and EOF analysis.In fact, the spatial correlation between the first EOF and the correlation pattern in Fig. 6a is 0.68, whereas the percent-ages of variance explained by these patterns are 35 % and 24 %, respectively.This variability mode is not only present in the simulation, but also contributes significantly to variability in observations.Indeed, we found that the percentage of variance explained by this pattern in an observational dataset for SLP in summer for the 20th century (Trenberth and Paolino, 1980) is up to 18 %.To show at what extent this pattern is sensible to forcings, we have projected this pattern onto the SLP field, considering not only the 990-yr simulation, but also including a continuation of the ERIK1 simulation under the SRES scenario A2 of the 4th IPCC.This calculation defines a time series of the intensity of this mode through the simulation.This series is shown in Fig. 6b, together with the 31-yr running mean (the index corresponding to the projection of the correlation pattern onto observations in the 20th century is also shown in green for comparison).It becomes clearly stronger under the climate change projection, when the external forcings is especially intense, and illustrates how this large-scale variability mode is stimulated by climate forcings.It is noticeable that the same mode is Fig. 6.Top: (a) correlation between summer (June-August) SLP and external radiative forcings in the experiment ERIK1.We show the results in the AOGCM field, since this area lies outside the simulated domain by the RCM.The calculations were performed with the smoothed low-frequency series.Bottom: (b) evolution of the climate index defined as the projection of the correlation map in (a) onto the summer SLP field (red line), together with its 31-yr running mean (black thick line).The projection of this pattern onto the observational dataset by Trenberth and Paolino (1980) is also shown (green line).Grey shadow indicates the time period where the prescribed forcings change from the paleosimuation to the future climate scenario SRES A2, with natural forcings kept constant.excited in the paleosimulation and in the climate change projection, although the forcings during the last millennium are mostly solar variability and volcanic activity, whereas in the A2 scenario these natural factors are kept constant and only changes in greenhouse gases are prescribed.
The results for winter are different.As discussed above, Fig. 3 shows how the correlation structure for winter is different from that of summer, which suggests that the underlying physical mechanism may also be different.Figure 7a, b and c are equivalent to Fig. 5a, b and c for winter (as before, the following argument is based on the MM5-ERIK2 experiment, although it also holds for the MM5-ERIK1 experiment).The evolution of SAT in winter is dominated by a homogeneous EOF, with associated PCs that are strongly correlated in the two experiments (correlation 0.68).However, an important difference between summer and winter is that in the latter, the leading EOF for precipitation (Fig. 7b) is very different from the correlation pattern between temperature and precipitation (Fig. 7c).The explanation for this difference has to be sought in the impact of large-scale circulation in the precipitation regime in winter over the IP.The North Atlantic Oscillation (NAO) is a variability mode of SLP in the North Atlantic area which strongly affects the winter precipitation amount in the IP, especially in the western parts (Trigo et al., 2004).The MM5-ECHO-G model is able to successfully reproduce this feature (Gómez-Navarro et al., 2011).Figure 7d illustrates this by showing the correlation between the NAO index, defined as the PC associated to the first EOF of winter SLP in the North Atlantic area, and the precipitation series in each grid point.It is apparent how the areas most affected by NAO (Fig. 7d) are those standing out in the first EOF of precipitation (Fig. 7b).This means that in winter the main variability mode of precipitation is dominated by NAO variations.However, the correlation of the low-frequency variations of NAO in the two simulations is only 0.17 (below 0.27, the significance level at the 95 % confidence level), indicating that this important circulation mode does not respond to the external forcing, but is dominated by internal variability in the AOGCM.This explains the generally lower impact of the driving forcing in the evolution of precipitation in winter.In winter, the fingerprint of the relationship between SAT and PRE illustrated in Fig. 7c has to be sought in the second EOF.This is shown in Fig. 7e and explains 18 % of the lowfrequency precipitation variance.We compare this second most important variability mode for winter precipitation with the correlation pattern between SAT and PRE.This mode of precipitation variability, once the influence of the non-forced NAO has been removed, responds to the forcings, as can be identified by comparing Fig. 7e and the correlation maps for winter in Fig. 3.In fact, the correlation between the PC associated to this precipitation mode and the PC associated with the SAT is 0.62 (0.66 in MM5-ERIK1).
The physical link between temperature and precipitation in winter described above, which is responsible for the resemblance between the evolution of the precipitation series during this season in the two simulations shown in Fig. 3, is not due to the response of the large-scale circulation to the external forcings, as in the case of summer.Instead, it is due to interactions between the large-scale circulation and the orography of the RCM. Figure 8a shows the correlation between the low-frequency evolution of precipitation series in the two simulations, together with the orography considered by the model (in green contours).The correlation is more intense near the main mountain systems such as the north side of the Pyrenees or the western part of the Iberian and Betic systems.This figure, together with the characteristic anti-correlation map between temperature and precipitation shown in Fig. 7c, suggests that there may exist a modulation in the condensation level, driven by temperature variations, which could affect especially the precipitation over mountain regions.To check this, the condensation level is estimated using the approximation of the difference between temperature T and dew point temperature T d CL ∝ (T − T d ).
(3) Both variables measured in the surface Lawrence (2005) present a modern review on the relations between moisture, temperature and how they are related.Figure 8b shows the correlation between temperature variations and the temperature differences in Eq. ( 3).There is a generally strong and positive relation between the height of the condensation level and the temperature (which is not an obvious result since the moisture content over the IP depends to a large extent on the evaporation rate in the Atlantic Ocean, which is lower in colder periods).The relation is stronger near the main mountain systems, although the correlation map also depicts sensitivity to the Atlantic flow on the windward side of the mountains.That is, in cold periods the condensation level sinks, especially on the windward side of the mountains, and this flavours the increase of precipitation in these areas.Hence, the physical link between precipitation and forcing in winter is through variations in the condensation level, directly modulated by variations in temperature.Surprisingly, despite the intensity of the noise due to internal variability in the winter precipitation, this mechanism is strong enough to leave an observable mark in the amount of precipitation in some areas over the IP, characterised by the orography.It is important to note that this mechanism can only be accurately reproduced within the context of a high resolution simulation capable of resolving the fine spatial scales involved.

Summary and conclusions
In this study we have compared the evolution of SAT and PRE in two millennial paleoclimate simulations performed with a RCM with a spatial resolution of 30 km for the IP.The comparison allows us to evaluate the importance and magnitude of the internal variability in the evolution of these variables relative to the influence of the reconstruction of external forcings used to drive the simulations.The underlying argument is that if internal variability dominates the evolution of a given varible, the temporal correlation of the series associated to it in both experiments would be negligible.
The results indicate that the long-term evolution of SAT is strongly affected by the external forcings driving the simulation.This variable responds homogeneously to the external factors over the IP at most temporal scales.The evolution of precipitation is, however, more strongly governed by chaotic variability at regional scale.In particular, there are few areas in the IP, the main mountain system in winter and the north and west areas in summer, where the precipitation is significantly driven by external forcings.However, in many parts of the domain the influence of the external forcing can not be detected in the evolution of precipitation.It is important to note that the significance of the correlation emerges at regional scales, and is blurred when a spatial average is performed.This stresses the importance of high resolution simulations in exercises comparing the model results with proxy reconstructions of precipitation.
The influence of the external forcing on precipitation is especially weak in winter.This is due to the nature of the winter precipitation over the IP, which is dominated by variations in the NAO.The NAO seems to be quite insensitive to the external forcing in the simulations at the investigated timescales.Once the NAO signal is removed from the precipitation series, the leading variability pattern corresponds quite well with the areas which are more clearly able to respond to forcings.Summer precipitation is overall more strongly affected by variations in the forcing.Its main variability mode matches well the areas where summer precipitation responds to external forcing.In fact, we have been able to demonstrate that this precipitation mode is dominated by modulation of the large-scale SLP by the external forcing in summer.This is in contrast with the stronger contribution of the internal variability of SLP in winter.During this season, there are still some areas where precipitation responds the forcing, but, in this case, it is not through modification in the large-scale flow but through the interaction between condensation level and orography.
Our findings regarding the impact of internal variability in the simulations may have a strong impact on how comparison between simulations and reconstructions are performed.In particular, we have been able to identify areas where we should not expect good agreement between the model and the reconstructions, even if both are perfect.On the other hand, there are areas, mostly in the main mountain systems, where mismatches between both approaches can not be argued to be due to internal variability.These results stress the importance of RCMs in paleoclimate studies, since we have demonstrated that the physical mechanism responsible for the response of precipitation to external forcings in winter can only be realistically reproduced by using high resolution simulations.
It is important to note that some of these findings may be model-dependent.Different global models develop slightly different circulation patterns.At regional scale, interactions with orography strongly depend on the set of parametrizations employed (this is especially true for precipitation (Fernández et al., 2007)).Thus, the mechanisms proposed here may suffer modifications if different model configurations are employed.However, it is hard to address a priori these important issues with only two simulations.A larger ensemble of runs with different driving global models and set of parametrizations would be required to reinforce our findings, but it is presently computationally prohibitive.Thus, it is beyond the scope of this study to address these uncertainties.
A further important comment has to be made regarding the reconstructions of solar forcing used in these simulations.The evolution of the TSI used in these simulations is from Crowley (2000).However, a more recent reconstruction of this variable (Krivova and Solanki, 2008) depicts a much smaller amplitude of the variations.In particular, these authors estimate a difference in total solar irradiance between the Late Maunder Minimum and late 20th century of 1.25 W m −2 (about 0.09 %), whereas the past solar irradiance used in these simulation changes by 0.3 %.Yet an even more recent reconstruction of TSI over the Holocene (Shapiro et al., 2011) again points to an even wider amplitude of TSI variations than those used in these simulations: 0.4 % change between the Late Maunder Minimum and late 20th century.It is beyond the scope of this study to analyse which of these reconstructions more realistically represents the past.The influence on our analysis of using reconstructions with lower amplitude in the simulation would be to reduce the term Var(f ) in Eq. ( 2), and thus reduce the correlation of the same variable between simulations.Higher-amplitude reconstructions of past TSI would have the opposite effect.
A similar argument could apply to volcanic forcing, for which the uncertainties are still also large.In addition, the implementation of volcanic forcing in these simulations, simply as a reduction in the effective solar irradiance, possibly precludes a more realistic simulation of the volcanic winter warming at mid and high latitudes due to the NAO response to differential effect of volcanic aerosols in the Tropics and high latitudes (Stenchikov et al., 2006;Fischer et al., 2007).According to this mechanism, winters after volcanic eruptions should experience a stronger NAO and thus the IP would tend to receive less precipitation.This mechanism tends to increase the influence of the volcanic forcing and it would phase-lock the simulated precipitation in both simulations in periods with intense volcanic activity more strongly than is simulated by our model set-up.
Finally, some of these conclusions can be extended with caution to the climate change projections.In a forced scenario, the SAT can be expected to be influenced by external forcings, and hence their projections present a reasonable degree of confidence.Evolution of precipitation is nevertheless less reliable since its behaviour at regional scale is governed by greater uncertainty due to the influence of internal variability.This drawback is an addition to the well known uncertainties characteristic of precipitation projections under climate change scenarios (Christensen et al., 2007).Nevertheless, it has to be taken into account that our findings depend on the intensity of the external forcings.In climate change projections the intensity of external forcings is stronger, and thus the overall role of internal variability can be expected to be lower.In future works, a similar study should be carried out with different runs for the same climate change projection.

Fig. 1 .
Fig.1.Evolution of the external forcings (first two panels), spatial average of SAT (two middle panels) and PRE (bottom panels) anomalies for winter (December-February) and summer (June-August).The MM5-ERIK1 and MM5-ERIK2 experiments are coloured red and blue, respectively.The eight series have been smoothed with a running mean of 31 yr, and the same vertical scale has been applied to both seasons to emphasise their different variability.

Fig. 2 .
Fig. 2. Correlation map of SAT series in winter (December-February, left column) and summer (June-August, right column) between the MM5-ERIK1 and MM5-ERIK2 experiments.The three rows represent the correlations calculated with the series smoothed by a running mean filter of 31 (top), 61 (middle) and 91 (bottom) years, respectively.Black circles denote grid points where the correlation is significant at the 95 % confidence level according to a bootstrap method (Ebisuzaki, 1997).

Fig
Fig. 3.Correlation map of precipitation series in winter (December-February, left column) and summer (June-August, right column) between the MM5-ERIK1 and MM5-ERIK2 experiments.The three rows represent the correlations performed with the series smoothed through a running mean of 31 (top), 61 (middle) and 91 (bottom) years, respectively.Black circles denote grid points where the correlation is significant at the 95 % confidence level according to a bootstrap method (Ebisuzaki, 1997).

Fig. 4 .
Fig. 4. Significance of the correlations between the MM5-ERIK1and MM5-ERIK2 series of SAT (top panel) and PRE (bottom panel) when using different running means to filter out the high frequency signal.Blue (red) color represents the results for December-February (June-August).The shaded area for each variable is the threshold for the correlation at the 95 % confidence level, obtained through a bootstrap method(Ebisuzaki, 1997).Solid bars represent the mean correlation for the domain, whereas dashed bars represents the percentage of grid points in the domain, which show a significant correlation.

Fig. 5 .
Fig. 5. Evolution of the summer (June-August) SAT and PRE during the last millennium in the two simulations.The three top maps show the first EOF of the standardised SAT (a) and PRE (b), as well as the correlation between these two variables (c).All maps are dimensionless, and correspond to the MM5-ERIK2 experiment, since the corresponding maps for MM5-ERIK1 are very similar (not shown).(d)showsthe normalised PCs associated with the two variables in the two experiments, as well as the external radiative forcing.All calculations were performed with the smoothed low-frequency series.Black circles denote grid points where the correlation is significant at the 95 % confidence level according to a bootstrap method(Ebisuzaki, 1997).

Fig. 7 .
Fig. 7.The three top maps show the first EOF of standardised series of winter SAT (a) and PRE (b), as well as the correlation between these two variables (c).(d) shows the correlation between the precipitation series and the NAO index, defined as the PC associated with the leading EOF of SLP in the North Atlantic area.(e) is the second EOF of winter precipitation.All maps are dimensionless and correspond to the MM5-ERIK2 experiment, since the corresponding maps for MM5-ERIK1 are analogous (not shown).The calculations were performed with the smoothed low-frequency series.Black circles denote grid points where the correlation is significant at the 95 % confidence level according to a bootstrap method (Ebisuzaki, 1997).

Fig. 8 .
Fig. 8. Physical response of winter precipitation over the IP.(a) shows the correlation between precipitation series in the MM5-ERIK1 and MM5-ERIK2 experiments, and it is the analog to Fig. 3 but using a low-pass filter instead of running means.(b) depicts the correlation between the temperature and height of the condensation level in the MM5-ERIK2 experiment.All calculations were performed using the low-frequency series.Green contours indicate the terrain height above 500 m spaced by 300 m.Black circles denote grid points where the correlation is significant at the 95 % confidence level according to a bootstrap method (Ebisuzaki, 1997).