Low-frequency variability in North Sea and Baltic Sea identified through simulations with the 3-D coupled physical – biogeochemical model ECOSMO

Here we present results from a long-term model simulation of the 3-D coupled ecosystem model ECOSMO II for a North Sea and Baltic Sea set-up. The model allows both multi-decadal hindcast simulation of the marine system and specific process studies under controlled environmental conditions. Model results have been analysed with respect to long-term multi-decadal variability in both physical and biological parameters with the help of empirical orthogonal function (EOF) analysis. The analysis of a 61-year (1948–2008) hindcast reveals a quasi-decadal variation in salinity, temperature and current fields in the North Sea in addition to singular events of major changes during restricted time frames. These changes in hydrodynamic variables were found to be associated with changes in ecosystem productivity that are temporally aligned with the timing of reported regime shifts in the areas. Our results clearly indicate that for analysing ecosystem productivity, spatially explicit methods are indispensable. Especially in the North Sea, a correlation analysis between atmospheric forcing and primary production (PP) reveals significant correlations between PP and the North Atlantic Oscillation (NAO) and wind forcing for the central part of the region, while the Atlantic Multi-decadal Oscillation (AMO) and air temperature are correlated to long-term changes in PP in the southern North Sea frontal areas. Since correlations cannot serve to identify causal relationship, we performed scenario model runs perturbing the temporal variability in forcing condition to emphasize specifically the role of solar radiation, wind and eutrophication. The results revealed that, although all parameters are relevant for the magnitude of PP in the North Sea and Baltic Sea, the dominant impact on long-term variability and major shifts in ecosystem productivity was introduced by modulations of the wind fields.


Introduction
Long-term variations and major changes in ecosystem dynamics occur throughout all trophic levels and have earlier been reported on in a number of studies for the North Sea and Baltic Sea system (Beare et al., 2004;Beaugrand and Ibañez, 2000;Clark and Frid, 2001;Lynam et al., 2017;Möllmann et al., 2000;Schlüter et al., 2008;Selim et al., 2016;Thurow, 1997;Weijerman et al., 2005;Wiltshire and Manly, 2004).A majority of those studies have focused on potential "regime shifts" (RSs, "changes in marine system function that are relatively abrupt, persistent, occurring at a large spatial scale, observed at different trophic levels and related to climate forcing"; deYoung et al., 2004).Such major changes throughout all trophic levels were reported for the North Sea and Baltic Sea system at the end of the 1980s (Alheit et al., 2005;Österblom et al., 2007;Weijerman et al., 2005).Beaugrand (2004) reviewed studies addressing RSs in the North Sea.He reported on studies considering temporal changes in single-species abundance and vital rates throughout all trophic levels, system productivity and species composition within trophic levels or feeding guilds.By com-Published by Copernicus Publications on behalf of the European Geosciences Union.
U. Daewel and C. Schrum: Low-frequency variability in the North Sea and Baltic Sea bining these studies with time series information on hydrometeorological conditions for the same time periods, Beaugrand (2004) hypothesized three different drivers for persistent changes in the North Sea ecosystem: (i) a change in the local hydro-meteorological forcing, (ii) a displacement of oceanic biogeographical boundaries, and (iii) an increase in oceanic inflow into the North Sea.Dippner et al. (2012) compared potential regime shifts in the North Sea and Baltic Sea and could associate the inter-annual variability and RSs in the Baltic Sea to changes in the atmospheric forcing only, while for the North Sea he found combined influences from the atmospheric and the Atlantic forcing to be most likely responsible for inter-annual variations in ecosystem dynamics.In fact, many studies could relate variations in the ecosystem to variations in atmospheric variables and indices, such as NAO, sea surface temperature (SST) and wind (Alheit et al., 2005;Beaugrand and Kirby, 2010;Edwards et al., 2010) but also to modification in the anthropogenic forcing such as fisheries or river nutrient loads (Österblom et al., 2007).Nonetheless, the identification of causal relationships and underlying processes is difficult based on in situ observations only, due to the complexity in identifying the relative relevance of single factors (Clark and Frid, 2001) and also the inhomogeneous characteristics of the datasets, which are often relatively short and lack the spatial diversity in regional ecosystem components.
However, understanding the relevance of environmental factors for ecosystem dynamics pioneers the identification of environmental indicators for long-term variations and RSs."Indicators are proxies for complex phenomena and can be used to reflect the provision of a service and how it is changing over time" (Hattam et al., 2015).Hence, the identification of potential indicators is of major relevance for both marine ecosystem understanding and management.Since bottom-up processes play a major role for long-term variations in functioning of many regional marine ecosystem, and the North Sea and Baltic Sea system in particular (Daewel et al., 2014;Frank et al., 2007), understanding processes impacting net primary productivity form the basis for indicator definition.To overcome the limitation of observationally based analysis, coupled physical-biological ecosystem models are valuable tools that provide spatially explicit long-term datasets of lower-trophic-level production (Daewel and Schrum, 2013).Additionally, these kinds of models allow a clear analysis of environmental factors and underlying mechanisms since the former are explicitly prescribed in the model formulation and set-up.Additionally, specific scenarios can be applied by artificially modulating the forcing parameters to test hypotheses and indicators.
Here we further analysed the 61-year simulation , which was earlier presented by Daewel and Schrum (2013).The length of the simulation period allows the identification of long-term changes in the environment and in primary production in the North Sea and Baltic Sea.Here, we aim at exploring key long-term variation in rele-vant environmental variables and the potential of different methods to derive environmental indicators describing the hydrodynamic and biogeochemical environment.We evaluate the potential of aggregated hydrodynamic, atmospheric and large-scale climatic indicators to explain modelled primary production variability.Finally, we utilize the model to simulate specific scenarios to understand the causal relationship between indicators and the low-frequency variability in simulated primary production.

ECOSMO II model description
ECOSMO II (ECOSystem MOdel; Daewel and Schrum, 2013;Schrum et al., 2006a) is a 3-D fully coupled physicalbiogeochemical model.The long-term simulation of lowertrophic-level ecosystem dynamics with ECOSMO II was presented and validated in Daewel and Schrum (2013).The hydrodynamic core of the coupled model system is a mature and validated (e.g.Janssen et al., 2001;Schrum, 2001) 3-D baroclinic coupled sea ice model based on the version of the HAMSOM (HAMburg Shelf Ocean Model) presented first by Schrum (1997) and Schrum and Backhaus (1999).The model is a free-surface model and allows for variable bottom layer thickness; hence, it resolves a realistic bathymetry.The model uses semi-implicit methods (Backhaus and Hainbucher, 1987), which allows for a relative large model time step of 20 min.In contrast to the earlier model version described by Schrum and Backhaus (1999), here we use a second-order total variation diminishing (TVD) scheme, namely the second-order Lax-Wendroff scheme, which was made TVD by a superbee limiter (e.g.Harten, 1997) for the advection of all scalar properties.Its implementation and the consequences for ecosystem dynamics are described in more detail by Barthel et al. (2012).The model equations are solved on a staggered Arakawa C grid for the North Sea and Baltic Sea, with a horizontal resolution of 6 nm (1 nm = 1852 m) and 20 vertical levels, whereof the upper 40 m have a 5 m resolution to resolve stratification.The model was used earlier to investigate seasonal and interannual to decadal variations in stratification and has been found to successfully reproduce the latter in the North Sea (Janssen et al., 2001;Schrum et al., 2000).
The biogeochemical processes in ECOSMO II were simulated using 16 state variables to resolve ecosystem dynamics by a functional group approach (Fig. 2).The model estimates two zooplankton functional groups; three phytoplankton groups; the nitrogen, phosphorus and silicon cycle; oxygen; detritus; biogenic opal; dissolved organic matter; and three sediment groups.The model equations, set-up and a model validation for a 61-year model hindcast integration were presented in detail by Daewel and Schrum (2013), who found the model able to reproduce temporal and spatial variability in primary and secondary production of the North Sea   (Daewel and Schrum, 2013).
and Baltic Sea on intra-and inter-annual timescales up to decadal timescales.The model was validated using nutrient observations only because of its good availability, reliability and comparability to model-estimated nutrients.
Atmospheric boundary conditions are required at the airwater interface and were taken from the NCEP/NCAR reanalysis data (Kalnay et al., 1996).Sea surface elevation, including the major tidal constituents as well as salinity and nutrients, was prescribed at the open boundaries to the North Atlantic (see Fig. 1).For the remaining ecosystem variables and temperature, a Sommerfeld radiation condition was applied at the open-ocean boundaries (Orlanski, 1976).Additionally, river runoffs and nutrient loads are given at the landocean boundary from a collection of different data sources.For more details on data sources and handling and a complete description of the simulation set-up, please consult Daewel and Schrum (2013).

Statistical methods
The advantage of model-derived data is their spatially and temporally explicit characteristics, which allows us to resolve the variability on various time and spatial scales.To identify major modes of variability, we apply a widely used method in climate and ocean science, the empirical orthogonal function (EOF) analysis, a statistical method to identify dominant modes of variability in multidimensional data fields (e.g.Storch and Zwiers, 1999;Venegas, 2001).Here the method is used to understand and compare major modes in the hydrographical and ecosystem components of the cou-U.Daewel and C. Schrum: Low-frequency variability in the North Sea and Baltic Sea pled marine system, namely for the mean winter (January-March) current field and net annual primary production, and to statistically compare these modes to potential driving environmental variables.
The method is comparable to the one used in Daewel et al. (2015), who gave the following brief introduction into the main elements of the analysis to clarify the terms used in the analysis.
The annual values of the spatially explicit variable field form a N × M matrix χ (N: number of years; M: number of wet grid points).The empirical modes are given by the K eigenvectors of the covariance matrix with non-zero eigenvalues.Those modes are temporally constant and have the spatially variable pattern p k (m = 1, . .., M) where k = 1, . .., K.The time evolution A k (t = 1, . .., N ) of each mode can then be obtained by projecting p k (m) onto the original data field χ such that χ (t, m) = K k=1 p k (m) A k (t).In the following we will refer to A k (t) as the principal components (PC) and to p k (m) as empirical orthogonal function (EOF).The percentage of the variance of the field χ explained by mode k is determined by the respective eigenvalues and is referred to as the global explained variance η g (k).Before using the method to analyse the spatiotemporal dynamics of the field, the data were demeaned (to account for the variability only) and normalized (to allow an analysis of the variability independent of its amplitude).The identified modes are not necessarily equally significant in all grid points of the data field.Thus, the local explained variance η local,k (m) could provide additional information about the regional relevance of an EOF mode and the corresponding PC in percent: where Var (X) = N t=1 X − X (t) 2 denotes the variance of the field X(t).
Note that in our study the data were additionally lowpass filtered using a 5-year running mean prior to applying the method.The principal modes of the EOF analysis are purely mathematical and not necessarily related to dynamical processes or physically interpretable.However, the use of a proper regional and temporal window encompassing the potential scales of variability in the targeted parameter improves the potential for several dynamically relevant modes (Schrum et al., 2006b).
Subsequent to the EOF analysis, the major principal components (PCs) were compared through correlation analysis to equally low-pass-filtered time series of environmental variables to identify potential environmental indicators and underlying processes.A Pearson correlation coefficient was estimated and tested against a t distribution to obtain a measure for significance (Storch and Zwiers, 1999).A list of tested environmental variables is given in Table 1.These variables were averaged in time (see Table 1) and space (North Sea and Baltic Sea) prior to analysis.

Scenario simulations -design
Three types of scenarios were designed to target the specific hypothesis deduced from the statistical analysis of model results and previously published hypotheses on processes behind ecosystem changes in the North Sea and Baltic Sea (see introduction).Here, we tested (i) the impact of short-wave radiation (SWR) as a parameter determining the season length and intensity of the annual primary production but which also plays a role in changes in water temperature and mixed-layer depth (MLD); (ii) the impact of the wind forcing, which affects not only the general current field and nutrient supply from the open ocean to the North Sea but also vertical mixing and upwelling and hence mixing of nutrients to the euphotic layer; and (iii) the ecosystem response to changes in the river nutrient loads.
Instead of just increasing or decreasing the magnitude of the forcing parameters by a certain percentage, we aimed at resolving the impacts of the multi-decadal variations in major shifts in the ecosystem dynamics.The first analysis identified that the 61-year simulation period covered two different 30-year periods, for which productivity was significantly different (Daewel and Schrum, 2013).To identify the driving mechanisms for this change, we divided the 61-year simulation period into two climatic sub-periods (TP1: 1948-1976and TP2: 1980-2008).Two climatic forcing variables were tested: SWR (sr) and wind stress (wi).For each of these two, scenario simulations were performed, for which all forcing variables but the target variable were kept unchanged with respect to the reference simulation.For the target variable, the forcing was repeatedly employed for both sub-periods (Fig. 3) such that in simulation 1 (sr1/wi1) the forcing from the TP1 was repeated in TP2, and in simulation 2 (sr2/wi2) the forcing from TP2 was also applied to TP1.
For the third set of scenarios we estimated average seasonal cycles for the river nutrient loads (NO 3 , PO 4 , SiO) in each of the six decades (Fig. 4) and performed a set of six simulations each forced by a different river load climatology.This enables us to explore the relevance of different persistent nutrient load situations and its relevance for abrupt changes in the system.The scenarios chosen include relatively high (80-89), intermediate (90-99) and low (00-08) nutrient loads, but also unusual N / P ratios in the forcing (70-79).
Table 1.Variables used for correlation analysis with principal components of the net primary production EOF analysis (Figs. 7 and 8).Both atmospheric and oceanic variables were averaged over the respective subregion (North Sea or Baltic Sea) for the analysis.

Name
Explanation  The set-up is valid for the short-wave radiation experiments (sr1 and sr2) and for the wind experiments (wi1 and wi2); ref denotes the reference simulation as described in Daewel and Schrum (2013).

Environmental indicators
To identify key long-term variations occurring in the North Sea and Baltic Sea system, we first investigated spatial averages of temperature, salinity and current speed for key regions.We focus here especially on the variations in the North Sea and present analysis in upper and lower water layers for the northern and southern North Sea (Figs. 5 and  6).Our analysis highlights several key characteristics related to long-term variations in hydrodynamics in the North Sea.Specifically, we find the following: an increase in temperature since the beginning of the 1990s was simulated for both the northern and southern North Sea SSTs and bottomwater temperature (Fig. 5).In the southern North Sea, trends in surface and bottom layer are similar.However, this is not the case in the northern North Sea where temperature varies independently for surface and bottom waters.Substantial multi-year variations superimpose the long-term trends in the North Sea temperature and are evident in both the surface and bottom layers.Additionally, surface water temperatures are also characterized by biennial periodicity.However, in the shallow southern North Sea the latter variations are also shown for the bottom layer, indicating a stronger coupling between surface and bottom in that region; the bottom layer of the deeper northern North Sea is largely uncoupled from these variations.Also, salinity patterns are dominated by long-term and decadal oscillations, whereof no long-term trend but rather multi-decadal variation is found in the northern North Sea.The southern North Sea, in contrast, features an increasing trend in surface salinity, accompanied by a slightly weaker increase in bottom-water salinity.Multiyear variations in salinity are comparable to those in temperature, but the strong biennial periodicity in surface temperature is not similarly evident for salinity, for which interannual and decadal to multi-decadal variability dominates.Current speed in the North Sea (Fig. 6) is dominated by a multi-decadal sinusoidal variation with low current speeds in the first 3 decades of the simulation period and higher current speeds in the later 3 decades.A contrasting trend is however found for the northern North Sea bottom layer, showing a period of minimum current speed in the intermediate simulation period .Here, a strong coupling between variability in surface and bottom layer is again identified.
The potential of statistical analysis to provide more detailed information on long-term variations in North Sea and Baltic Sea currents is explored through EOF analysis of current vectors.In Fig. 7 we present the mean (averaged over the 61-year time period) surface current field in the North Sea and Baltic Sea and the dominant mode from an EOF analysis over the anomalies to the mean current vector field for the winter season.The analysis indicates a substantial winter inflow anomaly in the North Sea, with current speeds from northwest to southeast during the last 2 decades.Contemporaneously, the Baltic Sea was characterized by a substantial cross-basin circulation anomaly from the Swedish coast towards the Polish coast that was likely related to a substantial ventilation of the Baltic Sea and nutrient transport from the lower layers to the euphotic zone as a consequence of enhanced coastal upwelling.This nutrient enhancement in the surface would foster the Baltic Sea primary production, a development that was indeed modelled (compare Fig. 10I and explanation below).Additionally, we find substantial decadal variability in the circulation.The first EOF thereby covers a significant part of the overall variability, with more than 60 % explained global variance.An additional EOF analysis performed for the scalar current speed further highlights the fact that this strong increase in strength of the northwest current component is connected to a general increase in current speed (Fig. 8c).The local explained variance in the first EOF mode (Fig. 8b) shows that this dominant mode of variability (Fig. 8a) is highly relevant in the central and northnorthwestern parts of the two main areas in the coupled North Sea and Baltic Sea system.However, it does not explain variability in the southern and eastern coastal regions nor in the Bothnian Bay and Gulf of Finland, indicating that the current speed variability in these areas differs substantially from the dominant pattern.

Ecosystem variability
As highlighted above, changes in environmental variables are hypothesised to play a crucial role in explaining long-term changes in North Sea and Baltic Sea ecosystem dynamics.Here, we aim at identifying hydrodynamic and atmospheric indicators, which could serve as potential predictors for spatially resolved primary production changes.A number of indicators were tested, covering large-scale climate, regional atmospheric and regional hydrodynamic indicators.The predictive potential of these indicators was tested and comparatively assessed through correlations to the major PCs of primary production estimates (Figs. 9 and 10).
In the North Sea (Fig. 9) the first and second EOFs explain the variability in the central North Sea and in the southern frontal areas (Fig. 9I and II), featuring substantially different temporal variability (PC 1 in Fig. 9Ic and PC 2 in Fig. 9IIc).While in the central North Sea a major shift in primary production was simulated around 1980 (PC 1 ), the production in the frontal regions passed through two major changes (around 1970 and around 1990) (PC 2 ).In general the signals (PC 1 and PC 2 ) were overlaid by a quasi-decadal variability, which is comparable but not identical (partly caused by the statistical filtering procedures) to the variability estimated for the wind field.
The correlation analysis (Fig. 9III) reveals that the potential indicators for production are very different for the two patterns (relevant in the different subregions).For the central North Sea, for which variability is mainly described by the first PC (PC 1 , Fig. 9Ic), changes in the North Atlantic Oscillation (NAO), changes in wind speed, specifically the western and southern wind components, and, associated with wind speed, changes in current speed show the highest correlations to the major mode of variability in primary production, although several other variables are also significantly (at the 5 % level) correlated to PC 1 (including SWR, winter vertical velocity, surface salinity, PO 4 and NO 3 ).The production changes in the frontal areas (PC 2 , Fig. 9IIc   the 25 considered environmental variables.The highest correlations in the frontal areas can be found for the Atlantic Multi-decadal Oscillation (AMO), air temperature, and precipitation.For the oceanic variables, the highest correlations can be found for SST and the stratification index early in the season (MLD_May).Despite the difference in regional and temporal variability, for both PCs the most significant indicators are linked to processes driving the surface nutrient concentration, which is meaningful in a system in which upper-layer primary production is limited by nutrient availability.Here, the two identified regions are influenced by different processes.(i) Processes related to EOF 1 /PC 1 : the longterm variability in the seasonally stratified central North Sea is mainly related to wind stress, which determines the nutri-ent inflow from the North Atlantic to the North Sea but also impacts vertical mixing and nutrient supply to the surface layer.(ii) Processes related to EOF 2 /PC 2 : in the frontal areas off the Danish and English coasts and at Dogger Bank, the long-term changes in primary production are negatively correlated to the AMO, air temperature and precipitation.The latter two variables impact the strength and timing of the seasonal stratification.Here the effect is inversely proportional: the warmer the temperatures, the stronger the stratification.
Especially in regions with intermediate depths, a strong stratification and an early onset of the stratification could substantially limit the nutrient supply to the euphotic zone.
In the Baltic Sea (Fig. 10), almost 70 % of the overall simulated variability in primary production is described by the first EOF mode and PC (Fig. 10I).Here, we see a clear increase in primary production for the time period 1950-1987 and an abrupt increase thereafter followed by an ever-soslight decrease in primary production (Fig. 10Ic).The steep increase at the end of the 1980s has been shown to statistically significantly differentiate two different periods (Daewel and Schrum, 2013) and clearly corresponds to the time ear-lier described for a regime shift in the Baltic Sea (Alheit et al., 2005).Daewel and Schrum (2013) showed that significant changes were evident for all three phytoplankton functional types but that changes in cyanobacteria and flagellate production contributed most to the overall change.Hence, it is not surprising that surface PO 4 shows the highest correlation (R = 0.97) to the production change (Fig. 10III), and thus processes impacting the latter must play a significant role in primary production in the Baltic Sea.Nonetheless, in contrast to the North Sea, the correlation analysis for the Baltic Sea PC 1 did not indicate a dominant factor or process that could serve as an environmental indicator for production since most of the considered parameters were found to significantly correlate to the main temporal changes in primary production (Fig. 10III).Additionally to the winter NAO, both wind speed and SWR are highly correlated to the major production pattern (PC 1 ).In contrast, the AMO was one of the few parameters with no significant correlation.The second EOF (Fig. 10II) is less distinct and explains only about 6 % of variability, mostly in some coastal areas and in the Gulf of Bothnia (Fig. 10IIb).For the related PC 2 (Fig. 10IIc) no clear relationships could be identified.

Causal relationships
Since correlation analysis can identify statistical relations but not causality, we compiled subsequent scenario experiments with the model to identify the role of variations in wind speed, SWR and river nutrient loads for production changes in the North Sea and Baltic Sea.Those parameters were chosen due to the high correlation we found between primary production and dynamic variables related to wind field changes (wind speed, wind components, current speed) and SWR.The latter showed particularly high correlation to Baltic Sea production variability.River loads were earlier hypothesized as one of the most relevant factors responsible for Baltic Sea system state changes from the late 1960s onwards (Thurow, 1997) and for production changes in the southern North Sea (Clark and Frid, 2001).To emphasize the changes in variability rather than magnitude, the temporal variability in the single forcing parameters were modified as described in Sect.2.3 (see Figs. 3 and 4).In Fig. 11 average low-pass-filtered time series for net primary production in the North Sea (southern North Sea in Fig. 11a and northern North Sea in Fig. 11b) and Baltic Sea (central Baltic Sea in Fig. 11c and Gulf of Finland and Gulf of Riga in Fig. 11d) are shown for the reference simulation and for the different scenario simulations.What becomes evident from this comparison is that the SWR forcing (sr1/sr2), although highly correlated to the Baltic Sea productivity and in addition to nutrient availability, one of the main limiting factors for primary production, changes surprisingly little of the low-frequency variability in both North Sea and Baltic Sea productivity.Despite some small changes in short-term variability, especially in the southern North Sea, the multi-decadal variability and the major shifts remain unchanged in all sub-areas.The wind forcing (wi1/wi2), on the contrary, can clearly be held responsible for structuring the long-term variation.Most notably, our results indicate that the appearances of major shifts in the system (around 1980 in the North Sea and at the end of the 1980s in the Baltic Sea) are mainly caused by changes in the wind field, while the quasi-decadal variations in the signal seem to remain largely unchanged.Note that we cannot exclude that the quasi-decadal variations in the newly compiled wind scenarios are coincidentally in phase with the variations in the reference forcing.Hence, this finding is no indication that the quasi-decadal variability is not attributed to wind field variations.However, in all four sub-areas the regime shifts in productivity are eroded or shifted in time when an alternative wind forcing is applied.This becomes most evident in the northern North Sea (Fig. 11b) and in the central Baltic Sea (Fig. 11c), where the long-term production variability quite closely follows the variability in the wind field and sea surface current speed (compare also Fig. 5 and correlation analysis in Figs. 8 and 9), and the major shift in experiment wi1 is displaced to the end of the 1990s following the wind forcing dynamics from the TP1.Similar to the SWR experiments, variation in the river nutrient loads does not substantially change the long-term variability in ecosystem productivity in either the North Sea or the Baltic Sea.However, it is shown that river loads clearly have an impact on the magnitude of the production in all areas, but especially in the Gulf of Finland and Gulf of Riga (Fig. 11d) region that features major river inflows.Clearly, nutrient loads from the 1980s are highest, resulting in higher system productivity.The comparison to the reference run shows that the river nutrient forcing does not cause major shifts in ecosystem productivity but can clearly amplify changes in the sys-  1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995  tem as seen in the two North Sea regions, where the production increase in the beginning of the 1980s is substantially enhanced by the high river nutrient loads in that decade.Interestingly, in the central Baltic Sea this effect is not similarly apparent.Here changes in nutrient loads aggregate and result in lower or higher production, with the changes increasing slowly over time.

Discussion and conclusion
We identified long-term multi-decadal variations in temperature, salinity, currents and primary production in the North Sea and Baltic Sea from a coupled biological-physical model simulation (Daewel and Schrum, 2013).While Daewel and Schrum (2013) already identified multi-decadal changes in simulated long-term dynamics of ecosystem productivity in the North Sea and Baltic Sea, the causes and underlying processes were only speculated on in their paper.One of the major advantages of coupled ecosystem models is the availability of all information relevant for the system dynamics, including physics and forcing variables.Thus, underlying process interactions can be obtained via statistical analysis and scenario simulations.
As already shown by Janssen et al. (2001), the model is able to simulate long-term dynamics in physical parameters.In this study we investigated average long-term changes in temperature, salinity and current speeds for the North Sea system.We also find here the long-term dynamics in temperature and salinity to cover average variability in observed temperature (Edwards et al., 2010) and salinity, by representing, for example, the "great salinity anomaly" as observed between 1977 and 1981 in the North Sea (Danielssen et al., 1996).In addition to temperature and salinity, current fields have been hypothesised to play a dominant role in ecosystem functioning.Here, average surface current fields for the northern and southern North Sea were identified to follow similar long-term dynamics with a clear increase in current speed starting already in the beginning of the 1970s.This pattern is a result of the changing wind forcing above the North Sea as shown by Siegismund and Schrum (2001), who reported an intensification of west-southwesterly wind directions, an almost linear increase in wind speed and a more frequent appearance of "strong wind" events since the early 1970s.The same authors reported "an extension of winterly wind climate towards February and March during the last (analysed) decade (1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997) 2015) also found the changes in the circulation to be highly correlated to changes in the NAO.Their analysis showed that under stronger and more frequent westerly wind conditions the North Sea inflow through the Fair Isle Passage was particularly enhanced, fostering a stronger southward flow of Atlantic water masses along the British east coast.Under opposing weather conditions, the circulation in the central and southern North Sea weakens and the inflow through the Fair Isle Passage follows the Dooley current.Thus, in that way, "effectively decoupling the water masses of the central and southern North Sea from the northern inflow" (Mathis et al., 2015).This process proves especially relevant for the central North Sea, which is, in contrast to well-mixed areas of the southern North Sea, neither strongly exposed to water inflowing from the English Channel nor to river run-offs and can hence serve as an explanation for the provided correlation between the first mode in North Sea primary production variability and the NAO and wind field.
Applying EOF analysis to primary production allows the identification of major modes of variability and their pattern together with a local indicator of explained variance.Here, the North Sea and Baltic Sea analyses led to very different results.While in the Baltic Sea we found one dominant mode that explains 67 % of the overall variability in primary production, the North Sea variability is spatially more diverse and we could identify at least two dominant modes of variability linked to specific spatial hydrodynamic features of the North Sea as described in Otto et al. (1990).Although commenting on the occurrence and relevance of actual regime shifts in the North Sea and Baltic Sea is beyond the scope of our model, the estimated primary production analysis indeed indicated major shifts for the times when regime shifts have been identified in the literature (e.g.Dippner et al., 2012;Weijerman et al., 2005).Hence, our findings can be considered relevant for explaining major indicators for RSs in the area.Clearly the results from our study indicate that analysing long-term variability in ecosystem dynamics for an average North Sea system is not sufficient.From the regime shifts detected in the North Sea, the change in 1978-1979 appears dominantly in the central North Sea (as indicated by the dominant mode of variability), while the second mode, rele-vant in the southern North Sea frontal areas, would at least show a stronger decrease in primary production around 1990 when the second regime shift is presumed.While the second mode was correlated to air temperature and precipitation, environmental variables that affect the oceanic mixedlayer depths, the first mode is clearly correlated to changes in the wind and current field and resembles the variability in average sea surface currents (compare Fig. 6 and explanation of North Sea circulation).As already described above, the main processes relevant for low-frequency variations in primary production of the North Sea and Baltic Sea are specifically those impacting nutrient supply in the euphotic zone.Although this is in line with what has been reported for the dynamics of the 1978-1979regime shift in Dippner et al. (2012)), the variability for the central North Sea was, in contrast to their explanations, not correlated to the AMO nor to changes in the air temperature.Our results also did not support the hypothesis that changes in salinity (Lindeboom et al., 1995) or changes in sunspot activity (results not shown) (Weijerman et al., 2005) caused changes in ecosystem dynamics.However, the identification of indicators for longterm variation assumes a priori that the indicator remains relevant for the entire time period, while studies tailored to regime shifts usually do not consider the impact on the longterm dynamics and hence might come to different results.
The Baltic Sea primary production dynamics is linked in almost the entire basin to changes in the wind field.This was particularly evident from the performed scenario runs showing that, although nutrient loads alter the magnitude of the primary production, the wind fields determine the timing and magnitude of long-term variations.In Daewel and Schrum (2013), it was already pointed out that the production variability is mainly seen in the flagellates and cyanobacteria bloom, while the analysis presented here indicates linkage to the winter current field (compare Figs. 7 and 8).In principle the underlying process can be explained by the cause-and-effect chain proposed by Janssen et al. (2004) and the preconditioning of the deeper water column phosphate concentrations through eutrophication and anoxic conditions (Rodhe et al., 2006), which is additionally mediated by atmospheric conditions (Schinke and Matthäus, 1998).Thus, our results would support the hypothesis that long-term changes in primary production of the Baltic Sea are a consequence of eutrophication, even though the latter does not serve as a respective indicator for abrupt regime shifts.A similar argument has been formulated in the regime shift analysis by Österblom et al. (2007).
Here, we can conclude that changes in the wind speed and/or changes in the east-west component of the wind field can serve as an indicator or maybe even as a predictor for changes in primary production in both targeted areas.Even in the southern North Sea the changes in wind fields explain more of the long-term production changes than variations in the nutrient forcing, which would, at least partly, contra-dict conclusions from Clark and Frid (2001) on the southern North Sea phytoplankton dynamics.
However, it should be pointed out that this analysis is performed to identify indicators for low-frequency variability; correlations are substantially weaker in unfiltered time series.Moreover, climatic conditions might change and the relevance of specific processes for inter-annual changes in production can alter due to changes in environmental and climate conditions.An example from our model study are variations in North Sea nutrient loads, which caused an amplification of the wind-induced variations in the 1980s in the northern North Sea as well as alterations of the primary production variability in the southern North Sea after 1990 when nutrient loads were substantially reduced.Other possible examples are changes in stratification and, at least in the Baltic Sea, sea ice retreat that could cause variations in primary production and become more relevant under future climate, in which case air temperature or short-wave radiation could become a more significant indicator than wind speed.

Figure 1 .
Figure 1.Model area and bathymetry.Black lines indicate the 30 and 60 m depth lines (the 60 m depth line separates the northern and southern North Sea; central BS includes all areas east of 14 • E excluding the gulf regions).

Figure 3 .
Figure 3. Schematic diagram for the scenario simulation set-up.The set-up is valid for the short-wave radiation experiments (sr1 and sr2) and for the wind experiments (wi1 and wi2); ref denotes the reference simulation as described inDaewel and Schrum (2013).

Figure 4 .
Figure 4. Decadal mean annual nutrient loads N tot (NO 3 +NH 4 )and PO 4 averaged for each of the six simulation decades for use in the scenario simulations.Note: SiO has also been modified but is not shown here.

Figure 5 .
Figure 5. Northern North Sea (left two columns) and southern North Sea (right two columns) temperature (a, b, c, d) and salinity (e, f, g, h) in the surface (left) and bottom layers (right).Displayed are monthly data as a 13-point moving average (black) and 61-point moving average (red).

Figure 6 .
Figure 6.Northern (a, c) and southern (b, d) North Sea surface (a, b) and bottom current speed (c, d).Displayed are monthly data as a 25-point moving average (black) and a 61-point moving average (red).

Figure 7 .
Figure 7. Mean surface current vectors in the North Sea (a) and Baltic Sea (b); EOF analysis of the anomalies in current vectors for the winter period January-March: current pattern for the first EOF (c, d) and first principal component (e).

Figure 8 .
Figure 8. EOF analysis of the anomalies in current speed for the winter period January-March: (a) current speed pattern for the first EOF; (b) local explained variance; (c) first principal component.

Figure 9 .
Figure 9. (I and II) (a) First and second empirical orthogonal functions for annual mean primary production in the North Sea (1948-2008); (b) local explained variance for the pattern for the corresponding EOF; (c) principal component (time variation) of the corresponding EOF.(III) Absolute values of the correlation coefficient between the principal components (PC1 & PC2) and an environmental variable stated on the x axis.

Figure 10 .
Figure 10.(I and II) (a) First and second empirical orthogonal functions for annual mean primary production in the Baltic Sea (1948-2008); (b) local explained variance for the pattern for the corresponding EOF; (c) principal component (time variation) of the corresponding EOF.(III) Absolute values of the correlation coefficient between the principal components (PC1 & PC2) and an environmental variable stated on the x axis.

Figure 11 .
Figure 11.Estimated net primary production for the reference run (ref) and the scenario simulations concerning short-wave radiation (sr1 and sr2) and wind (wi1 and wi2) (a, b) and river nutrient load (Cl) (c, d) for two subregions in the North Sea (southern and northern North Sea) and two subregions in the Baltic Sea (central Baltic Sea and Gulf of Finland and Gulf of Riga).

Figure 2. Schematic diagram of biological-chemical interactions in ECOSMO II
• E excluding the gulf regions).