Towards an integrated forecasting system for pelagic fisheries

First results of a coupled modeling and forecasting system for the pelagic ﬁsheries are being presented. The system consists currently of three mathematically fundamentally di ﬀ erent model subsystems: POLCOMS-ERSEM providing the physical-biogeochemical environment implemented in the domain of the North-West European 5 shelf and the SPAM model which describes sandeel stocks in the North Sea. The third component, the SLAM model, connects POLCOMS-ERSEM and SPAM by computing the physical-biological interaction. Our major experience by the coupling model subsystems is that well-deﬁned and generic model interfaces are very important for a successful and extendable coupled model framework. The integrated approach, sim-10 ulating ecosystem dynamics from physics to ﬁsh, allows for analysis of the pathways in the ecosystem to investigate the propagation of changes in the ocean climate and lower trophic levels to quantify the impacts on the higher trophic level, in this case the sandeel population, demonstrated here on the base of hindcast data. The coupled forecasting system is tested for some typical scientiﬁc questions appearing in spatial 15 ﬁsh stock management and marine spatial planning, including determination of local and basin scale maximum sustainable yield, stock connectivity and source/sink structure. Our presented simulations indicate that sandeels stocks are currently exploited close to the maximum sustainable yield, but large uncertainty is associated with determining stock maximum sustainable yield due to stock eigen dynamics and climatic 20 variability. Our statistical ensemble simulations indicates that the predictive horizon set by climate interannual variability is 2–6 yr


Introduction
In recent years evidence of degradation of marine ecosystems by overfishing, bycatch, climate change, eutrophication and chemical pollution from land runoff, coastal development, habitat destruction and other human activities has become undisputable (Levin and Lubchenco, 2008).Why are we still uncertain about how best to manage the crisis of ecosystems, especially fish communities, triggered directly or indirectly by human impacts?Forecasting fish stocks akin to the ubiquitous weather forecasts has proved to be an elusive target for fishery science (Huse and Ottersen, 2003;Megrey et al., 2005).This is fundamentally related to the character of the ecosystem dynamics where the processes are poorly constrained, quantified and partially still poorly understood and involve ranges of temporal/spatial scales that span several orders of magnitude, generally non-linearly linked.A key problem here remains the insufficient field sampling of ecosystem state variables to parameterize processes, define its initial state and boundary conditions as well as insufficient spatial resolution of ecosystem models (Fulton et al., 2003).Consequently fundamental difficulties in analyzing ecosystem behavior arise and hypotheses are numerous and difficult to discriminate (Carpenter, 2002), further complicated by mathematical issues typical of non-linear systems like instabilities, chaos and regime shifts.These pessimistic prospects have spurred some interests into alternative routes to ecosystem forecasting, like artificial neural networks (Chen and Ware, 1999;Huse and Ottersen, 2003;Suryanarayana et al., 2008).However.ecosystems are facing regime shifts and major trophic reorganization these years (Beaugrand et al., 2002;Alheit, 2009;Kenny et al., 2009;Moellmann et al., 2009), which renders unguided machine learning from recent history dubious; there is no safe way around improving a process-oriented representation of the marine ecosystems; additionally in this paper we advocate data assimilation as a path for reducing model error.
Despite not being on the mark yet, many advances in process-oriented ecosystem modelling have been presented recently (Fulton et al., 2011;Hinrichsen et al., 2011).Introduction

Conclusions References
Tables Figures

Back Close
Full Many of these advances are supporting operational systems of physical regional circulation models and associated biogechemical models describing the lower trophic ecosystem levels (Gallego et al., 2007;Brasseur et al., 2009) which provide access to high quality input data for fish stock models.The process of integrating operational oceanography products with fish stock model has been accelerated by projects like MyOcean (MyOcean, 2009(MyOcean, -2012)).
To achieve a better understanding of the full system including the higher trophic levels it is attractive to focus on fish species with comparatively simple life cycles and then, once under control, move to fish species displaying more complex traits as well as interacting fish stocks.Here sandeel and other sedentary species with well-characterized habitat requirements appear ideal, because stock migration is a very challenging issue to model (Kishi et al., 2011).Ecologically, sandeel is an important mid trophic wasp-waist species, capable of exerting both bottom-up and top-down control, in that it constitutes a significant part of the fish biomass (∼25 %) in the North Sea ecosystem.Therefore predicting and understanding the dynamics of this stock has attracted some interest (van Deurs et al., 2009;Arnott and Ruxton, 2002;Lewy et al., 2004), however most of these works are statistical approaches.Lewy et al. (2004) found that the predictability of the all North Sea sandeel stocks together was better than that of individual sub stocks; our work will challenge this finding.The short life cycle of the species makes it especially susceptible to spatial and temporal variability in ambient physical conditions, because the biomass and recruitment are controlled by one or two offspring cohorts (Arnott and Ruxton, 2002).For instance, in 2010 the proportion of fish just one year old in the catch was more than 90 % and such high proportion has been observed in other years as well (WGNSSK, 2011).Sandeel stocks display strong interannual fluctuations, often biomass changes by more than a factor two between years.
In addition to this, North Sea sandeel biomass have historically experienced an apparent regime shift (Grandgeorge et al., 2008;van Deurs et al., 2009), going from a high average abundance level (∼2 Mt) to a lower average abundance level (∼1 Mt) around 1998-1999(WGNSSK, 2011)).The exact reason for this has not unambiguously been Introduction

Conclusions References
Tables Figures

Back Close
Full resolved, even though fishing pressure increased prior to the regime shift.However, other fish stocks predating on sandeels (Beaugrand et al., 2003) as well as the North Sea zooplankton composition (Beaugrand et al., 2002) experienced preceding regime shifts, pointing in the direction of trophic cascades.Sandeels are fished for fish meal and fish oil production.They bury in sandy sea bed when not feeding and therefore these stocks are localized to areas with sandy bottom substrate.These habitats have been mapped in detail from fishery log books (Jensen and Rolev, 2004), see Fig. 1.The temporal and spatial variability of the North Sea sandeel stocks makes this case study a rigorous test of the current capability of forecasting fish stocks.Such a variable resource poses a great challenge to an economical sector requiring a relatively constant supply to operate with low cost levels.Predicting stock dynamics on a short time scale of 1-2 yr would be of great value to fishermen and industry, allowing to anticipate inherent stock fluctuations and reallocate catch effort in time to other opportunities.Current fish stock assessment is data driven with limited spatial resolution (i.e.regional or basin scale).Marine Spatial Planing often require knowledge about ecosystem properties and responses at subregional scales (<50 km) and there is a pressing need to know how fish stocks respond to management actions and human activity on these smaller scales.
In this paper we present a spatial fish stock forecasting system addressing the needs above that can run in either hindcast/forecast mode or in reanalysis mode with flexible levels of data assimilations with the purpose of providing short term forecasts including the effects of seasonal and interannual variability; providing a flexible tool for assessing effects of alternative scenarios of stock management, involving both local scale and basin scale management actions; down scaling regional scale stock assessments to provide access to stock variability on spatial scales <50 km which is not accessible through regional scale stock assessments.Introduction

Conclusions References
Tables Figures

Back Close
Full 2 Forecasting model system The forecasting system presented here contains three elements as depicted in Fig. 2: at the bottom, ocean physics and biogeochemistry are provided by the POL-COMS+ERSEM model; at the top fish stocks and fisheries are described by the SPAM model.The top and the bottom are glued together by the SPAM model, describing early life stages of fish stocks.All three model components are described below.In this way all of the pathway from climate to fish landings is described by the system.We illustrate the forecasting model system on North Sea sandeel stocks.

Physics and biogeochemistry
The hydrodynamic and biogeochemical environment for the sandeel population is given by an implementation of the POLCOMS-ERSEM modelling system (Allen et al., 2001) implemented on the NW-European shelf (Butensch ön et al., 2012).This model is a further development of the MyOcean v0 model (Siddorn et al., 2007) for the NW-European seas that provides an enlarged domain corresponding to the v1 model (Edwards et al., 2012), allowing for a more realistic representation of the impacts of the shelf exchange processes on the shelf waters (Holt et al., 2012), and uses consolidated versions of the POLCOMS and ERSEM building blocks.POLCOMS is a primitive equation model for the coastal ocean using sigma-coordinates and a Arakawa B-grid with a piecewise parallel scheme for advection processes.Turbulence closure is achieved through the Mellor-Yamada model with an algebraic mixing length.Details on this part of the modelling system and extensive validation of its forecasting quality are given in the above mentioned references.From a model hindcast of the years 1960-2004 only the years 1970-2004 were used for this work to exclude spin-up effects.Figure 3 shows spatial and interannual temperature distribution to illustrate a major environmental impact factor on the sandeel habitats (see Fig. 1).The POLCOMS-ERSEM modelling system is coded in Fortran 90 to ensure high computational performance.Introduction

Conclusions References
Tables Figures

Back Close
Full

SLAM (Sandeel Larval Analysis Model)
SLAM (Sandeel Larval Analysis Model) is an individual-based model for describing transport and survival probability of early life stages (fish eggs and larvae) of sandeel (Christensen et al., 2008).Sandeels have pelagic larvae and their physical transport is described by passive advection-diffusion in this work.Larval hatching date is set to 20 February, based on biological observations (Wright and Bailey, 1996;Jensen, 2001).Larvae (try to) settle on a sandy habitat, when they reach the length of metamorphosis (L m = 40 mm).Larval growth is described by a temperature-controlled growth model (model 3 in Table 2 in Christensen et al., 2008).This (and alternative growth models) were compared extensively in previous work (Christensen et al., 2008) and results were found to be relatively insensitive to model details.The SLAM model addresses density independent growth and mortality processes.The model SPAM (see below) corrects for density dependence on larval growth and mortality (van Deurs et al., 2009).By hatching larvae from all suitable 585 habitat cells in Fig. 1, a 585 × 585 transport matrix T y is generated by the SLAM model for each year "y" describing habitat connectivity; a database T y was generated for y = 1970-2004 from POLCOMS-ERSEM output to span climatic variability.These transport matrices are the input for the stock model SPAM (see below) used to calculate spatial recruitment of new generations.The SLAM model is coded in modular Fortran 90 to ensure high computational performance.

SPAM (Sandeel Population Analysis Model)
SPAM (Sandeel Population Analysis Model) is a spatial stock model for settled sandeel that follows cohorts (generations) through their life cycle.The SPAM model was previously parameterized for a regional-scale setup (Christensen et al., 2009) whereas this work presents a high resolution setup on a 10 × 10 km grid, coincident with the physical grid used for the POLCOMS-ERSEM hindcast, so that spatial resolution of all three model components is alligned.Introduction

Conclusions References
Tables Figures

Back Close
Full The state variables in the SPAM framework are the abundance N y i ,a and average length L y i ,a in each habitat cell i in Fig. 1 for each generation with age "a" each year "y".The state variables describe conditions at 1 January each year, so they represent time snapshots and not averages over the year.State variables are updated from 1 January to 1 January by integrating processes over the elapsed year.The SPAM model and its parameterization were thoroughly discussed in (Christensen et al., 2009) and we adhere to the notation established there and provide a summary along with developments in Appendix A. The SPAM model can run as stand-alone life cycle model in hindcast mode or in reanalysis mode with various degrees of data assimilation.The integration step in the SPAM model is where Z = F + M + Z 0 is the total mortality composed of fishing pressure F , predation M and background (other) mortality Z 0 .δ is the Kronecker symbol.R is the recruitment which is the future number of (somewhere) settled offspring that a new juvenile can expect to produce in its lifetime, given that it settles in habitat cell i at year "y", and σ ξ i ,a is the probability of surviving to age "a" at year ξ, given settlement in habitat i .r y i is a function of fishing mortality F ξ>y i ,a (in addition to state matrices N ξ>y and L ξ>y ).In relation to stock management, it is interesting at what fishing pressure r = 1, i.e. when a newly settled juveniles closes its life cycle with one new recruit on average This defines a local maximum sustainable fishing pressure.For sinks, Eq. ( 6) will not have a solution (thus identifying sinks, using this definition).
The SPAM model is written as a Python class library.Simulating stock dynamics at the present spatial resolution takes less than half a second per year on a standard laptop, thus enabling extensive ensemble simulations at high spatial resolution.

Coupling models
It is important to realize that the three elements in the present forecasting system are fundamentally very different and address different temporal and spatial time scales: POLCOMS+ERSEM is an Eulerian model based on the continuum hypothesis solving the balances for mass, momentum and energy in form of partial differential equations, while SLAM is a Lagrangian framework and SPAM implements discrete difference equations of stock dynamics.Therefore at a technical and mathematical level the coupling of the three elements is by no means trivial and requires well-defined inter- pressure in POLCOMS+ERSEM, even though it may potentially be strongly spatially and temporally heterogeneous.

Stock data and data assimilation
Recently ICES Working Group on the Assessment of Demersal Stocks in the North Sea and Skagerrak (WGNSSK, 2011) moved from basin scale to regional scale assessment of sandeel stocks in ICES area IV, based on statistical analysis of catch landings, catch sampling and biological surveys.ICES WGNSSK adopted a regional division of the sandeel habitat network, see Fig.   to enable reanalysis runs.These WGNSSK fish stock assessments are considered as pseudo observations that can be assimilated independently in SPAM in reanalysis mode.Assimilation is performed multiplicatively to conserve sub-regional relative variability, i.e. if the average observation X A applies to area A, model state is rescaled as where i refers to a habitat cell.We will show later that within each of these regional assessment areas there is considerable spatial variability.

Future climatic variability
The computational cost of a full decadal ensemble forecast of the lower trophic levels of the marine environment that is sufficiently constrained to quantify all the uncertainties involved in such a system is currently untractable in an operational system.The Introduction

Conclusions References
Tables Figures

Back Close
Full question is then what to use as environmental forcing when future predictions of stock development must be assessed, i.e. what climate forcing to use computing future matrices T y needed for stock forecasting.We use a simple statistical extrapolation of the environmental forcing into the forecasting period, as this requires no explicit projection of the environmental state, that would come at extremely high computational cost without leading to significant improvement in the results of the sandeel prediction, leading to a computationally equilibrated forecasting system.Each statistical ensemble each year drew a random member of the set T y y∈ [1970;2004] to generate the response envelope to seasonal and interannula variability.

Model validation
To assess the skill of the framework at various levels of data assimilation we compared model output with results from statistical analysis based on fish landings and biological sampling data as reported by ICES stock assessment by WGNSSK (2011).The model skill for an output property X is quantified by the conventional cost function (Allen et al., 2007) cf where M X k and O X k is the kth out of N matching pairs of model output and observations on property.σ(O X ) is the observed RMS on property X .We have applied the absolute value of each model residual to accentuate a potential time bias or model drift, if present.The performance level is conventionally ranked as cf < 1: very good, cf < 2: good, cf < 3: reasonable and cf > 3: poor (Radach and Moll, 2006).In Table 1 we show the model skill for biomass (TSB) prediction for each WGNSSK assessment area for 1983-2004 from hindcast mode ("-" = no assimilation) to full reanalysis (assimilation of TSB).Similarly Table 2 shows the model skill for recruitment (R) prediction for each WGNSSK assessment area for 1983-2004 from hindcast mode ("-" = no assimilation) Introduction

Conclusions References
Tables Figures

Back Close
Full to full reanalysis (assimilation of R).In all cases, from no to moderate data assimilation, the model performs good or very good.Of course the last columns in Tables 1 and 2 are trivial, because the testing observations are assimilated in the model runs.
The fact that model performance in hindcast mode is quite decent reflects a fair mechanistic representation of biological processes in the coupled model framework, however for more operational purposes we recommend reanalysis mode to reduce model error further.

Results
In this section we provide some examples of model output that illustrate how the forecasting system is able to address typical questions in relation to the scientific basis for spatial fish stock management and marine spatial planning.
Each simulation starts with a spin-up period of 50 yr to relax age, size and spatial distribution of the fish stock.In the spin-up period we apply area resolved time averages for the period 1990-2011 of stock drivers F , M, Z 0 and T (the latter at 10 km resolution); we refer to this as reference conditions.After this may follow a reanalysis period, where observations (WGNSSK, 2011) are imposed as described in Sect.2.6.For all ensemble runs the same initial conditions (state matrices N and L) are applied to all ensembles, corresponding to system state at the end of the spin-up/reanalysis period.

Fish stock forecasts
In relation to fish stocks, short term forecasts means 1-5 yr, whereas long term forecasts means 10+ yr.After spin-up, a reanalysis run was performed for the years 1990-2004 with assimilation of T y , F y k and TSB y k ), after which the system state was copied to 100 identical ensembles.Forcings were set to reference conditions in the forecasting period, apart from T, which was assigned as described above.Introduction

Conclusions References
Tables Figures

Back Close
Full In Fig. 4 we show 20 yr forecasts for total stock biomass in WGNSSK areas A = 1-4 as defined in Fig. 1.Some features are very apparent: stock dynamics display a prominent two-year auto correlation, which is most clear in the forecast period 2004-, especially for area 1, but also clearly visible in the assimilation period 1990-2004, where predation and fishing pressure fluctuates between years and masks the two-year auto correlation.The two-year auto correlation effect on the stock biomass is observed for the stock (Arnott and Ruxton, 2002;van Deurs et al., 2009) and is a population density effect which also is an emergent feature in the SPAM model parameterization.This population density effect comes from food competition between adult fish and offspring or cannibalism (on offspring) within the stock.It is also seen that the two-year auto correlation is not in phase between areas or ensembles, reflecting weak larval exchange between major regions.An interesting property in relation to the validity of forecasts is the climatic envelope decorrelation time scale, i.e. the time it takes for interannual fluctuations to delete information of the stock past history, which sets an upper limit for medium/long-term predictions.This time scale is an intrinsic property of the stock dynamics model coupled to prevalent climatic forcing distribution T. We have previously shown that transport matrices T y of subsequent years have little correlation structure, if any (Christensen et al., 2008).If the stock starts with a certain biomass B 0 , the distribution of biomasses will evolve as P t (B) with where the attractor P ∞ (B) is the is distribution of predicted biomasses in the far future (if the stock is forced with same climatic variability observed presently for long time).
In Fig. 4 we see that the attractors P For the simplest reference system that can be solved analytically, diffusion with an attractive point B 0 , the estimator 10 gives the exact time scale in the analytic solution, where diffusivity D and advection speed k are parameters.
The fact that transport matrices T y of subsequent years can be considered independent makes τ insensitive to which year we start.If the estimator 10 is applied to data in Fig. 4, we find that areas 1-4 have τ = (6, 2, 3, 10+) yr respectively, and the total biomass have τ = 5 yr.Thus, unlike Lewy et al. (2004), we do not find that regional decorrelation time scales are consistently lower than basin scale decorrelation, and one may speculate that their conclusion is related to their statistical approach.After this period, the system returns to a dynamic attractor P

Spatial distributions and variability
An essential piece of information in spatial stock management is the spatial and temporal distribution and variability of the harvested biomass.In Fig. 7 we plot the biomass distribution (of 1 January) of contrasting years: 1996,1997,2004.The contrasting years are chosen so that 1996 and 1997 are before the apparent stock regime shift (1998-1999) and 2004 after the apparent stock regime shift.The coupled framework was run in reanalysis mode for the period 1983 to 2004, with assimilation of T, F and SSB.We see that the biomass gets more concentrated in the Dogger bank region after the apparent stock regime shift, whereas it previously was more evenly distributed over the habitat network in Fig. 1.The subsequent years 1996 and 1997 were chosen to elucidate potential spatial patterns in the two-year autocorrelation in the stock biomass discussed above in Sect.3.1; here we see that "blinking" is not in phase between areas.We also see that there is significant spatial variability within each main area in Fig. 1, stressing the importance of spatially downscaling stock assessment results.
In Fig. 6 we show spatial distribution of stock recruitment intensity, based on a reanalysis for the period 1983 to 2004, with assimilation of T, F and SSB. Figure 6a shows recruitment intensity plotted by source, i.e. which habitats contributes most to the overall stock recruitment, and Fig. 6b shows areas of successful initial settlement of juveniles in 2004 (regardless where recruits came from).We see that there is significant differences between habitats that contribute and receive recruits on small scales, i.e. the stock dynamics on small scale is strongly dependent on exchange of larvae between nearby habitats.In a previous work (Christensen et al., 2008(Christensen et al., , 2007) ) we demonstrated that typical exchange distances are ∼50-100 km.Such source/sink analysis allows to identify sensible areas (sources) and areas that can sustain heavy exploitation (sinks, which are automatically recolonized).Introduction

Conclusions References
Tables Figures

Back Close
Full

Local stock sensitivity
A way to access local stock sensitivity to harvesting is to plot and analyze the recruitper-recruit derived maximum local fishing pressure ψ i , Eq. ( 6).This is plotted in Fig. 8.The field ψ i was generated by solving Eq. ( 6) with same fishing pressure for 1+ age groups and otherwise reference conditions; r i in Eq. ( 4) was evaluated at N = 0 (and L(N = 0)) which gives the upper limit for ψ i , since mortality and growth processes are increasing and decreasing monotonously, respectively, with population density.Quite surprisingly, we see that some of the northern parts of the habitats are able to sustain highest local fishing pressure, but on the other hand northern parts also have variable values for ψ i : we see in Fig. 8 that quite a fraction of the north-eastern habitats have ψ i = 0 due to the fact that r i < 1 for any fishing pressure, suggesting that these areas are sinks.Another clear feature in Fig. 8 is that boundary habitats in the major areas have lower values for ψ i than inner habitats: this is simply because of higher hydrographic loss of offspring for boundary habitats (offspring have lower probability of getting advected to a suitable settling habitat).

Stock management scenarios
A central issue in the scientific basis for stock management is at what fishing pressure does the stock give maximal catch and at which fishing pressure does the stock (risk to) collapse.In Fig. 9 we show the predicted sustained catch for each WGNSSK area 1-4 along with the total catch, as function of spatial and temporal constant fishing pressure F for all 1+ age groups.The catches are 20 yr averages after a 50 yr relaxation period (to relax state variables, average out fluctuations and avoid stock depletion).The bars in Fig. 9 (only shown for the total catch) display the RMS on expected yearly catch that comes from inherent stock eigen fluctuations generated by population density effects, indicating that large uncertainty is associated with determining stock maximum sustainable yield due to stock eigen dynamics and climatic variability, and that optimal long-term harvest is inherently associated with large interannual variability in catch.Introduction

Conclusions References
Tables Figures

Back Close
Full Symbols in Fig. 9 show actual historic catch and fishing mortality (plainly averaged for the period 1990-2011) for each area (no fishing mortality were available for area 4).We see that catches are slightly overestimated by the SPAM model in hindcast mode; however, the figure also suggests that the sandeel stocks are currently exploited close to the maximum sustainable yield (MSY) point on average in all areas 1-3, even though catch and fishing mortality displays large fluctuations in the period 1990-2011.We see that different areas have different optimal yield fishing pressures, ranging from F ∼ 0.4 to 0.8 yr −1 , with F ∼ 0.7 yr −1 as North Sea wide optimum.Stock collapses occur gradually, at different fishing pressures for different areas, between F ∼ 0.8 yr −1 (area 4) and F ∼ 2.0 yr −1 (areas 1, 3).This upper limit is consistent with the local estimates based on recruit-per-recruit analysis, see Fig. 8.We see that the route of stock collapse is gradual when increasing fishing pressures: the stock concentrates in the most productive habitats, while extinguishing in lesser productive habitats; this pattern needs not be generic to other sedentary fish species.For F > F MSY catch fluctuations increase, until close to stock collapse, where fluctuations become very small: this is due to the fact that stock density effects disappear, because the stock does not utilize the habitat carrying capacity.Another way to appreciate this comes from the observation that fish stock models essentially map to a discrete logistic map (Murray, 2002): increasing F drives fish stock models reversely through bifurcation cascades, eventually collapsing the stock; the decrease in catch fluctuations for decreasing F comes along because catch is the product of F and biomass, which saturates at lower F .

Conclusions
We have presented an integrated forecasting system describing a single fish species, which is based on process-oriented simulation of ecosystem dynamics, from physics to fish.The forecasting system consisted of three major building blocks: the Eulerian POLCOMS-ERSEM system for physics and biogeochemistry, the Lagrangian model SLAM for early life stages and SPAM, the discrete box model for spatial fish Introduction

Conclusions References
Tables Figures

Back Close
Full population dynamics.All these three model systems are developed independently in different coding languages and are mathematically very different and address different time-scales.These aspects stresses the importance of well-defined interfaces, e.g.POLCOMS-ERSEM and SLAM is coupled by exchanging daily/subdaily time frames of full three-dimensional oceanographic data, whereas the coupling from SLAM to SPAM (and POLCOMS-ERSEM to SPAM) is via seasonal index matrices like T (and S and Q) produced by SLAM using oceanographic data from POLCOMS-ERSEM.
The three models can in principle run in online mode, exchanging data via files (or data pointers), or offline, exchanging data via files only.Times series for full threedimensional oceanographic data are relatively big (0.1-10 TB) for 10-50 yr needed for spanning climatic variability, but this becomes feasible with present day storing capacity and data exchange rates.Seasonal index matrices, like T, require negligible storage once computed.The offline mode is definitely of advantage during interactive development, iterative calibration runs and for the statistical extrapolation of past states into the future applied in this study.
Online coupling may be relevant for production runs, but, as mentioned above, the three model components are not load balanced.However, online mode would become necessary to include the grazing feedback of the fish/fish larvae on the plankton biomass of the lower trophic level model.The relative model load is set by the integration time step, which is smallest for POLCOMS-ERSEM (in the order of seconds for the external mode and around 15 min for the internal mode), in the middle for SLAM (30 min) and largest for SPAM (∼seasons-1 yr).Consequently, it becomes attractive for iterative fish stock simulations (e.g.parameter calibration runs over same time window) to run offline, whereas the addition of the SLAM/SPAM component has neglectable impact on the POLCOMS-ERSEM system in terms of computational resources, unless the SPAM simulation includes an excessive number of ensembles.
We used the fish species sandeel as case study, however there is nothing that prevents application of the framework to other species, either isolated (i.e. a single species model) or many species together (i.e. a multi species model).The most interesting Introduction

Conclusions References
Tables Figures

Back Close
Full species that could be covered with little reparameterization is other small pelagics like sprat, anchovies and sardines, or other sedentary species like flat fish.The major issue in relation to multi species models for fish is the complexity of species interaction, both due to lack of data and uncertainty about many aspects of species interaction, e.g.switching and preference in predator-prey relations, which is further complicated by density-dependence.When parameterizing the fish stock model component SPAM, we have not imposed a stable equilibrium state, but rather included the limit cycle (attractor) of the biological model which is consistent with biological observations, displaying negative autocorrelation of biomass between subsequent years.This of course poses extra challenges when parameterizing a model, but this has been accomplished by focusing on statistical properties of the limit cycle of the biological model, as elaborated in Appendix A.
Facing that the computational cost of a full decadal ensemble forecast of the lower trophic levels of the marine environment requires excessive computer resources, we have applied a simple "pick-a-random-past-year" scheme for future climate forcing that reproduces statistical properties of a historical time series by construction; however alternative schemes to generate synthetic T time series are conceivable, e.g.convex interpolation in a historic time series like where s z is a random vector with z s z = 1.This sampling scheme has smaller variance < T * T > − T * T than the random-year sampling scheme applied in this paper.Also, this scheme may blur a potential covariation and clustering in the historical set T y .Alternatively principal-component correlation with teleconnection indices (Barnston and Livezey, 1987), like NAO, could be applied.We believe that a scheme for generating a synthetic T time series should satisfy certain minimal requirements like T > 0 and reproducing basic statistical properties of a T y .We found that the decorrelation time τ of biomass dynamics due to stochastic climatic variability is of order 2-6 yr, with some variations between areas in the modelled region; this is the time it takes for random

OSD Introduction Conclusions References
Tables Figures

Back Close
Full kicks from the climatic variability to wash out the past biological state of the fish stock; this time frame is somewhat encouraging in that it allows for some short to medium time scale economical planning.When the time horizon exceeds ∼6 yr, the principal model output is the asymptotic probability distributions (attractors) of the system properties, e.g.P ∞ (B) the probability distribution of biomasses in the far future.These attractors are also very useful for management scenarios, they limit the ignorance of the future and provide an average trend as well as variability on the stock state, when subjected to alternative management scenarios.
The complexity of life processes increases as one moves up the trophic chain and therefore many assumptions have to be made for early life stages (the SLAM model) and adult fish (the SPAM model).In the SLAM model, e.g.active motion of larvae has not been included, as well as grazing feedback on zooplankton from sandeel larvae.The latter is partially justified by the fact that sandeel larvae and adults on average only constitute a fraction of the zooplankton grazing biomass, even though locally sandeel may dominate the zooplankton grazing (A. Rindorf, personal communication, 2011).This isolated grazing feedback from sandeel on zooplankton is partially represented by the density dependence of S, see Eq. (A10).Anyway, all these uncertainties are absorbed in the product of T and S in the integrated framework.In the SPAM model the major short comings are the lack of spatial and temporal heterogeneity of predation pressure on the fish stock.This will decrease the decorrelation time scales τ found in Sect.3.1.In the present paper, a spatial and temporal average predation pressure was applied, as available from stock assessment work (WGNSSK, 2011).In the future, it is conceivable that temporal variability can be modelled from multi species stock assessment time series, in conjunction with predator-prey interaction matrices.Additionally, the habitat carrying capacity did not include temporal variability, due to lack of data.(Jensen et al., 2011), triggered by daily feeding cycles.In Appendix A we show for high resolution setups like the one in the present paper that this will create a minor spill-over effects between adjacent habitats, but will not affect results on mesoscale, ∼50 km.Even though fish stock assessment results are associated with significant uncertainty model skill benchmarks in this paper indicate that numerical representation of biological oceanographic processes underlying recruitment still needs to be improved, even though skill currently is quite good.
We have presented a few illustrative examples of using our forecasting framework on questions relevant to management of fish stock, e.g.providing short term forecasts including the effects of climatic variability and down scaling regional scale stock assessments.We have demonstrated that the integrated forecasting system is a flexible tool for assessing effects of alternative stock management scenarios and providing biological insight into stock dynamics.The spatial explicit nature of the forecasting system allows for very rich and highly policy relevant output, e.g.warning signs of over fishing and imminent stock collapse from spatial signatures of how habitats are depleted, as well as risk of stochastic extinction of local fish stocks.

Appendix A SPAM parameterization
In this Appendix we briefly summarize the parameterization of the biological processes in the SPAM model.Most parameters and process representations of the high resolution (10 km) setup are the same as in the previous regional scale (∼100 km) setup, see (Christensen et al., 2009), where a more elaborate derivation and model discussion can be found, along with a parameter sensitivity test.We adhere to the notation established there.In the following i designates habitat cell in question, "y" year for time varying properties, and "a" age of a cohort.Introduction

Conclusions References
Tables Figures

Back Close
Full The larval conditional (given transport) survival S is the product of the predation avoidance S p and starvation avoidance S s .The former is parameterized as where t larv is the duration of the larval drift phase and M pl the average predation risk, computed from size spectrum theory as where τ is the expected life time of a newly hatched larvae before it get eaten on average.The prefactor η is slightly species specific and depends on average growth speed and length-at-settlement. ξ i is an area specific weighting factor of order unity (explained below).
which for sandeel gives η ∼ 0.403, where −ν ∼ −3/4 is the standard length-mortality scaling (West et al., 2001).Starvation avoidance S s is parameterized as Finally, stock catches are estimated from the narrow-season limit as where t c is duration from the end of hibernation to the middle of the (narrow) catch season.w

A4 Reference conditions
Reference conditions refer to average values (WGNSSK, 2011) (over the reference period 1990-2011) applied of stock drivers: predation mortality Mi,a and fishing mortality Fi,a .For predation mortality, only a North Sea average value were provided (WGNSSK, 2011), Ma , and the area specific predation mortality was constructed as so that larval and adult predation risk are correlated (reflecting spatial heterogeneity in density of predators).Similarly, under reference conditions, T was averaged over the reference period 1990-2011 T = < T y > y∈ [1990−2011] ; reference carrying capacity was set as outlined below.

A5 Inferred biotic/abiotic parameters
For the remaining model parameters, there were insufficient empirical basis for direct parameter estimation, and these parameters are estimated from SPAM model output by minimizing a residual of remaining observations.These parameters are mostly area weighting factor representing spatial heterogeneity and these are given in Table 4.No spatial heterogeneity were allowed for α and τ to avoid over-parameterization.We also see that spatial variability of C i and ξ are in a sensible regime, creating confidence to the robustness of the parameterization.The residual used to infer (α, τ, C i , ξ) were the sum of residuals for <TSB>, RMS (TSB) and average growth for the reference period 1990-2011 in each area as provided by (WGNSSK, 2011).Introduction

Conclusions References
Tables Figures

Back Close
Full  Full  Full  Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

(
of a new generation) which is assembled in SPAM from the transport matrix T passed from SLAM and fecundity Q and conditional survival chance S, both of which depends on N and L. The primary forcings of the model are (F y i ,a , M y i ,a , T y i j ) The indirect forcings are the habitat carrying capacity C y i , see Appendix A for model parameterization details, where also the length-specific growth parameterization g i is described.A useful diagnostic quantity that we will discuss later is the recruit-per-recruit number r Discussion Paper | Discussion Paper | Discussion Paper | faces.Data exchange between POLCOMS+ERSEM and the SLAM model occurs in daily 3-D time frames of average ocean fields; the SLAM-SPAM interface exchanges yearly transport matrices T y and S y .There is currently no dynamic feedback (zooplankton grazing) from SLAM+SPAM on POLCOMS+ERSEM and this makes sense, because sandeel stocks only constitute part of the planktivorous biomass in the North Sea.The zooplankton grazing feedback is currently applied as an average grazing Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

∞A
(B) are different for each area A, where area 2 has the broadest range.Also we note that area 4 is predicted to be prone to years of low biomass abundance, due to climatic variability.The climatic envelope decorrelation time scale τ is then the speed at which P t (B) converges to P ∞ (B).A rough estimator for τ is

∞A
(B) characteristic for each area.The asymptotic total biomass attractor for North Sea sandeel stocks P ∞ NS (B) predicted by the SPAM model is shown in Fig. 5.It is important not to confuse the climatic envelope decorrelation time scale τ with climatic change time scales; the former addresses how the stock responds to climate variability, i.e. interannual climatic changes (considered stochastic) on short time scales not related to global warming (or other long-term time-scales related to the Earth system).Finally we notice that the forecasts 2004 -predict stock returns to a higher average level, consistent with recent years observationsDiscussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | The current parameterization however did not indicate very strong gradients in habitat carrying capacity in the North Sea for sandeel (see Appendix A), but spatially coherent interannual fluctuations are quite likely.Finally post settlement adult migration was not included; adult sandeel are believed to display only small scale migration ( 28 km) Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | y i ,a (t c ) is the average weight at catch, computed from Eqs. (A3) and (A2).Migration between adjacent habitats is not included in the present setup, as it is believed to just imply some local smearing of population densities for sandeel.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | mate Change on North Atlantic Fish Stocks, Bergen, Norway, 11-14 May 2004, 2005.1439 Moellmann, C., Diekmann, R., Muller-Karulis, B., Kornilovs, G., Plikshs, M., and Axe, P.: Reorganization of a large marine ecosystem due to atmospheric and anthropogenic pressure: a discontinuous regime shift in the Central Baltic Sea, Glob.Change.Biol., 15, 1377-1393, doi:10.1111/j.1365-2486.2008.01814.x,2009.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Table 2 .
SPAM model cost function "cf" (Eq.8) on recruitment prediction R y A = i ∈A,j R T, F , SSB T, F , TSB T, F , SSB, R

Table 4 .
Inferred biotic/abiotic parameters in the SPAM model obtained from residual minimization.