The importance of external climate forcing for the variability and trends of coastal upwelling in past and future climate

The eastern boundary upwelling systems, located in the subtropics at the eastern boundary of the Atlantic and Pacific oceans and mainly driven by the trade winds, are the major coastal upwelling regions. Previous studies have suggested that the intensity of upwelling in these areas in the past centuries may have been influenced by the external radiative forcing, for instance by changes in solar irradiance, and it will also be influenced in the future by the increasing atmospheric greenhouse gases. Here, we analyse the impact of the external climate forcing on these upwelling systems in ensembles of simulations of two Earth system models. The ensembles contain three simulations for each period covering the past millennium (900–1849) and the 20th century (1850–2005). One of these Earth system models additionally includes the near future (2006–2100). Using a set of simulations, differing only in their initial conditions, enables us to test whether the observed variability and trends are driven by the external radiative forcing. Our analysis shows that the variability of the simulated upwelling is largely not affected by the external forcing and that, generally, there are no significant trends in the periods covering the past and future. Only in future simulations with the strongest increase of greenhouse gas concentrations the upwelling trends are significant and appear in all members of the ensemble.


Introduction
Eastern boundary upwelling systems (EBUSs, including the California, the Canary, the Benguela, and the Peru upwelling systems) are highly productive coastal ocean areas where nutrient-rich, cold water upwells by the action of favourable winds.The easterly trade winds which dominate these subtropical regions cause Ekman transport from the coast to the open ocean, perpendicular to the wind stress forcing, leading to upwelling at the coast.In addition, the wind stress curl between the jet of the trades and the coast, induced by to the relaxation of wind speed towards the coast, causes upwelling further offshore.
It has been suggested (Bakun, 1990) that changes in the external climate forcing, such as greenhouse gases and solar activity can influence the intensity of coastal upwelling.Bakun's hypothesis states that surface temperature over land should warm faster than over the oceans under increased radiative forcing, leading to an intensification of the subtropical continental lows and the oceanic highs, to an intensified cross-shore air pressure gradient, and to a strengthening of the upwelling-favourable winds (Bakun, 1990).
Some observations over the 20th century (Narayan et al., 2010) and simulations of the 21st century (Wang et al., 2015) have been interpreted, according to this hypothesis, as indicative of upwelling intensification due to stronger external climate forcing in the recent past and future.Also, coastal sediment records covering the past millennium and indicative of upwelling (McGregor et al., 2007) have been interpreted as a response to past variations in the external climate forcing, mainly solar irradiance and volcanism, with weaker upwelling during the Little Ice Age (centuries around 1700 AD) and stronger upwelling during the Medieval Warm Period (around 1100 AD).Although not totally ascertained, it is generally assumed that the Little Ice Age was mainly driven by a weaker solar activity and increased volcanic activity, which lead to reduced radiative forcing (Hegerl et al., 2006).For the Medieval Warm Period, this attribution is not as clear, but proxy records indicate a stronger solar activity than during the subsequent Little Ice Age (Usoskin et al., 2004).
For the past millennium, sediment cores are used as a proxy for upwelling due to the lack of direct long-term observations (McGregor et al., 2007).However, the upwelling itself is not directly derived from the sediment cores, but rather indirectly inferred on the basis of the water temperatures, assuming that cooler temperatures are indicative of stronger upwelling.This assumption could lead to a misinterpretation because temperature changes could also have been driven by the direct effect of the radiative forcing itself, in addition to any possible influence of the radiative forcing on upwelling.
For the 20th century, a metadata analysis has found a significant intensification of wind stress in three of the four major coastal upwelling systems (California, Benguela, and Humboldt but not Canary) (Sydeman et al., 2014).Previous studies using observations, reanalysis, and model data have detected increasing upwelling intensity over the past century (Di Lorenzo et al., 2005;Gutiérrez et al., 2011;Santos et al., 2012), but others have not (Pardo et al., 2011).A possible cause may be due to insufficient data homogeneity in longterm wind station records and in meteorological reanalysis, which may blur the identification of long-term trends (Sydeman et al., 2014).Also, long-term records of upwelling intensity are indirect, and sometimes even based on wind records themselves (Bakun et al., 2010).
The effect of increasing greenhouse gas concentrations on upwelling regimes in future scenarios has been investigated by several authors.Wang et al. (2015) and Rykaczewski et al. (2015) analysed the simulated trend of upwelling in several simulations included in the Climate Model Intercomparison Project (CMIP5) (Taylor et al., 2012) driven by the representative concentration pathways (rcp) 8.5 scenario, a scenario with an increase of the globally averaged external radiative forcing of 8.5 Wm −2 by the year 2100.They found that in most of the EBUSs, models tend to simulate an intensified upwelling and longer upwelling seasons under climate change.In contrast, the studies of Mote and Mantua (2002) (21st century simulation) and Hsieh and Boer (1992) (double carbon dioxide simulation) found only a weak trend or no trend in future coastal upwelling.
The evidence for the effect of radiative forcing on upwelling is, therefore, still not conclusive.For instance, the Benguela upwelling system does not exhibit a long-term intensification in the recent decades (Bakun et al., 2010).This has been explained by the counteracting influence of El Niño-Southern Oscillation (ENSO) on humidity in the Peru and Benguela EBUSs regions, which would also influence the regional radiative forcing and modulate the land-ocean thermal contrast (Bakun et al., 2010).During an El Niño event, the Benguela region is less humid, which causes a reduction of the greenhouse effect due to less water vapour, the most important atmospheric greenhouse gas (Bakun et al., 2010).According to this explanation, the upwelling trend due to stronger external radiative forcing would be biased by a higher frequency of ENSO events in the recent past.
The verification of Bakun's hypothesis in the recent past and future is of major importance for the EBUSs.Even if Bakun's hypothesis is correct, it is not clear whether past variations in external forcing could have been strong enough to drive upwelling intensity beyond the range of variations caused by internal chaotic climate variability.Hence, an apparent agreement between the sign of predicted and observed trend cannot be considered a conclusive proof of Bakun's hypothesis until those observed trends cannot be attributed to the external forcing.The analysis of the recently available ensemble of climate simulations with CMIP5 models over the recent past can shed light on this question.Firstly, they comprise different climate periods in which the external forcing has initially varied little in the past millennium, then more strongly presently, and much more strongly in the scenario simulations.Secondly, the use of ensembles of simulations with the same model offers the advantage of a much easier identification of the effect of the external forcing than the analysis of a single simulation.If the upwelling trends and multi-decadal variability is mainly externally driven, therefore determining the upwelling variations, all simulations should show a similar evolution of upwelling over time, since all simulations have been driven by the same prescribed external forcing.In contrast, if upwelling variations and trends are mainly a result of stochastic internal dynamics, the evolution of upwelling in the different simulations of the ensemble will tend to be uncorrelated in time, and the long-term trends will show different signs in the different simulations.This reasoning is more formally explained in Sect.2.3.
Bakun's hypothesis includes a physical mechanism by which upwelling may be driven by the external forcing, namely by the intensification of the sea level pressure (SLP) difference between the continents and the oceans adjacent to the upwelling regions.We also test in the climate model simulations whether the variability and trends of this air pressure difference, and of the associated alongshore wind stress, is correlated across the members of the ensemble of simulations.Additionally, another mechanism by which the external forcing could affect upwelling is through the ocean stratification.Global warming would lead to a warmer ocean surface, generally increasing the stability of the water column, with a shallower thermocline.This would hinder the upwelling (Hsueh and Kenney, 1972).
In summary, the main objective of the analysis is to test whether the external forcing prescribed in a set of ensembles of climate simulations is strong enough to drive upwelling in the four EBUSs, and whether, accordingly, there are significant long-term trends in upwelling over the last millennium, the last 156 years, and the near future that appear common to all simulations.This is obviously closely related to Bakun's hypothesis.
It is not our goal here to investigate the detailed mechanistic chain by which external forcing drives upwelling mecha-nisms in the models used.This task would be in any case not meaningful, since the main result of our study is that the external radiative forcing has not varied strongly enough over the past millennium, including the 20th century, to drive the upwelling variability over this period.Even the simulations for the future do not indicate that a moderate greenhouse gas emission scenario would be strong enough to unmistakably drive upwelling intensity.
In the following section, we present the data and the methods used for this analysis.In the main sections we compare some model results with the few available, mostly indirect, observations of upwelling.Later, we analyse the link between external forcing and the upwelling itself and with the immediate drivers of upwelling.We finalize the study with a discussion and conclusion section.
2 Data and methods

Satellite data, models, and simulations
We analyse here ensembles comprising three simulations (r1, r2, and r3) of two different Earth system models: the Max Planck Institute -Earth system model (MPI-ESM) and the Community Earth System Model -Last Millennium Experiment project (CESM-LME).
The MPI-ESM includes the ECHAM6 model (Stevens et al., 2013) for the atmosphere, the MPI-OM model (Jungclaus et al., 2013) for the ocean, the JSBach model (Reick et al., 2013) for the vegetation, and the HAMOCC5 model (Ilyina et al., 2013) for the marine biogeochemistry (Giorgetta et al., 2013).The simulations of the MPI-ESM cover the periods 900-1849 (past1000, MPI-ESM-P), 1850-2005 (historical, MPI-ESM-MR), and 2006-2100 (future, MPI-ESM-LR).For the future, we analyse here three scenarios with different strength in greenhouse gas forcing, representative concentration pathways rcp2.6,rcp4.5, and rcp8.5,where the numbers indicate the anthropogenic radiative forcing in Wm −2 reached by the year 2100 (Taylor et al., 2012).The horizontal resolution of the atmospheric model of the MPI-ESM is about 1.9 • (spatially varying) with 95 levels for the historical period and 47 levels for the past1000 and future periods, whereas the ocean model resolution is approximately 1 • (past1000 and future) and 0.4 • (historical), including 40 ocean layers.The climate model MPI-ESM participated in the CMIP5 project (Giorgetta et al., 2013), contributing three simulations to the historical ensemble, three simulations to the future ensemble, and one simulation for the past1000 period (r1).
The simulations in each ensemble were driven by almost identical external forcing.Only the volcanic scheme in the past1000-r1 differs from the ones in the past1000-r2 and past1000-r3 experiments in the prescribed distribution of the volcanic aerosol size.This distribution is assumed to be lognormal, with a standard deviation of 1.2 µm in r1 and 1.8 µm in r2 and r3 (Jungclaus et al., 2014).The forcings prescribed in the past1000 (Masson-Delmotte et al., 2013) follow a standard protocol defined by the project CMIP5 and comprise orbital forcing (Berger, 1978), variability in solar irradiance (Vieira and Solanki, 2010), volcanic activity (Crowley and Unterman, 2013), greenhouse gas concentrations (Flückinger et al., 2002;Hansen and Sato, 2004;MacFarling Meure et al., 2006), and land-use changes (Pongratz et al., 2009).The external forcing of historical simulations comprise the same variables and are described in IPCC, Annex II (2013).The reconstructed changes in the solar forcing have been disputed over the past decade, but the CMIP5 protocol agreed that small solar variations are likely more realistic.In these simulations they are rather small, with an increase of about 0.15 % from the Maunder Minimum (around the year 1700) to present time.In the scenario simulations, only changes in anthropogenic forcing were prescribed according to the rcp scenarios (Taylor et al., 2012).
To test the sensitivity of our results with respect to the choice of model, we also analysed simulations with the CESM-LME model (Otto-Bliesner et al., 2016).This model comprises the Community Atmosphere Model version 5 (CAM5) for the atmosphere, the Parallel Ocean Program version 2 (POP2) for the ocean, the Community Land Model version 4 (CLM4) for the land, and the Community Ice CodE (CICE4) for the sea ice (Hurrell et al., 2013).The simulations with CESM-LME were conducted after the simulation time of the CMIP5 project and do not belong to that set of simulations.The CESM-LME ensemble comprises a large amount of simulations with different configurations of the external forcing and different initial conditions.Here, we analyse simulations driven by all external forcings, natural and anthropogenic (denoted as "allforcings" in the CESM-LME project (Otto-Bliesner et al., 2016) comprising the same set of external forcings as used in the MPI-ESM simulations (https://www2.cesm.ucar.edu/models/experiments/LME)).The CESM-LME covers the period 850-2005 with a horizontal resolution of the ocean of 1 • and a horizontal resolution of the atmosphere of 2 • (Otto- Bliesner et al., 2016).Unfortunately, no future simulations with this model are publicly available yet (Kay et al., 2015).
A recent paper by Small et al. (2015) investigated the relevance of model resolution for the question of the warm bias that global climate models usually display in the eastern ocean basins (Richter, 2015).This bias may be related to the simulated upwelling intensity, but also to other physical processes, in particular to model clouds.The analysis by Small et al. (2015) indicates that a high atmospheric horizontal resolution is of larger importance than the horizontal resolution of the ocean model.The most realistic wind stress curl and upwelling was simulated with a high-resolution nested regional atmospheric model coupled with an eddy-resolving ocean model, resulting in more intense values of wind stress and shifted towards the coast.
The model resolutions of the two Earth system models used here are similar to the resolution of the CMIP5 models used by Wang et al. (2015).The strength and location of the jet seems to be impacted by the atmospheric model resolution (Small et al., 2015).This study also indicated that an atmospheric resolution of at least 1 • would be preferable for the representation of the jet.Unfortunately, no higher-resolution version of these Earth system models is available.Thus, although ocean processes of spatial scales of a few kilometres are only imperfectly resolved, the resolution of these two Earth System Models is plausibly fine enough to realistically represent the basic relationship between the large-scale wind forcing and upwelling dynamics, including the shorelineparallel winds and the wind stress curl.However, fine-scale ocean structures, like coastal filaments and wave dynamics, cannot be realistic represented.Recalling that Bakun's hypothesis is formulated as the effect of external forcing on the atmospheric dynamics, the analysis of these simulations is an informative test of Bakun's hypothesis, until higherresolution, multi-centennial simulations with coupled models become available.
Within each ensemble of simulations (past1000, historical, and future), the simulations differ in their initial state.The initial conditions for each simulation have been randomly chosen from a long pre-industrial control run, after which the model is allowed to run for several centuries driven by the (constant) forcing of the year 850 until a state of quasiequilibrium is attained.From this point onwards the transient simulation starts with the prescribed temporally variable external forcing.
To validate the MPI-ESM simulations, we compared the simulated sea surface temperature (SST) with remotely sensed SST at 4 km resolution from the Advanced Very High Resolution Radiometer (AVHRR; Casey et al., 2010) version 5.0, a radiometer onboard the National Oceanic And Atmospheric Administration (NOAA) satellites.

Upwelling index
Upwelling intensity in each of the four upwelling regions is defined by the corresponding model variable vertical mass transport (wmo) at the ocean model layer at 52 m (MPI-ESM) and the vertical velocity (wvel) at 50 m (CESM-LME) (close to the modelled mixed layer depth), and spatially averaged over each of the upwelling regions: Benguela (8-30 • E, 15-28 • S), Peru (80-70 • W, 20-10 • S), California (130-110 • W, 20-50 • N), and Canary (20-10 • W, 20-34 • N).Upwelling indices are here defined as the seasonal means over the strongest upwelling season in each region.Seasonal averages are calculated using the standard seasons definition: December-February (DJF), March-May (MAM), June-August (JJA), and September-November (SON).In the MPI-ESM past1000 and historical simulations, the main upwelling season in Benguela is September-November and June-August in all other upwelling regions.In the future simulations, the main upwelling season in Benguela is also June-August.Strongest upwelling in Benguela takes place between July and October (Tim et al., 2015).The main upwelling season could be both of the standard seasons (JJA or SON), since the mean difference is small.In the CESM-LME simulations, the main upwelling season is September-November in Benguela and Peru and June-August in California and Canary.
As briefly mentioned in the introduction, upwelling directly at the coast is caused by the alongshore wind stress, whereas further away off the coast upwelling is driven by the wind stress curl.In the simulations analysed, the regions selected to define the upwelling indices present spatially coherent upwelling, i.e. the vertical velocities at 52 m depth in each grid cell are almost all positively correlated with upwelling at the grid cells most closely located to the model land grid cells.Only at the oceanward fringes of these boxes this correlation vanishes.Therefore, each geographical box encompass the set of grid cells that show a coherent coastal upwelling, although the physical mechanism that drive upwelling is, strictly considered, not the same.This coherency is, however, reasonable since the strength of the alongshore wind stress directly at the coast should be correlated with the wind stress curl further offshore, both being generally driven by the ocean-to-land sea level pressure gradient, as also assumed in Bakun's hypothesis.

Methodology
When not explicitly stated, all correlations have been calculated with linear long-term detrended series.Low-pass time filters of 10 years for the historical and future data and lowpass time filters of 30 years for the past1000 data have been applied before correlating time series to filter out the high frequencies to focus on variations of decadal and multi-decadal timescales.For the past1000 simulations, the filter is chosen to highlight the frequencies of variations in the external forcing.For the shorter historical period we are more interested in the long-term trends, since the external forcing, mostly anthropogenic, increases continuously.Applying a 30-year low-pass filter for a time series covering 156 years would leave out too many degrees of freedom to establish the significance of these trends.
For the analysis of the statistical significance of the longterm trends in the upwelling indices and the linear correlations between the ensemble members a two-sided significance level of p = 0.05 was adopted.All time series, and in particular the ocean time series, are affected by serial correlations that preclude the application of the usual statistical tests to establish the significance of trends or correlations.Here, we have applied a Monte Carlo-based method (Ebisuzaki, 1997).This method generates copies of the original time series with the same serial correlation properties, but uncorrelated in time with the original series.For instance, to establish the significance of correlation coefficients, 1000 surro-gate pairs of series are generated and used to generate an empirical distribution of correlation coefficients under the null hypothesis that the correlation is zero.In the case of linear trends, 1000 series are generated with zero long-term trend but the same serial correlation as the original series, providing an empirical distribution of trends under the null hypothesis of zero trend.
Here, we show results from three of the available simulations from the CESM-LME ensemble, as three simulations already confirm the results obtained with the MPI-ESM model.These two models used here are the only ones that provided ensembles of simulations with basically the same forcing and with different initial conditions, thus being suitable to detect the possible imprint of external forcing.
After assessing the realism of the connection between atmospheric drivers and simulated upwelling, we further evaluate whether the upwelling intensity in these regions displays any imprint of the external climate forcing.To estimate the ratio of forced variability to total variability, we build a simple statistical model of the variations of a climate record, which decomposes its variability into a sum of a forced component, related to the external climate forcing, and an internal component: (1) The same model can be applied to describe climate records simulated in two climate simulations driven by the same forcing and started with different and random initial conditions: where the forced components are by construction equal in both simulations (the prescribed forcing is equal) and the internal components are uncorrelated in time.Since the response to external forcing could be theoretically lagged by a number of years, the time step (t ) of the external forcing does not have to be necessarily the same as for the internal variability.The ratio between the variance of y f and the variance of y 1 , equal to the variance of y 2 , can be shown to be equal to the correlation between y 1 and y 2 : where x, y = t x(t)y(t).
The equation shows that, in two simulations that have the same forcing but different initial conditions, the ratio between the variance of the forced component and the total variance is the same as the correlation between both time series.If the correlation between y 1 and y 2 is weak, the variance of y f is much weaker than the variance of y i , meaning that the effect of internal variability is much stronger than the effect of external forcing.This is true because of the assumption that the variance of the internal components of both time series are uncorrelated.The forcing is externally prescribed and, thus, cannot be influenced by the initial conditions.This mathematical description is the base of the statistical analysis presented here.Correlating the three simulations of an ensemble shows us whether the external variability is an important component of the total variability.This statistical reasoning provides the rationale for the calculations of the correlations between the indices of different simulations within an ensemble, which is profusely used through the whole paper.

Representations of upwelling and its atmospheric drivers in the climate model
The patterns of mean SST in the four main upwelling regions in their season with maximum upwelling simulated in one of the historical simulations with the model MPI-ESM-MR are compared with the corresponding data derived from the AVHRR in Fig. 1, both calculated for the overlapping period .The SST patterns display, in general, colder temperatures directly along the coast, with warmer temperatures offshore.The climate model is able to replicate these SST gradients remarkably well, in spite of the coarse spatial resolution.However, in the cases of Benguela (Fig. 1e and f) and Peru (Fig. 1g and h), the climate model shows an evident warm bias relative to the satellite data (Richter, 2015).This indicates that the main mechanisms of upwelling are likely reasonably represented in the model, though other important mechanisms are not.The shape of the modelled seasonal cycle of upwelling (Fig. 2) generally agrees with what is known from observations, with upwelling being more intense in the boreal warm half year (Fig. 7 in Chavez and Messié, 2009).The mean upwelling is, however, more intense by about a factor of 2 in the lower-resolution past1000 simulation (horizontal resolution of the ocean of 1 • (past1000) compared to 0.4 • (historical)), indicating that ocean model resolution is important to simulate the correct mean upwelling intensity, which may be relevant for the SST bias in the EBUSs simulated by many climate models (Wang et al., 2014).A reason for this difference is not clear cut but it is also not essential for this analysis.
The link between each upwelling index and wind stress, calculated as the patterns of correlations between the upwelling field and spatial averaged alongshore wind stress (Fig. 3 right panels), is also very realistic, indicating in each case that the alongshore wind stress favours Ekman transport, the main mechanism causing coastal upwelling (Tomczak and Godfrey, 2003).The alongshore wind stress was calculated as the wind stress over the ocean in the upwelling areas (for Benguela between 15-40 • S covering the whole Benguela upwelling region, North and South Benguela), with an angle of the shoreline of 15 • for Benguela and 45 • for the other regions.The relative standard deviation of the upwelling (Fig. 3 left panels) indicates the low variability of the more offshore regions.Therefore, the negative correlations in the more offshore regions, especially in the Canary upwelling system, do not account for a large portion of the upwelling index.A principal component analysis confirms these results.Correlation coefficients between the leading principal component and the spatial mean of the upwelling are between 0.80 and 0.92 in the four upwelling regions.Thus, the first principal component which explains most of the variability of the upwelling (more than 60 %) is highly correlated with the spatial mean.The correlation patterns between each upwelling index and the SLP realistically show regions with negative correlations over land and positive correlations offshore (Fig. 4), indicating that the across-coastline SLP gradient, with an intense oceanic subtropical high and an intense continental low, is conductive for alongshore wind stress through the geostrophic relation.In the last millennium, the external climate forcing drove changes in the climate.The main external climate drivers over this period were volcanic forcing, solar variability, greenhouse gases, and land-use changes (Hegerl et al., 2006) www.ocean-sci.net/12/807/2016/Ocean Sci., 12, 807-823, 2016 and, at millennium timescale, the slowly varying orbital configuration of the Earth (Laskar et al., 2004).These forcings had an effect on the global mean surface temperatures, both in climate simulations driven by these forcings and in proxy-based climate reconstructions (Fernández-Donado et al., 2013).Schurer et al. (2014) found that the variability of the surface air temperature of the Northern Hemisphere was mostly driven by the greenhouse gas concentration and volcanic eruption in the last millennium.There was a period of relatively high temperatures, the Medieval Warm Period with its maximum at around 1000-1300 AD.These high temperatures could have been caused by high solar and low volcanic activity but there is no consensus about the driver of the Medieval Warm Period, and internal variability could have played a role too.The Medieval Warm Period was followed by the Little Ice Age, a period of low temperatures caused by low solar and high volcanic activity.After the Little Ice Age, global mean temperatures rose mainly as result of higher concentration of greenhouse gases (Fernández-Donado et al., 2013).As illustration, Fig. 5 displays the multi-decadal global mean temperature evolution simulated in the ensembles of simulations analysed in this study.These temperature evolutions illustrate the common temperature response to the external forcing with the different climate sensitivities of both models, and also the internal variability of the global mean temperature at multi-decadal timescales, since the simulations within each ensemble display slightly different temperature paths.
For both ensembles of simulations of the past, the decadally smoothed time series of regional upwelling index extracted from each of the three simulations display different time evolutions (Fig. 6) and do not show a centennial evolution comparable to the global mean temperatures or the global mean external forcing as just described.The correlations between these series for each region are correspondingly low (Fig. 7) and mostly not statistically significant.These correlations remain low and non-significant for stronger low-pass time filtering up to 50 years (Fig. 8a).Based on Eq. ( 4), this indicates that the upwelling variance shared by all simulations, which could only be due to the common external forcing, is also small.These simulations, therefore, do not support any significant imprint of the external forcing on upwelling in any of the EBUSs up to multidecadal timescales over the past millennium.
This finding is consistent with previous studies derived from observed records, which identified a strong influence of the Pacific North American pattern on California upwelling.The Pacific North American pattern is a well-known mode of internal climate variability, unrelated to external forcing, and probably linked to the dynamics of the tropical and midlatitude Pacific Ocean (Macias et al., 2012).Thus, the influence of the Pacific North American pattern on the California Ocean Sci., 12, 807-823, 2016 At centennial timescales, the externally forced climate variability is a large portion of the total climate variability, as the random internal decadal variability is filtered out.Additionally, the increase in the external climate forcing over the past 156 years is stronger than the variations in the forcing prescribed in the past1000 simulations.For instance, the amplitude of decadal variations of the pre-industrial external forcing is of the order of 0.3 Wm −2 (Schmidt et al., 2011), whereas the increase in the external forcing over the past 250 years is 1.6 Wm −2 (Myhre et al., 2013).
Thus, the long-term influence of the external forcing -the strongest being anthropogenic greenhouse gases in the historical simulations -could be more easily detected in the form of common long-term upwelling trends.However, the simulated upwelling trends are small, mostly not statistically significant, and generally display opposite signs in the different simulations of the MPI-ESM ensemble (Fig. 6).The are only two cases out of eight, California and Canary in the historical simulations, in which all members exhibit a negative trend (opposite to the expected strengthening), of which only two trends of the California upwelling are statistically significant.Analysing the imprint of external forcing on upwelling with the CESM-LME confirms the results with the MPI-ESM.The simulated upwelling velocity in three simulations are weakly correlated with each other (Fig. 9) and the   trends are either not coherent in the simulations or are not statistically significant.

Scenarios
Looking at the future development of the upwelling in the EBUSs, we analyse the future simulations conducted with the MPI-ESM model under three different emission scenarios.An effect of the external forcing should be reflected in consistent centennial trends in upwelling in all members of each ensemble.Consistent significant trends in all three simulations of an ensemble only occur in the rcp8.5 scenario, the one with the strongest external forcing (Fig. 10).The trends are negative in California and Canary, and positive only in Benguela.In Peru no significant trends can been seen.This supports our results obtained for the past millennium, indicating that the external forcing has not been intense enough in the past to force the upwelling significantly.Even the much stronger external forcing of the rcp2.6 and rcp4.5 scenarios compared to the one of the past1000 and historical simulations does not leave an imprint on the upwelling.Furthermore, the envisaged positive upwelling trend only occurs in one of the four regions, underlining the local differences among the EBUSs.
These results do not fully support the ones obtained by Wang et al. (2015), agreeing on the positive trends in Benguela and the negative trend California.For Peru (Humboldt in their analysis) and Canary they found an intensification in upwelling.However, they analysed the trend for another period  and only in the rcp8.5 scenario.They concluded that a clear connection between external forcing and upwelling could be found.However, when analysing the trends simulated in other scenarios, it turns out that in the weaker greenhouse gas scenarios this link is weaker, and it cannot be easily identified.As Wang et al. (2015) used a model ensemble with one simulation of each model and we use ensembles of simulations of two models, differences could be due to the methodology.This would indicate that either the trend differ between the models or that using only one simulation of each model could distort the result.5 Importance of external climate forcing for the drivers of upwelling in the past1000 and historical simulations In this section, we want to determine why the imprint of the external forcing on coastal upwelling is weak over the past periods, by looking at the imprint of the external forcing on processes involved in the functioning of upwelling.These processes are, on the one hand, the atmospheric drivers of the upwelling and, one the other hand, the oceanic factors influencing the upwelling.

Atmospheric drivers
Upwelling in all four investigated regions is mainly driven by the wind stress curl that modulates upwelling in offshore regions and by the alongshore wind stress that drives the coastal upwelling.The upwelling in the four EBUSs is related to the SLP gradient between land and ocean (Tomczak and Godfrey, 2003).The strength of this gradient impacts the intensity of the trades (Bakun, 1990).We investigate whether the time evolution of this SLP gradient is consistent across the simulations and whether the lack of correlations between the upwelling indices across the ensemble may be due to www.ocean-sci.net/12/807/2016/Ocean Sci., 12, 807-823, 2016 a weak influence of the external forcing on wind stress and on the SLP gradient.The across-coastline SLP gradients in the simulations are calculated by subtracting the averaged air pressure over land from the averaged pressure over ocean in the regions identified as most closely correlated to the upwelling indices (Fig. 4).These areas are Benguela: 15 As in the case of the upwelling indices, the correlations between the time series of SLP gradients and alongshore wind stress across the simulations in each ensemble are all low and not significant (Tables 1 and 2).These correlations also remain low and not significant regardless of the time filtering (Fig. 8b).Correlating the time series of both drivers, the SLP difference between ocean and land and the alongshore wind stress, without low-pass filter shows their close relation with significant correlations for both periods, past1000 and historical, with correlation coefficients between 0.34 and 0.70.Regarding the longest timescales captured in the simulations, the long-term trends of the wind stress and of the SLP gradient within each ensemble have either inconsistent signs, or are not statistically significant, or are incompatible with the expected effect of the external forcing varying at centennial timescales.These expected trends are negative in past1000 due to the long-term orbital forcing and positive in historical due to the Bakun hypothesis (Bakun, 1990) based on the increase in greenhouse gas forcing.Therefore, even if the connection between the atmospheric drivers and upwelling might not be totally realistic in the Earth system model, the lack of common time evolution of wind stress or SLP across the simulations clearly shows that, from the at- mospheric perspective alone, the simulations do not support a discernible influence of the external forcing on the atmospheric drivers of upwelling over the last centuries.This implies that such an influence could not have been found even if the ocean model perfectly represented the real upwelling dynamics.Regardless of the quality of the upwelling representation in the coupled model, one could not detect any impact of external forcing on the upwelling.For the future period, significant trends of the alongshore wind stress in all three simulations of an ensemble occur in the rcp8.5 scenario with the same sign as for upwelling, negative for California and Canary and positive for Benguela.Additionally, there are significant negative trends in the whole ensemble in the rcp4.5 scenario for the California upwelling region.This indicates that, again, the external forcing has to be very strong to be discernible.

Stratification
Another possible mechanism by which the external climate forcing could influence upwelling involves the stratification of the water column.In periods with a stronger external forcing, the temperatures at the surface should warm more rapidly than in the deeper layers, increasing the stability of the water column, and hindering the mechanical effect of the alongshore wind stress (Hsueh and Kenney, 1972).The amount of variability in the SST that can be attributed to the variations in the external forcing can also be estimated by the correlation between the grid cell SST series simulated in each member of the simulation ensemble (Fig. 11).The correlation patterns indicate that the SST variability is more strongly driven by the external forcing in the tropical belt and tends to be weaker in the middle and high latitudes.This occurs despite the strongest response of high latitudes to external radiative forcing, known as the Arctic amplification (Taylor et al., 2013).The reason is that at high latitudes the internal variability is also larger than at low latitudes (Tett  , 2007).The ratio between both external forcing signal and internal variability, which is encapsulated by the correlation patterns shown in Fig. 11, is therefore highest in the tropics, a feature which has been so far mostly overlooked but that has been found in previous analysis of paleoclimate simulations (Tett et al., 2007).In the EBUSs, the correlations between the simulated SSTs are all positive, of the order of 0.2-0.3 after 10-year low-pass filtering (Fig. 11).These relatively low correlations indicate that the external forcing has some influence on the SST variability in these regions but most of the multi-decadal variability is internally generated.
The impact of the external forcing on the stratification in the future has not been analysed.There might be an impact due to much stronger increase in greenhouse gas concentrations in the 21st century, especially in the rcp8.5.

Discussion and conclusions
In this paper, the importance of the external climate forcing on the upwelling and its drivers in the four eastern boundary upwelling systems over the past millennium and the 21st century is investigated.Regarding Bakun's hypothesis (Bakun, 1990), the increase in radiative forcing linked to increased greenhouse gas concentrations would lead to an intensification of coastal upwelling.
The analysis of three simulation ensembles with the Earth system model MPI-ESM over the past millennium and the future are in contrast with the hypothesis of a discernible influence of the external forcing on coastal upwelling intensity.
Analysing ensembles of simulations enables us to distinguish between variations driven by internal and external variabilities.
For the periods covering the past, the past1000 and the historical simulations, no significant trends in all three simulations of an ensemble could be found in the upwelling itself, nor in its atmospheric drivers.Furthermore, correlations of the ensemble members are low, indicating the missing imprint of external forcing on the upwelling and its drivers.The internal variability dominates the external one for both periods in all four EBUSs.These results are in conflict with Bakun's hypothesis of an intensification of upwelling due to climate change (Bakun, 1990).The work of Narayan et al. (2010) shows that the results of a trend analysis of observed upwelling differs with the analysed period, data set, and analysed variable.
For the future, the effect of external forcing on the EBUSs can be identified when greenhouse gas concentrations are assumed to follow the strongest scenario, rcp8.5, among the three representative concentration pathways analysed here.All three simulations in the ensemble display consistent trends, but these trends are not always consistent with the expected intensification of upwelling predicted by Bakun (1990), with only Benguela showing an intensification, California and Canary showing a weakening, and Peru showing no significant trend.
Our results partly agree with the ones obtained by Wang et al. (2015) on the influence of a strongly increased future greenhouse gas forcing on upwelling intensity in the EBUSs, agreeing in an intensification of the Benguela upwelling system and a weakening of the California upwelling system.However, our results indicate that the conclusion obtained from the analysis of only the rcp8.5 scenario cannot be extended to weaker scenarios of future greenhouse gas forcing, nor to the trends observed over the 20th century, nor the evolution of upwelling over the past millennium.Differences in the results obtained from Wang et al. (2015) compared to our study could be due to several factors.The analysed time period in their analysis (1950-2099) differs from ours.They analysed only one simulation with each model and not an ensemble of simulations with the same model.Using an ensemble of simulations allows to better identify the presence or absence of an externally forced signal.Furthermore, longterm trends in their analysis (extended data, Fig. 3) show that there are positive and negative upwelling trends (blue and red squares in their figure) for all regions, so that clearly not all models agree in the sign of the upwelling trends.In addition, for some regions, e.g.Humboldt, this discrepancy in the trend may not only be due to the different model structure but also due to internal variability, as we find in the MPI-ESM model.
Compared to Rykaczewski et al. (2015), where again only the upwelling in the EBUSs in the rcp8.5 scenario in one simulation per model is analysed and there is one simulation per model, our results match theirs only for Benguela and California.They found significant positive trends in the upwelling-favourable winds in the Canary, Humboldt, and Benguela systems and negative trends in the California system, with trends in Benguela being less consistent among www.ocean-sci.net/12/807/2016/Ocean Sci., 12, 807-823, 2016 the models.This underlines the differences between the upwelling regions and illustrates the complexity of upwelling by the different results obtained depending on the analysed model, but also highlight the impact of the external forcing on upwelling detected in the rcp8.5 scenario.For the rcp8.5, there seems to be a difference between the EBUSs in the Northern Hemisphere (showing significant negative trends) and the Southern Hemisphere (showing significant positive or not significant trends).Differences could be caused by changes in the Hadley cell, which is expected to expand poleward as a consequence of global warming (Lu et al., 2007).The temperature differences between the equator and the subtropics will decrease in the Northern Hemisphere, due to the larger fraction of land on the Northern than on the Southern Hemisphere, leading to a weakening of the Hadley cell and thus to weaker trades (Ma and Xi, 2013).This effect may be smaller in the Southern Hemisphere due to a weaker warming of mid-latitudes and high latitudes (Ma and Xi, 2013).
A possible reason why upwelling does not show the signal of the variations of external forcing in the simulations, with the exception of the rcp8.5 scenario, may be that the atmospheric circulation has a large internal variability, in contrast to other variables, like temperature, which are more directly related to the external forcing.In general, detection and attribution studies of climate change detect weaker signal in atmospheric dynamical variables like SLP than in thermal variables.
Uncertainties still remain.For instance, the magnitude of the external forcing variations over the past millennium is still not well established (Schmidt et al., 2011), and larger variations than hitherto assumed may cause a tighter connection between forcing and upwelling than the one found in the simulations analysed here, driven by relatively weak variations of the external forcing.Over the past 156 years, however, the trends in external climate forcings are much more certain and over this period the historical simulations do not show any consistent sign of intensification or weakening in an ensemble.
It has to be kept in mind that our results are based on the realism of the analysed Earth system models.The relatively low model resolution of the atmosphere and of the ocean components could result in an unrealistic representation of the upwelling itself and/or its drivers.As stated by Small et al. (2015), especially the resolution of the atmospheric model may have the strongest influence on the simulated coastal upwelling.Furthermore, the current global coupled climate models still display a strong SST bias in the EBUSs.The cause of this bias is not completely understood, and it may be related to a deficient representation of coastal upwelling but it may also have other causes, for instance related to biases in the stratocumulus clouds (Richter, 2015).This caveat, nevertheless, also affects the recent studies by Rykaczewski et al. (2015) and Wang et al. (2015), since they are also based on the CMIP5 models.
The definition of the upwelling index can be sensitive to the region over which the vertical velocities are averaged.Redefining the subregions of the eastern boundary upwelling systems to the latitudinal extend of the regions used in Wang et al. (2015) (extended Canary region (18.5-10.5 • W, 16.5-42.5• N), south Benguela (8-30 • E, 28-40 • S), and Chile (80-70 • W, 20-40.5 • S)) does not change the main results of this study.The correlation between the three simulations remains low for upwelling.Significant trends with the same sign in all three simulations do not occur, neither in the simulations of past periods nor in the future simulations.
Analysing ensembles of simulation of the Earth system model CESM-LME over the past millennium supports the results of the MPI-ESM.
Thus, we conclude that the circumstantial evidence linking the recent observed trends in EBUSs upwelling to external climate forcing is in conflict with the present analysis of ensembles of simulations with two state-of-the-art climate models.

Data availability
All model data used in this study are publicly available.
The data of the simulations with the MPI-ESM-P (past1000, r1), the MPI-ESM-MR (historical), and the MPI-ESM-LR (future) are part of the Climate Model Intercomparison Project Version 5 (CMIP5).These data are openly available and can be downloaded from any of the nodes of the Earth System Grid Federation, after registration.A menu allows to select the type of simulation (past1000, historical or any of the RCP scenarios), the variable in question and the time resolution (e.g.daily or monthly).The data download can be performed interactively by clicking on the links of the individual files or by downloading a Unix script "wget", which can be locally run on a Unix computer.All files are written in "netcdf" format.Two additional MPI-ESM-P simulations (past1000 r2 and r3) were later conducted and were kindly provided by the Max Planck Institute for Meteorology (Zanchettin et al., 2012) (contact: johann.jungclaus@mpimet.mpg.de).
Detailed information about the Last Millennium Ensemble Project with the model CESM (Otto-Bliesner et al., 2016) can be found on line http://www2.cesm.ucar.edu/models/experiments/LME.The data are available under this link http://www2.cesm.ucar.edu/models/experiments/LME/data-sets.After the section of the desired variable by the user, the data need to be processed and made ready for download.The user is then electronically informed when the netcdf files are ready to be downloaded.discussions in the Cluster of Excellence Integrated Climate System Analysis and Prediction (CliSAP).The Max Planck Institute for Meteorology kindly provided the model data.There is no potential conflict of interest of the authors.We thank Dennis Bray for his editorial assistance.
The article processing charges for this open-access publication were covered by a Research Centre of the Helmholtz Association.
Edited by: M. Meier

Figure 1 .
Figure 1.Mean sea surface temperature in the four upwelling regions in the corresponding main upwelling season for 1985-2005 simulated by the Earth system model MPI-ESM compared to satellite data from AVHRR.California (June-August) AVHRR (a) and MPI-ESM historical r1 (b); Canary (June-August) AVHRR (c) and MPI-ESM historical r1 (d); Benguela (September-November) AVHRR (e) and MPI-ESM historical r1 (f); and Peru (June-August) AVHRR (g) and MPI-ESM historical r1 (h).Not the different temperature scale for Benguela and Peru.

Figure 2 .
Figure 2. Monthly mean vertical velocity at 52 m depth in four main coastal upwelling regions simulated in one past1000 (900-1849) (a) and one historical (1850-2005) simulation (r1) (b) of the MPI-ESM.The values for California have been re-scaled for better visibility.Note the different y axis scale for each simulation.As variations of the annual cycle are very small between the simulations of an ensemble, the annual cycle of r1 represents the simulated annual cycle of all simulations of the ensemble.

Figure 3 .
Figure 3. Correlation patterns between the seasonal (June-August, except for Benguela in September-November) upwelling field in each eastern boundary upwelling system and the simultaneous spatial mean alongshore wind stress (right panels) and the relative standard deviation of the upwelling field (left panels) simulated in one of the past1000 simulations (r1) (900-1849) with the Earth system model MPI-ESM.

4
Importance of external climate forcing for upwelling variations in the recent past and future 4.1 Past1000 and historical simulations

Figure 4 .
Figure 4. Correlation patterns between the seasonal (June-August, except for Benguela in September-November) upwelling index in each eastern boundary upwelling system and the simultaneous seasonal mean sea level pressure field simulated in one of the past1000 simulations (r1) (900-1849) with the Earth system model MPI-ESM.

Figure 5 .
Figure 5. Global mean near-surface air temperature over the past millennium simulated in two ensembles of simulations with the Earth system models MPI-ESM and CESM-LME.The time series represent 31-year running mean of the temperature anomalies relative to the 20th century mean.

Figure 6 .
Figure 6.Time series of the simulated upwelling indices in each upwelling region simulated in two ensembles of climate simulations (three members each), denoted past1000 (900-1849) and historical (1850-2005), with the MPI-ESM model.The plus and minus signs in each panel indicate the sign of the long-term trend, their value being included whenever statistically significant.The series have been low-pass filtered, with a 30-year filter for past1000 and a 10year filter for the historical ensembles.

Figure 7 .
Figure 7. Frequency histogram in bins of 0.05 width showing the distribution of across-ensemble correlations between the upwelling indices simulated in each upwelling region for the ensembles past1000 (a) and historical (b) of the MPI-ESM, after low-pass filtering with a 30-year (past1000 simulation) and a 10-year (historical simulation) filter, respectively.The vertical black lines indicate the 5-95 % significance bounds taking the time filtering into account.

Figure 8 .
Figure8.Correlations between the simulated upwelling indices in all four EBUSs in three past1000 simulations with the MPI-ESM model after filtering the time series with low-pass filter of increasing period (a).The same as (a) but for the time series of sea level pressure difference between land and ocean averaged over the areas most closely correlated with the simulated upwelling in this region (b).The 95 significance levels depend on which pair of series are being tested, since the test takes into account the serial correlation of each series.The figure shows, as an example, the significance levels for the correlations r1-r2 for Benguela (bold black line), but all correlations shown here are not statistically significant at the 95 % level.

Figure 9 .
Figure 9. Frequency histogram in bins of 0.05 width showing the distribution of across-ensemble correlations between the upwelling indices simulated in each upwelling region for the ensembles past1000 (a) and historical (b) of the CESM-LME, after lowpass filtering with a 30-year (past1000 simulation) and a 10-year (historical simulation) filter, respectively.The vertical black lines indicate the 5-95 % significance bounds taking the time filtering into account.

Figure 10 .
Figure 10.Time series of the simulated upwelling indices in each upwelling region simulated in three ensembles of climate simulations (three members each) for 2006-2100, rcp2.6,rcp4.5, and rcp8.5, with the MPI-ESM model.The plus and minus signs in each panel indicate the sign of the long-term trend, their value being included whenever statistically significant.The series have been low-pass filtered with a 10-year filter.

Figure 11 .
Figure 11.Correlation pattern between the global skin temperatures (temperatures of ocean and land surface) simulated in two past1000 simulations (r1 and r2) (900-1849) with the Earth system model MPI-ESM in the June-August season after applying a 10-year lowpass filter.

Table 1 .
Cross-correlation coefficients of the three simulations of each ensemble (r1, r2, and r3) of the cross-coastline pressure gradient of all upwelling regions for the past1000 and historical simulations with the MPI-ESM model after 30-and 10-year filters, respectively.All correlations are not statistically significant at the 95 % level.

Table 2 .
Cross-correlation coefficients of the three simulations of each ensemble (r1, r2, and r3) of the alongshore wind stress of all upwelling regions for the past1000 and historical simulations with the MPI-ESM model after 30-and 10-year filters, respectively.All correlations are not statistically significant at the 95 % level.