Decadal variability and trends of the Benguela upwelling system as simulated in a high-resolution ocean simulation

Detecting the atmospheric drivers of the Benguela upwelling systems is essential to understand its present variability and its past and future changes. We present a statistical analysis of a high-resolution (0.1) ocean-only simulation driven by observed atmospheric fields over the last 60 years with the aim of identifying the large-scale atmospheric drivers of upwelling variability and trends. The simulation is found to reproduce well the seasonal cycle of upwelling intensity, with a maximum in the June–August season in North Benguela and in the December–February season in South Benguela. The statistical analysis of the interannual variability of upwelling focuses on its relationship to atmospheric variables (sea level pressure, 10 m wind, wind stress). The relationship between upwelling and the atmospheric variables differ somewhat in the two regions, but generally the correlation patterns reflect the common atmospheric pattern favouring upwelling: southerly wind/wind stress, strong subtropical anticyclone, and an ocean–land sea level pressure gradient. In addition, the statistical link between upwelling and large-scale climate variability modes was analysed. The El Niño–Southern Oscillation and the Antarctic Oscillation exert some influence on austral summer upwelling velocities in South Benguela. The decadal evolution and the longterm trends of simulated upwelling and of ocean-minus-land air pressure gradient do not agree with Bakun’s hypothesis that anthropogenic climate change should generally intensify coastal upwelling.


Introduction
The Benguela upwelling systems (BUS; North Benguela and South Benguela) is one of the four major Eastern Boundary Upwelling Systems (EBUs) of the world (Shannon, 1985;Leduc et al., 2010) and among the most productive oceanic regions (Bakun et al., 2010;Leduc et al., 2010). Nutrient-rich coastal upwelling in the EBUs are mainly driven by wind patterns that cause offshore Ekman transport that cannot be balanced by the horizontal advection of water (Bakun and Weeks, 2004). Further offshore, the wind stress curl causes upwelling by Ekman pumping. These areas are thus rich in pelagic fish biomass and important for coastal fisheries (Bakun et al., 2010). The BUS, off Angola, Namibia and South Africa (Blanke et al., 2005), has its northern boundary at the Angola-Benguela Front (between 14 and 17 • S) (Shannon and Nelson, 1996;Veitch et al., 2010). The southern boundary of the BUS is defined by the Agulhas retroflection at the southern tip of the continent where Agulhas Current water penetrates into the southern Benguela at 37 • S (Shannon and Nelson, 1996). The main driver of the BUS is thought to be the wind stress off southwestern Africa (Nelson and Hutchings, 1983) while north of the front the upwelling is more strongly related to the larger-scale equatorial upwelling (Hardman-Mountford et al., 2003). The strongest upwelling takes place near Lüderitz (27 • S) (Shannon and Nelson, 1996), resulting in intense cold sea surface temperatures (SSTs) that persist throughout the year (Parrish et al., 1983). This upwelling cell naturally divides the BUS into a northern and a southern part (Shannon and Nelson, 1996;Hutchings et al., 2009). In the southern Benguela the upwelled water occurs near the coast, whereas in the northern Benguela it ex-tends farther off the coast, westward up to about 150-250 km (Shannon and Nelson, 1996;Fennel et al., 2012).
The upwelling intensity is highly seasonal in the temperate latitudes (30-34 • S), with generally higher intensity in boreal spring and summer, but more uniformly distributed across seasons in regions closer to the Equator (15-30 • S) (Parrish et al., 1983;Mackas et al., 2006;Chaigneau et al., 2009;Chavez and Messié, 2009). However, a clear picture of the upwelling seasonality is not established yet. Bakun and Nelson (1991) stress that the strongest upwelling takes place in the austral summer seasons, Blanke et al. (2005) widen this seasonal frame from October to March. In contrast, Hagen et al. (2001) and Hagen et al. (2005) argue that the main upwelling season is rather the austral winter to spring. Veitch et al. (2010) clarify that upwelling peaks in the southern Benguela in austral spring and summer, whereas in the north, with weaker seasonal variations, it does in austral autumn and spring. Central Benguela upwelling near Lüderitz does not a show clear seasonal cycle (Shannon and Nelson, 1996). This distinction supports the division in two distinct regimes (Shannon, 1985) and highlights the complex nature of the BUS, which warrants a closer look at the large-scale climate drivers.
The Benguelan upwelling is thought to be driven by several factors. The coastal topography and the climatological winds frame the areas of upwelling (Shannon, 1985;Chavez and Messié, 2009). The region is dominated by southerly and southeasterly winds (Hagen et al., 2001;Risien et al., 2004) that are influenced by the high pressure system over the South Atlantic, by the cyclones moving westward over the southern Benguela and by the pressure over southern Africa (Shannon and Nelson, 1996;Risien et al., 2004). Seasonal trade winds also influence the dynamics of the upwelling (Shannon, 1985;Chavez and Messié, 2009). In addition, the BUS displays some particular characteristics, compared to other EBUs, determined by its geophysical boundaries: the Agulhas retroflection, the Angola-Benguela Front and the passage of westerly winds in the south (Shannon and Nelson, 1996). Zooming more closely into the BUS, Lachkar and Gruber (2012) characterised the Benguela subregions as having a shallow mixed layer depth, but with central Benguela dominated by strong upwelling (net primary production) due to a wide shelf and low eddy activity, whereas southern Benguela is characterised by weaker upwelling, a wide shelf, and moderate to high eddy activity. In the northernmost Benguela, the upwelling is moderate in a narrow shelf and moderate to high eddy activity.
The upwelling intensity is not constant, but presents intraseasonal variations. In general, the variability of the Benguela system is dominated by the intraseasonal (eddy) variability (Chavez and Messié, 2009) and coastal-trapped waves (Rouault et al., 2007). Lachkar and Gruber (2012) identified factors that inhibit the net primary production and thus have indirectly been linked to upwelling dynamics, namely strong eddy activity, a narrow continental shelf and a deep mixed layer. These factors help to provide a spatial characterisation of the Benguela system and its spatial heterogeneity (Lachkar and Gruber, 2012).
From a larger-scale climatic perspective, climate modes seem to have an impact on the interannual variability of upwelling through modulation of the local conditions. Hagen et al. (2001) mentioned a possible influence of the Quasi-Biennial Oscillation (QBO), a mode of variability that describes quasi-oscillations of low stratospheric winds. The St. Helena Island Climate Index(HIX, first mode of the empirical orthogonal function (EOF) analysis of air temperature, sea level pressure (SLP) and precipitation) is positively correlated with the southeasterly trades and thus should modulate the strength of the upwelling (Hagen et al., 2005). Furthermore, the Antarctic Oscillation (AAO), which is related to the strength of the circumpolar westerly flow (Jones and Widmann, 2004), could also influence the upwelling intensity.
In addition to the mentioned climate modes, the impact of El Niño-Southern Oscillation (ENSO) on the Benguela upwelling could, according to some authors (e.g. Dufois and Rouault, 2012), be of major importance, even stronger than the one of the AAO (Rouault et al., 2010). According to this study, El Niño and La Niña events tend to weaken and strengthen upwelling, respectively. The proposed physical mechanism is related to an equatorward shift in the high pressure over the southern Atlantic which leads to weaker upwelling-favourable trades (and conversely for La Niña events) (Dufois and Rouault, 2012). The results obtained by Shannon and Nelson (1996) support this by identifying the relation to ENSO as the most significant one on interannual timescales. Deutsch et al. (2011) also presented changes in the upwelling intensity and the depth of the thermocline associated with ENSO. Some authors (e.g. Huang et al., 2005) suggest that ENSO, affecting the whole tropical belt around the Earth, also has an influence on the so-called Benguela Niños. These extreme events are a type of Atlantic Niños that reach the BUS when warm waters extend into the eastern Atlantic around 600 km further south than usually (Shannon et al., 1986). The Benguela Niño events have been explained by two different causal chains. According to the first, SST anomalies along the southwest African coast develop due to zonal wind anomalies in equatorial South America and the equatorial western Atlantic (McCartney, 1982;Wang, 2006). Kelvin waves crossing the basin are induced due to the relaxation of the trades along the Equator (Rouault et al., 2007). The unusual sea surface height reaches the Angola-Benguela Front which leads to increased intrusion of tropical water into the Angola-Benguela upwelling system (McCartney, 1982;Rouault et al., 2007). According to a second scenario (Richter et al., 2010), the main drivers of Benguela Niños are alongshore wind anomalies, which occur due to the weakening of the South Atlantic High that develops 2-3 months in advance of Benguela Niños and 5-6 months before the  (Rouault et al., 2010;Dufois and Rouault, 2012). The warming events have been found to be weaker and less frequent than the Pacific ENSO (Shannon et al., 1986;Latif and Grötzner, 2000). In contrast, Hardman-Mountford et al.
(2003) estimate a much higher frequency of warm events in the Atlantic than in the Pacific, but these authors found that the warming events penetrate more rarely into the Benguela system. For this to occur, southerly winds at the Equator are required to blow simultaneously with the southward water flow (Hardman-Mountford et al., 2003). Shannon et al. (1986) find that, in any case, the effect on the southern Benguela is small, whereas in the northern Benguela the intrusions of warm water persist for at least 6 months, leading there to extreme events in precipitation in Namibia and Angola (Rathmann, 2008), and a weaker or shorter (about 2 months) upwelling season (Hagen et al., 2001). Hagen et al. (2001) found a periodicity of the occurrence of Benguela Niños of around every 11 years  and between 1974 and 1999 every 5 years, and interpreted this shift as a tendency of these extreme events to occur more often. Another ongoing discussion hints to a more complex interaction from the tropical Pacific into the tropical Atlantic as a simple unidirectional influence from the Pacific to the Atlantic. Rodríguez-Fonseca et al. (2009) mentioned that ENSO events are preceded by Atlantic events of opposite sign and that the Atlantic could strengthen ENSO. Ham et al. (2013) agree and explain that easterlies over the equatorial far-eastern Pacific are considerably weaker under Atlantic Niña conditions. The same study adds that the northern tropical Atlantic is more strongly connected with the Pacific ENSO than with the Atlantic Niño. In contrast, Wang (2006) detected no correlation between both El Niños (Pacific and Atlantic) but found that anomalous SSTs of the two equatorial oceans can induce an inter-Pacific-Atlantic SST gradient that leads to anomalies in surface zonal wind over parts of both oceans and equatorial South America. These anomalies of the zonal wind constitute the bridge conveying the interaction between the two oceanic basins. Wang (2006) stated that the SST gradient between the equatorial Pacific and Atlantic is of higher importance than the individual ocean SST anomalies as driver of the atmospheric circulation across northern South America. Contrasting to the results of Rodríguez-Fonseca et al. (2009), Latif andGrötzner (2000) mentioned that the equatorial Atlantic responds to ENSO with a lag of 6 months, this time being needed by the equatorial Atlantic to adjust to low-frequency wind stress variations. Also Colberg and Reason (2006) found that ENSOinduced wind anomalies play a major role in driving upper South Atlantic temperatures by an atmospheric teleconnec-tion with a lag of one season. Thus, even if the direction of teleconnection is unclear, the influence of ENSO and the Benguela Niños on the Benguelan upwelling could be of major importance when investigating the atmospheric drivers due to their effects on the Benguelan ecosystem.
The long-term trend of the Benguela upwelling brought about by anthropogenic climate change is another important reason to identify its large-scale atmospheric drivers. In his landmark paper, Bakun (1990) put forward his hypothesis that can be summarised as follows: as a consequence of rising greenhouse gas concentrations, the surface temperature over the continents warm faster than oceans. This leads to a strengthening of the continental lows and oceanic highs. The land-sea pressure gradient is thereby enhanced, causing a strengthening of alongshore winds and enhancing upwelling. In addition, the warming ocean surface results in more humidity in the otherwise very dry atmosphere over land, which reinforces the greenhouse gas effect (Bakun, 1990). Observations do indicate that over southwestern Africa surface temperatures over land have increased more rapidly than over the ocean from 1980 onwards, although in previous periods the statistical significance of the trend is compromised due to a poorer data quality (Hartmann et al., 2013, Fig. 2.22). In agreement to Bakun's hypothesis, Narayan et al. (2010) found decreasing trends in coastal SSTs in four major EBUs over 1960-2000 and attributed them to meridional wind stress intensification. This was supported by Demarq (2009), who found a positive trend in southerly winds in the Benguela EBU for 2000-2007. In addition, Rathmann (2008) reported an intensification of the trade winds due to SLP trends. However, Narayan et al. (2010) recognised the spread in estimation of the wind stress trend using different data sets (reanalyses and observations) and Belmadani et al. (2014) could not identify the mechanisms linking the land-sea thermal contrast and upwellingfavourable winds in the Peru upwelling system in future climate model simulations. Bakun et al. (2010) also expressed some doubts about the robustness of the estimated upwelling trends in Benguela. To explain this apparent lack of clear upwelling trends in BUS, Bakun et al. (2010) later reasoned that the influence of ENSO on atmospheric humidity in Benguela, combined with the succession of recent strong ENSO, events may have stifled a long-term intensification of upwelling that should have occurred as a result of the anthropogenic greenhouse effect. However, climate models still give an unclear picture about the long-term effect of anthropogenic greenhouse gas forcing on ENSO (Vecchi and Wittenberg, 2010), but at least 9 of 22 models of the Coupled Model Intercomparison Project phase 5 (CMIP5) show an increasing trend in ENSO amplitude before 2040 and a decreasing trend thereafter (Kim et al., 2014). In addition, global climate models still display clear systematic errors in simulating the SST in the eastern ocean basins (Richter and Xie, 2008;Echevin et al., 2012;Grodsky et al., 2012), the origin of these errors still being under investigation.  Against the backdrop of these competing hypotheses about the role of large-scale climate modes on Benguela upwelling and the long-term evolution of upwelling itself, the purpose of this paper is to analyse the atmospheric drivers of the interannual variability and the long-term trends of upwelling in the BUS as simulated in a high-resolution (about 0.1 • long-lat) ocean simulation with the ocean model MPI-OM developed at the Max Planck Institute for Meteorology in Hamburg, and spanning the last 6 decades driven by the NCEP/NCAR (National Centers for Environmental Prediction/National Center for Atmospheric Research) meteorological reanalysis.
The focus of this analysis lies on the longest timescales possible for such high-resolution simulations. This focus requires a trade-off between a possibly higher quality atmospheric forcing data provided by most recent meteorological reanalysis, on the one hand, and the period covered by those data sets, on the other hand, which is considerably shorter than those of the NCEP/NCAR reanalysis. Although the ocean simulation was driven by the NCEP/NCAR reanalysis, some analyses were confirmed with the more recent and higher-resolution ERA-Interim reanalysis. However, through this study, it has to be borne in mind that the conclusions obtained may be dependent on the realism and homogeneity of the NCEP/NCAR product. The resolution of the wind stress provided by this meteorological reanalysis may be too coarse to fully represent the effect of the atmospheric forcing on the ocean dynamics. The results of this simulation must, therefore, be considered as conditional on the realism of this wind forcing and should be contrasted against future simulations that will eventually use higher-resolution atmospheric forcing. For the time being, however, due to the high computational costs required for a high-resolution atmosphere, other simulations to estimate the recent and future evolution of upwelling have also used either meteorological reanalysis or atmospheric models of comparable horizontal resolution (Serazin et al., 2015;Wang et al., 2015). It has to be kept in mind that the results of this analysis are based on a simulation and that it is very difficult to establish the relia-bility of simulations when the direct observations are fragmentary and other gridded observation data sets are themselves based on models. Therefore, this study does not claim to represent the truth, but rather our goal is to use a state-ofthe-art model simulation and compare the conclusions that can be derived from it with predictions about long-term environmental changes based in more general principles like Bakun's hypothesis. This approach should contribute to improve our understanding of those changes and also of possible model deficiencies in situations in which long-term observations are incomplete.
The paper has the following structure. The data and methods are described in Sect. 2, followed by a basic validation of the model output and the NCEP/NCAR reanalysis regarding the upwelling. The results are presented in the following two sections. First, the analysis of the upwelling index, its trends and variabilities, and correlations with drivers are shown. Second, the results of the investigation on the SLP gradient are presented. Finally, these results are discussed in Sect. 6, with the most important conclusions in Sect. 7.

Data and methods
For this study we analyse climate variables derived from observational data sets, atmospheric data from meteorological reanalyses and from a high-resolution global ocean-only simulation. A summary of the characteristics of these data sets is presented in Table 1.

Gridded observational data
Atmospheric variables (10 m-wind, wind stress, and SLP) were obtained from the NCEP/NCAR reanalysis 1. This data set is available from 1948 onwards with a horizontal resolution of 2.5 • (Kalnay et al., 1996). Monthly means were downloaded and averaged to seasonal means of the area 50 • W-40 • E, 0-50 • S. An analysis of the gradient of the SLP field was performed with the SLP of the ERA-Interim reanalysis (available from 1979 onwards) (Dee et al., 2011) as well as of the NCEP/NCAR reanalysis 1. Here we chose In addition to the reanalyses data, the observational gridded sea surface temperature data set HadISST1 with a grid spacing of 1 • and data between 1870 and 2012 (Rayner et al., 2003) was used in this study for the validation of the model output. For the same purpose, 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, was used to detect the possible SST bias of the model data described in the following. However, it is known that this version of pathfinder has a warm bias of up to 5 • C close to the coast in austral summer in the southern Benguela .

Model data
The ocean simulation analysed here is a global simulation (hereafter denoted STORM) with a high-resolution version of the state-of-the-art model MPI-OM (Marsland et al., 2003) covering the period 1950-2010 driven by the global NCEP/NCAR meteorological reanalysis (von Storch et al., 2012;Li and von Storch, 2013). The horizontal resolution of this model version is 0.1 • and it has 80 levels. The model is driven in this simulation by fluxes of heat, fresh water and momentum derived from the meteorological reanalysis. Note that in this stand-alone simulation, there is no feedback of the simulated state of the ocean onto the driving fluxes. In par-ticular, the prescribed wind stress flux does not consider the effect of the simulated ocean currents.

Upwelling indices
To identify the large-scale atmospheric drivers, an upwelling index had first to be defined. The definition may depend on the available data, since not always are direct data of the vertical water mass transport readily available. Chen et al. (2012) defined an index derived from the SST and mentioned that it represents well the upwelling intensity in terms of the spatial variation while the index defined derived from the offshore Ekman transport represents the temporal variations. The same study found that the index constructed from the SST is not able to correctly provide the upwelling intensity during warm water intrusions of the Angola Current.
In general, SSTs are not only affected by upwelling but also by a complex interaction between horizontal ocean advection, ocean-atmosphere heat fluxes and vertical mixing (McCabe et al., 2015). Therefore, the simulated vertical mass transport at 52 m depth, close to the modelled mixed layer depth in the Benguela region, is used as upwelling index. All data sets were seasonally averaged using the standard seasons definition: December-February (DJF), March-May (MAM), June-August (JJA) and September-November (SON). The index derived from the vertical velocity is denoted here as up_wvelo. Previous studies (Shannon and Nelson, 1996;Hutchings et al., 2009) suggested that there exists a strong difference in the behaviour of the upwelling in the northern and southern Benguela. Thus, the area was divided in two regions with the border at Lüderitz (28 • S) (Fig. 2a , b). The South Benguela region covers around 2.7 times the oceanic region of North Benguela, because the south coast of South Africa with the Indian Ocean is included too. There, upwelling takes place as well (Shannon and O'Toole, 2003). Thus, the region denoted by South Benguela contains the southern BUS and additionally the south coast upwelling. Our results are not sensitive regarding the chosen upwelling region. A southern Benguela region restricted close to the coast does not change the results of this paper.

Climate indices
The upwelling indices were not only correlated with the atmospheric variables but also with climate indices. The influence of the Multivariate ENSO index (MEI; Wolter and Timlin, 1993), the Antarctic Oscillation (AAO; Marshall, 2003), and the tropical Atlantic (ATL-3, 20  off southern Africa were examined. The ATL-3 values were calculated with the data sets of HadISST1 and STORM. For the analysis of the statistical significance of the longterm trends in the upwelling indices and the linear correlations between the indices, a two-sided significance level of p = 0.05 was adopted. The p value of the correlation and of the trend significance indicates the significance of the correlation (or trend) under the usual assumption applied in ordinary least-squares regression of serially uncorrelated and normally distributed regression residuals. The series of upwelling in North Benguela clearly display decadal variability, so that the test of significance of a linear trend is not correct using the usual assumptions of uncorrelated regression residuals. This is why we additionally used another, more so-phisticated, method based on Monte Carlo generation of surrogate series by resampling the original series and on phase randomisation to preserve the autocorrelation structure of the original series (Ebisuzaki, 1997). Using this method, the estimated linear trend in North Benguela is no longer statistically significant.
3 Long-term means, annual cycle and link between upwelling and sea surface temperature Figure 2 shows the long-term annual means of some important variables obtained from the STORM simulation in the South Atlantic and in the southeast Atlantic: the ver- tical velocity at 52 m depth, SST (compared to the corresponding fields derived from HadISST1) and surface currents. Figure 2a and b display the upwelling cells at the western and southern coast of southern Africa. The SST pattern of STORM (Fig. 2c, d) and HadISST1 (Fig. 2e, f) look very similar, with colder SST at the coast due to upwelling. The SST of the STORM simulation shows lower values than the observational data set HadISST1. This may be due to relatively coarse resolution of the HadISST1. In Fig. 2g and h one can see the realistic representation of the equatorial currents, the Agulhas Current, Benguela Current and the westward drift along the west coast caused by the trade winds. The simulated SST in STORM compared to that of AVHRR is shown in Fig. 3. To have the same spatial resolution in both data sets, the AVHRR data set was spatially smoothed by averaging two grid points in both spatial directions. If the two grid-points include one missing value, the missing value is ignored in the mean. The STORM data set was interpolated onto the grid of AVHRR. The STORM simulation displays a warm bias in the upwelling region, as it is well-known for general circulation models. The issue of SST bias in climate models is being vigorously investigated, since it is pervasive across many climate models and may be related to the global climate sensitivity as well. The reason is not clear, as it has been recently reviewed by Richter (2015). Both atmospheric and oceanic processes may be involved: stratiform clouds in the boundary layer, too weak winds, and remote oceanic influence from the tropics. The cold bias located along the coast may be also partially caused by the bias of the AVHRR data mentioned by . The stronger bias of AVHRR underlines the more realistic cold SSTs simulated by STORM along the coast compared to HadISST1.
A presentation of NCEP is displayed in Fig. 4. The climatology of the alongshore wind stress, the wind stress curl and SLP show a realistic pattern for the South Atlantic Ocean and the Benguela region. Nevertheless, cyclonic wind stress curl takes place in the wind drop-off zone near the coast. This contribution to the coastal upwelling might be underrepresented in the NCEP/NCAR reanalysis data set due its coarse resolution (Chelton et al., 2004). The long-term trends of NCEP winds over the whole simulation period are small and not significant (Fig. 5). The validation of these trends is more problematic due to the lack of in situ observations in this area. Comparing NCEP reanalysis to the World Ocean Circulation Experiment (WOCE) shows an underestimation of the pressure systems in our regions and thus of the nearsurface winds (Smith et al., 2001). The comparison of NCEP reanalysis to QuickSCAT displays the deficits of NCEP in capturing small-scale features of the wind stress (Risien and Chelton, 2008), although it is difficult to assess whether these limitations also influence the long-term trends in the wind stress The results of the STORM simulation are admittedly conditional on the quality of the NCEP forcing, but unfortunately there are no clear alternatives for the study of the longterm variability of the atmospheric forcing. Satellite data as well as more recent reanalyses with higher resolution surely represent an improvement compared to NCEP reanalysis but the longer temporal coverage is a key element for our analysis, and this is missing in the higher-resolution data sets.
The annual cycle of the upwelling index in North and South Benguela is displayed in Fig. 6. The up_wvelo index in South and North Benguela display clear differences.
Whereas in the north upwelling shows a broad maximum in JJA, in South Benguela it tends to be stronger in DJF, due to the position of the subtropical Atlantic high (Mackas et al., 2006). This is in broad agreement with observations, considering the discrepancies already indicated in the introduction (Sect. 1).
The relationships between the upwelling index up_wvelo and the SST are displayed in Fig. 7 as a map of the correlation coefficient between the upwelling index and the SST in each model grid cell. These correlations are negative near the coast and more pronounced in North Benguela. The correlations between the up_wvelo of South Benguela and the SST field in STORM are negative along the entire coast in DJF, with stronger correlations in the south. In JJA the correlations are negative only at the southern coast of South Africa.   ranges) expressed as the change of upwelling intensity over the period 1950-2010 relative to their long-term mean. The long-term mean values themselves may depend on the depth, area, and quite presumably on model resolution, but they are in the reasonable range of tens of metres per year, which agrees with estimations using simplified models of upwelling dynamics (Rykaczewski and Checkley, 2008). The time evolution of the upwelling indices (deviations from their long-term mean) through the period of the STORM simulation (Fig. 8) indicates that the variations of upwelling in the north and in the south do not go hand in hand. The significant linear correlation between up_wvelo in North and South Benguela is 0.30 in DJF and −0.52 in JJA, indicating that the interannual variations of upwelling during JJA (the season with maximum upwelling in North Benguela) is out of phase in both regions. These results strongly suggest that the Benguelan upwelling cannot be considered as a homogeneous system, in accordance with previous observational studies (see Sect. 1).
The linear trends of the upwelling indices in the STORM simulation are positive in North Benguela in MAM, JJA and DJF and positive in South Benguela in MAM, SON and DJF (Table 2). Upwelling in North Benguela is about twice as intense as in South Benguela (Figs. 6, 8) and it is less variable relative to its long-term mean (   1950 1960 1970 1980 1990 2000 2010 1950 1960 1970 1980 1990 2000 2010 year the north than in the south ( Table 2). The second characteristic is that the decadal variations are much more coherent across seasons in the north than in the south (Fig. 8). There is obviously some persistent factor in the tropical upwelling region that coherently drives the upwelling decadal variability in all seasons. In the south, in contrast, the seasonal upwelling indices display a more incoherent behaviour across the four seasons. In addition, the timescales of variations are longer in the north than in the south, as it will be more clearly shown later with a spectral analysis of these indices.
The analysis of the significance of the trends shows that only in North Benguela the positive trends in MAM and DJF are nominally significant (i.e. assuming Gaussian and serially uncorrelated trend residuals; see confidence intervals in Table 2). However, Fig. 8 clearly shows that a linear trend in the north is not a good description of the long-term evolution of upwelling. In general, upwelling intensity remains stable from 1950 to 1985 and undergoes a decadal surge and subsequent slowdown until the end of the simulation. This behaviour is not what one would expect from a steady increase in external climate forcing over the second half of the 20th century. These long-term changes over the whole simulation period are of the same order as the coefficients of variation (standard deviation divided by the long-term mean), if not weaker, which further confirms that the decadal variations are dominating the temporal variability. A more sophisticated analysis of the statistical significance considering the serial correlation of the trend residuals (see Sect. 2) indicates that both trends are not significant at the 95 % level. In the south, the trend in DJF is close to the p = 0.05 significance level, but the long-term evolution, as displayed in Fig. 8, does show a steady, albeit weak, increase through the simulation. The trend residuals are in this case not serially correlated. The central estimate of the trend in South Benguela upwelling in DJF implies a linear increase of about 40 % from the start through the end of the simulation. We will return to this point when discussing the connection of upwelling to the largescale climate modes.

Correlations with large-scale climate fields
Correlations with atmospheric fields were calculated for the up_wvelo index. Figure 9 depicts the correlation patterns between the upwelling index and wind stress (defined as pos-itive when directed into the ocean) for DJF and JJA. This figure not only illustrates the expected relationship between upwelling and wind stress but also the seasonality in this relationship. In North Benguela, wind stress variability is more important for upwelling during JJA -the season with stronger upwelling -whereas in South Benguela the seasonality is not very marked. In both regions, the direction of upwelling-favourable winds does not change with season. These wind stress correlation patterns can also explain the mutual correlations between the north and south upwelling indices. As indicated before, the correlation between the interannual variations in upwelling in both regions in JJA is negative (r = −0.52). North Benguela upwelling is statistically connected to winds south of the southern tip of Africa that actually inhibit upwelling in South Benguela, and vice versa. In DJF, by contrast, this does not happen, and North Benguela upwelling is statistically linked to winds south of South Africa that are also upwelling-favourable for South Benguela, and vice versa, explaining the weak but positive correlation between both upwelling indices in this season. In the STORM simulation, upwelling is strongly correlated with upwelling-favourable wind stress in North Benguela. South Benguela upwelling in DJF is strongly driven by wind stress at the southern coast and not so pronouncedly by the wind stress at the southern west coast. This could be due to our definition of the region chosen as South Benguela, which includes parts of the coast of the Indian Ocean as well. Figure 10 shows the corresponding correlation patterns to the wind stress curl. The patterns also display the expected negative correlation between upwelling in the boxes defining North and South Benguela and the wind stress curl, as the Coriolis parameter is negative. The exception here is North Benguela in the DJF season, when upwelling is weakest. However, further away from the coast, the correlation with the wind stress curl changes sign for both regions and both seasons, indicating that near-coastal upwelling driven by the wind stress curl and upwelling further offshore in the ocean would include anticorrelated contributions. Widening the geographical boxes in which regional upwelling is defined to include more open-ocean regions would thus lead to smaller upwelling variability than in the coastal and open-ocean regions considered individually.
Summarising, the up_wvelo in North Benguela is driven mainly by the expected wind pattern. In contrast, the correlations with South Benguela show upwelling-favourable wind stress that is more pronounced at the south coast than at the west coast south of Lüderitz. The correlation pattern with the SLP supports this with a positive pressure anomaly south of the continent instead of an intensified subtropical high located in the central South Atlantic (not shown here).

Correlations with climate modes
To identify the connection between climate modes and the upwelling index up_wvelo, the correlations of this index with ENSO (MEI), the AAO, the ATL-3, HIX, the QBO, and the AMM were also calculated ( Table 3). The upwelling in North Benguela in STORM is significantly positively correlated with MEI in MAM, with the AAO in MAM and DJF and HIX in DJF. Significant negative correlations were found with AAO in JJA and ATL-3 in DJF (Table 3). In contrast, the upwelling in South Benguela is significantly negatively correlated with ENSO in MAM and DJF and positively with AAO in all four seasons and AMM in DJF (Table 3). The correlations with the QBO are all insignificant.
The correlation pattern between wind stress and the ENSO index in DJF (Fig. 11) clearly explains the stronger negative correlations with South Benguela upwelling and the negligible correlation in the north. The wind stress anomalies associated with ENSO tend to weaken the climatological winds in the south and are very small in the north. The pattern of wind stress anomalies related to ENSO have a large-scale character, spanning almost the whole South Atlantic at midlatitudes, being clearly strongest in South Benguela. The correlation between ENSO and wind stress in North Benguela is quite weak, consistent with the weak correlation found between ENSO and North Benguela upwelling. The wind stress correlation pattern indicates that, if anything, ENSO may contribute to drive upwelling anomalies in South Benguela but not in the north. Even if the correlation is weak, the wind stress pattern favoured by ENSO shows the upwellingunfavourable influence of a positive phase of ENSO. We suggest that these correlation patterns, which display systematically large values at midlatitudes and very small values in the tropical regions, arise as a result of a large-scale dynamical teleconnection with ENSO via the midlatitudes, rather than due to a tropical bridge between the tropical Pacific and the tropical Atlantic, as can be seen in Fig. 11. Also, a link between ENSO and Benguela upwelling resulting from the modulation of smaller-scale humidity dynamics in the Benguela region, as discussed by Bakun et al. (2010), does not seem to be required, since the correlation patterns display a large-scale teleconnection structure. Similarly, the AAO index displays suggestive teleconnection patterns with the wind stress field in DJF and JJA (Fig. 11). These patterns indicate that the influence of the AAO is more pronounced in South Benguela and suggests the reason why the correlation between the North Benguela upwelling and the AAO is weaker. In both seasons, the AAO is associated with upwelling-favourable winds in South Benguela, and this link is somewhat stronger in DJF.

Spectral analysis
To identify the more important timescales for upwelling variability, we perform spectral analyses on the upwelling index up_wvelo.
The spectra of indices (Fig. 12) highlight their different behaviour in North and South Benguela. In the north, the spectra are red through the whole frequency range, with some superimposed broad and weak peaks. In the south, the spectra flatten for frequencies lower than 2×10 −1 (periods longer than 5 years). The uncertainty range in the estimation of the spectral density (Jenkins and Watts, 1968) is also included in the figure panels. In both regions the spectra show some enhanced variability that is not strictly statistically significant, but which are nevertheless documented here as it supports some of the results of the correlation analysis and as it can be useful for future analysis with longer simulations in which these signals may be better revealed with a larger sample size.
The up_wvelo index of North Benguela in STORM (Fig. 12) displays a very broad spectral peak at periods of around 5 years in the season DJF and 4 years in JJA. In the season SON, the variabilities with periods of 3 years are slightly enhanced, whereas in MAM this occurs for slightly shorter periods of 2.5 years. In South Benguela, the spectra hint at an enhanced variability in the period range of 3-5 years, with other maxima at 2.5 years. However, whereas the spectra of the up_wvelo index in both North and South Benguela display enhanced variability at the typical ENSO periodicities, correlations to the ENSO index differ between the two regions ( Table 3). The index in North Benguela shows mainly weak correlations to ENSO. The index in the south displays stronger and significant correlations in DJF and MAM.
In summary, the red spectral background in the north is probably a result of tropical dynamics, whereas the more white spectra in the south may reflect the atmospheric forcing at midlatitudes, which generally has a flatter spectrum itself. Both upwelling spectra display enhanced variability at periods of around 2.5, 3.3 and 5 years. The periods of 5 years in the DJF season suggest a link to ENSO, as it was found in the correlation analysis between the upwelling indices in South Benguela and the ENSO index. The 2.5-year period could indicate an influence of the QBO. However, the QBO and up_wvelo of STORM are insignificantly correlated (Table 3). Possible reasons for these variations have to be further investigated. Bakun (1990)  although the upwelling evolution cannot be well described by a steady linear trend. In the south, the long-term evolution is described by a linear trend better than in the North, but its magnitude lies just below the 95 % significance level. The mechanism proposed by Bakun (1990) involves a stronger increase of near-surface temperatures over land than over the ocean, which would lead to an intensification of the continental thermal low pressure relative to the ocean SLP. Consequently, the SLP gradient between land and ocean would tend to increase, the SLP difference between ocean and land becoming more positive, leading to a strengthening of the upwelling-favourable wind and hence to an increase of upwelling. We have tested Bakun's hypotheses concerning the link between the SLP gradient and upwelling in the framework of the ocean simulation (STORM) in Benguela, and concerning the long-term trend in upwelling in this simulation. We defined two regions, one over land (12-25 • E, 10-20 • S) and one over ocean (20 • W-10 • E, 15-35 • S), and calculated the difference (ocean minus land) of their respective area-averaged SLP. Analysing the trend of this SLP difference in NCEP/NCAR reanalysis and ERA-Interim reanalysis does not support Bakun's hypothesis in this region. Bakun's hypothesis would predict a positive trend in this difference -land SLP decreasing relative to ocean SLP -whereas the NCEP/NCAR and ERA-Interim SLP provide significantly negative trends for the MAM seasons (Fig. 13, Table 4), with the trend in other seasons being not significant.

SLP gradient and long-term trends in upwelling
The trends of the SLP over land and over ocean are positive for both data sets and all seasons, except for the SLP in MAM of NCEP (Table 4). This could indicate a missing imprint of increasing greenhouse gas concentrations on the SLP, especially over land, which is may be due to the relatively short temporal period covered by the reanalyses data or the reanalyses data itself.
The correlation at interannual timescales yields a clear relationship between the upwelling intensity and the oceanland SLP gradient. Bakun's hypothesis, applied at interannual timescales, predicts a positive correlation between upwelling and this gradient. The correlations with the up_wvelo are positive and significant in North Benguela in all four seasons (with NCEP/NCAR and ERA-Interim). South Benguela upwelling is significantly correlated with the SLP difference in MAM, JJA and DJF calculated from the NCEP/NCAR SLP gradient, but not in all seasons when correlating with the ERA-Interim SLP gradient. This is probably due to the relatively short period covered by the ERA-Interim. The ERA-Interim covers only 32 years, whereas NCEP spans twice that period. Although the strength of the correlation is similar in both data sets, the smaller sample size in ERA-Interim renders the correlation not statistically significant.
Therefore, the SLP difference between ocean and land derived from these two meteorological reanalyses do not support Bakun's hypothesis of a long-term intensification due to increasing greenhouse gas concentrations, although there is a relationship at interannual timescales as Bakun (1990) and Bakun et al. (2010) envisaged. The lack of a significant trend does not necessarily indicate that Bakun's hypothesis is not correct. The reason could lie in insufficient quality  of the reanalysis SLP in this region to identify long-term trends or that the long-term trend in the SLP gradient over the last decades is still overwhelmed by other factors and therefore has not emerged from the background noise yet. Bakun et al. (2010) indicated that this possible factor which blurs the long-term trend in upwelling in BUS is its relationship to ENSO, to which we briefly direct our attention now.
According to the results of the correlation between upwelling and climate modes, and the correlation patterns between the climate modes ENSO and AAO with the wind stress, it appears that the influence of both climate modes is more pronounced in South Benguela than in North Benguela.
We attempt to estimate the contribution of ENSO and the AAO to the long-term trend in upwelling in South Benguela in DJF, the season in which the interannual correlations are stronger and the wind stress patterns also display a clearer signal. For this, we set up a linear regression model between upwelling as predictand and ENSO and the AAO as predictors: where up_wvelo represents the upwelling index in South Benguela in the DJF, t represents the year, α enso and α aao are the regression coefficients and represents the unresolved variance in upwelling. Using the data from the period covered by all three data sets (ENSO, AAO andSTORM, 1957-2010) we can estimate the values of α enso and α aao . The contribution of ENSO and AAO to the long-term upwelling trend is then estimated as α enso trend enso and (2) α aao trend aao , respectively, where trend enso and trend aao are the long-term trends in the ENSO and AAO indices, respectively. It turns out that the ENSO trend in DJF is by far not statistically significant, whereas the AAO trend in DJF is positive and clearly statistically significant over the p = 0.05 level. The South Benguela upwelling in the STORM simulation increases linearly over 1957-2010 by 39 % of its long-term mean value. The estimated contribution of ENSO to the longterm trend in upwelling in South Benguela in DJF is −3 % of the mean upwelling and the estimated contribution of the long-term trend in the AAO is 41 % of the mean upwelling. This means that, under the assumptions of this simple statistical analysis, the trend in upwelling in South Benguela in DJF simulated in the STORM simulation could be almost entirely explained by the long-term trend in the AAO, with the ENSO contribution remaining negligible.

Discussion
The large-scale atmospheric drivers of Benguela upwelling in this ocean-only simulation were identified using an upwelling index based on the vertical velocity. An important result is the disagreement between our analysis of the SLP gradient and the hypothesis of Bakun (1990). There are several possibilities for this disagreement. The results of the correlation analysis between the upwelling in-dices and the SLP gradient do show the expected relation of wind forcing driving the upwelling on interannual scales, but neither the upwelling nor its drivers, the upwellingfavourable winds and the SLP gradient show significant trends (except the SLP gradient in MAM). A comparison to measurements of the upwelling itself is not possible due to the lack of such observation on multidecadal timescale. Therefore, we compare our results to other studies using reanalysis data, sediment cores or SST observations as upwelling indices. Our results do agree with the findings of Narayan et al. (2010), who did not find a significant trend in upwelling-favourable wind of ERA40. Our results, however, disagree with the positive trends in COADS (Comprehensive Ocean-Atmosphere Data Set) and NCEP/NCAR reanalysis. Nevertheless, the trends in the SST derived from HadISST mentioned by Narayan et al. (2010) strongly depend on the analysed period. The trends of upwelling shown by Narayan et al. (2010) based on sediment records are dominated by multidecadal variability rather by a long-term trend, although analysing only the last 30 years would probably lead to a positive trend. The long-term trend in the SLP may be more strongly burdened by uncertainties in the SLP data from NCEP/NCAR reanalysis in the Southern Hemisphere. For instance, Marshall (2003) found considerable deviations in the magnitude of SLP trends between those derived from the NCEP/NCAR reanalysis and those from station data in the Southern Hemisphere, although the sign of the trend seems to agree in both data sets. Since STORM was driven by the NCEP/NCAR reanalysis, this uncertainty is a strong caveat when estimating long-term upwelling trends. Unfortunately, it is difficult to ascertain which data set (reanalysis, station wind records) is closer to reality, as wind station records may also be strongly influenced by relocation of the stations and instrumental inhomogeneities. Nevertheless, the coarse resolution of NCEP/NCAR reanalysis has to be kept in mind when interpreting our results, not only when looking on the atmospheric variables (SLP and wind stress) but also for the output of the ocean model.
Another oceanic factor that may influence the upwelling is the meridional gradient of the sea surface height (SSH) (Colas et al., 2008). We calculated this gradient in the STORM simulation, defined as SSH in the box 10 • N-10 • S, 50 • W-30 • E minus SSH in the box 30 • S-40 • S, 50 • W-30 • E. This gradient is negatively correlated with the Benguela upwelling. The correlation with North Benguela upwelling is significant with a coefficient of −0.42. The correlation with South Benguela upwelling is not significant with a coefficient of −0.17. This gradient has a positive but not significant long-term trend in the STORM simulation in both seasons JJA and DJF. Therefore, although the sign of the trend of the SSH gradient opposes a positive trend in upwelling by an onshore transport in North Benguela and thus could theoretically contribute to mask the long-term trend in upwelling caused by anthropogenic forcing in this region, its magnitude seems to be too small. In addition, it does not appear related to upwelling in South Benguela.
The SLP mechanisms later augmented by Bakun et al. (2010) involve small-scale dynamics of humidity transport, which may not be properly represented in the coarse resolution atmospheric model used to generate the NCEP/NCAR reanalysis. Unfortunately, no long instrumental SLP records spanning this region are available, so this must remain an open question for the moment. Anther possibility is that, although the effect of the external forcing on the SLP field as envisaged by Bakun (1990) may be correct, the amplitude of internally generated wind variability at decadal and multidecadal timescales blurs the long-term signal. Finally on interannual timescales, as also proposed by Bakun et al. (2010), the influence of other large-scale climate modes, such as ENSO, on BUS upwelling may mask the forcing by anthropogenic greenhouse gases. Our correlation analysis sheds some light on this last question, as explained below.
The statistically significant negative correlation between ENSO and upwelling in South Benguela during the DJF has been reported by previous studies. Rouault et al. (2010) stated that El Niño leads to reduced upwelling whereas La Niña enhances upwelling. The same study argues that the upwelling-favourable winds are more pronounced during La Niña and weaker during El Niño, which explains the SST anomaly in the BUS during ENSO events. The connection between extreme SSTs (in the False Bay, South Africa) and ENSO is also supported by the study by Deutsch et al. (2011) and Dufois and Rouault (2012). The former authors also found significant correlations between El Niño and the first principal component of austral summer, autumn and partly spring SST (strongest with a 4-month lag, which is contained in the seasonal resolution of our analysis).
The spectral analysis emphasises the connection between the upwelling and ENSO in DJF. The source of the other timescales in the spectrum (2.5 and 3.3 years) could be due to ENSO too.
The link between upwelling in South Benguela and ENSO has been invoked to explain the parent lack of a long-term trend in BUS upwelling. Earlier studies, for instance Vecchi et al. (2006), indicated a weakening of the trade winds across the tropical Pacific (i.e. a tendency towards El Niño state) as a response to anthropogenic greenhouse gas forcing. A trend towards more intense or more frequent El Niño phases would thus tend to weaken Benguelan upwelling, counteracting the direct effect on upwelling of the anthropogenic forcing. However, more recent studies of the response of ENSO to anthropogenic greenhouse gas forcing recognise much larger uncertainties in the predictions of the future evolution of ENSO (Vecchi and Wittenberg, 2010), although one recent multimodel study has found that the frequency of extreme El Niño events may increase as a result of anthropogenic greenhouse gas forcing , at least in the first 40 years of the 21st century (Kim et al., 2014). In the observations analysed here, the trend in the ENSO index since 1950 is weak and not significant, so that it is unclear how the long-term trend in ENSO could contribute to robust weakening of upwelling counteracting the greenhouse gas forcing in the past decades. Our analysis rather shows that the climate modes ENSO and the AAO could at most have influenced upwelling in South Benguela and to a lesser extend North Benguela. The observed long-term trend in the indices of these two climate modes could explain almost 100 % of the upwelling trend in DJF in South Benguela in the STORM simulation, but this contribution stems almost entirely from the AAO. In North Benguela, the influence of ENSO and the AAO is very weak and the evolution of the simulated upwelling there cannot be described well by a linear steady trend (Fig. 8) associated with the observed evolution of greenhouse gas forcing.

Conclusions
The large-scale atmospheric drivers of the Benguelan upwelling system, as simulated in a global high-resolution ocean simulation (STORM), driven by a prescribed atmosphere in the last 6 decades, have been described in this paper. The major results are summarised as follows.
-The BUS is better described by two subsystems, North Benguela and South Benguela, as their mean seasonality, their time evolution, and correlations with atmospheric drivers and large-scale climate modes differ.
-As general characteristic for upwelling-favourable atmospheric conditions in both subsystems are an intensified South Atlantic High, strong and southerly winds/wind stress along the coast, and a SLP contrast between land and ocean.
-There is some evidence of the influence of ENSO and the AAO on upwelling in South Benguela. The El Niño phase in the tropical Pacific tends to weakly hinder upwelling, whereas a stronger AAO tends to reinforce upwelling in this region. In North Benguela, there is no clear influence of these climate modes.
-The long-term trends of the simulating upwelling over approximately the last 50 years do not clearly support the hypothesis put forward by Bakun (1990) that anthropogenic greenhouse gas forcing should lead to more vigorous upwelling. The analysis of the trend of the SLP gradient between land and ocean as suggested by Bakun (1990) does not support an influence of the sea-land SLP contrast on the long-term behaviour of upwelling. The estimated influence of the trends of the large-scale climate modes on Benguela upwelling do not indicate that these climate modes may have disturbed the hypothesised connection between anthropogenic greenhouse forcing and upwelling in this region in the last decades.

500
N. Tim et al.: Variability of the Benguela upwelling system -The atmospheric forcing and thus the model results are dependent on the uncertainties inherent in the NCEP/NCAR reanalysis data set. The influence of the atmospheric resolution on the long-term trends in the atmospheric fluxes needs to be reassessed when higherresolution data covering the past decades become available.