Variability of plankton production during the spring bloom in NW Iberia

: A coupled physical (Regional Ocean Modeling System, ROMS)−biogeochemical model (N 2 PZD 2 ) at a horizontal resolution of 3.5 km was implemented for N and NW Iberia, an area of high productivity associated with upwelling. The physical model has been the object of previous studies and has proven its capability to perform well in reproducing the main oceano-graphic features in the area (upwelling, river plumes, slope currents, fronts, filaments), which is fundamental to properly representing the variability and distribution of the biogeochemical variables. The biogeochemical model was set up to account for the main nutrient inputs in the area: upwelling and rivers. Upwelling input required proper characterization of the nutrient content variability of the Eastern North Atlantic Central Water, which was achieved by using a tempera-ture−nitrate relationship obtained from observations to impose nitrate at the open boundaries. The resulting biophysical model accurately reproduced the timing and interannual variability of the spring bloom compared with satellite chlorophyll (chl a) observations. A comparison with the Instituto Español de Oceanografía’s in situ spring-monitoring Pelacus cruises (which include plankton) revealed that the model was able to reproduce the variability at shorter scales (days) and demonstrated its ability to complement the observational data and reveal the variability in the area around the spring transition. In this respect, both the model and observations showed that productivity on this narrow shelf is affected by seasonal upwelling that results from the interplay of wind, river plumes and light intensity, all varying at interannual, seasonal and event scales.


INTRODUCTION
Primary production in the Iberian Atlantic shelves, located north of the Canary-Portugal eastern boundary ecosystem, is largely influenced by seasonal upwelling. Two main seasons, driven by the prevailing wind regime, are distinguished in the area: an upwelling season from April to September and a downwelling season during the rest of the year (Wooster et al. 1976). Interannual variability exists in the time of transition between seasons and their intensity (Ruiz-Villarreal et al. 2006). In addition, a significant part of the variability is concentrated in short-term wind events (3−14 d), meaning that up welling pulses occur in the downwelling season and downwelling pulses and relaxation occur in the up welling season (Blanton et al. 1984). In the spring transition from downwellingto upwelling-favourable winds, the input of nutrients promotes plankton proliferation (Varela 1992, Casas et al. 1997. The associated high levels of primary production sustain a wealthy fisheries and aquaculture industry that represents an important percentage of the economy in the area.
Since 1986, the Instituto Español de Oceanografía (IEO) has carried out a series of annual spring surveys (Pelacus cruises) covering the shelf from Porto (Portugal) to the southeastern Bay of Biscay. The objective is to monitor the pelagic ecosystem from hydrography to fish (Bode et al. 1996, Bernal et al. 2007, Santos et al. 2013. The monitoring of environmental conditions around the spring transition performed during the Pelacus cruises has shown intense variability in the volume of freshwater on the shelf, the intensity of the Iberian Poleward Current (IPC) and its penetration into the Cantabrian Sea (Bode et al. 2002, Ruiz-Villarreal et al. 2006) and the frequency and intensity of upwelling pulses around the spring transition. These processes affect oceanographic conditions like temperature, frontal structures and stratification, and hence plankton productivity, which has been observed to vary at different temporal and spatial scales (Bode et al. 1996, Gil et al. 2002, Ruiz-Villarreal et al. 2006, Picado et al. 2016. The Pelacus surveys provide information on the variability of plankton, but interpretation of the observations is complex because sampling is not synoptic (i.e. different areas are sampled on different dates) and there is strong variability in short-term events. Satellite imagery provides synoptic information on the distribution of phytoplankton/chlorophyll (chl a), but this data is restricted to the sea surface on clear days. In addition, the estimation of chlorophyll from ocean colour satellite images in coastal waters is challenging because in addition to phytoplankton, total suspended matter and/or coloured dissolved organic matter also affect the optical properties of water (IOCCG 2000). These concerns regarding temporal and spatial coverage of phytoplankton observations could be solved by the use of adequate biogeochemical models in the sense that they provide time-varying, 3-dimensional (3D) information on several biogeochemical variables, including plankton. However, it is still a challenge for state-of-the-art regional biogeochemical models to reproduce the variability of primary production in coastal areas. Most of the difficulty arises from the strong influence of hydrodynamics on the ecosystem, which requires the use of a reliable physical model that adequately represents the variability in oceanographic conditions (Doney 1999, Anderson 2005, Skogen & Moll 2005. The complexity of the biogeochemical model may also play a role, with the optimal complexity of models being a focus of research (Friedrichs et al. 2007, Ward et al. 2013, Kwiatkowski et al. 2014, Xiao & Friedrichs 2014. In this contribution, we describe the implementation of a coupled hydrodynamical−biogeochemical model configuration for the Iberian Atlantic shelf and slope. The biogeochemical model (N 2 PZD 2 ) is coupled to an existing configuration of the Regional Ocean Modeling System (ROMS) hydrodynamical model for the same area, which was developed by the modelling group of the IEO over the last years to simulate hydrodynamical variability at scales relevant for assessing the effect of the environment on the ecosystem. The physical model configuration focuses on high-resolution shelf and slope processes (e.g. upwelling, river plumes, slope currents) and adequately accounts for the variability of these processes in response to wind events and tidal variability , Marta-Almeida et al. 2013.
The objective of this study is to combine the coupled hydrodynamical−biogeochemical model results with observations (in situ and satellite), combining the strengths of each method to better characterize and understand the oceanographic conditions that (1) trigger the phytoplankton blooms on and off the shelf and (2) affect the distribution of plankton around the spring bloom at the proper scales of variability (e.g. tides, meteorological events, seasonal, interannual). We examine with special emphasis the interaction of upwelling/downwelling conditions on the shelf with variable river runoff and nutrient input, using the model to gain insights into the influence of river plumes on the variability of phytoplankton biomass.

Physical model ROMS
We used a configuration for West and North Iberia (see the configuration summary in Fig. 1) of the 3D, free-surface, hydrostatic, primitive equation ocean model ROMS (Shchepetkin & McWilliams 2005) (Rutgers version 3.5; www.myroms.org) for the hydrodynamic calculations. ROMS uses stretched terrainfollowing coordinates in the vertical, also known as s-coordinates, and curvilinear coordinates in the horizontal. Details of the configuration regarding the horizontal and vertical discretization of the computational domain, grid smoothing, numerical methods used and the turbulent closure are summarized in Table 1.

Biogeochemical model
The nitrogen-based Fasham-type (Fasham et al. 1990) biogeochemical model in the version implemented in ROMS by Fennel et al. (2006) was used for our ecological modelling component. This N 2 PZD 2 model comprises 7 state variables: one group of phytoplankton, one group of zooplankton, 2 nutrients (ni-trate [NO 3 ] and ammonium [NH 4 ]), 2 groups of detritus (small detritus, which is the detrital pool smaller than 10 μm, and large detritus, which results from the aggregation of small detritus and phytoplankton) and phytoplankton chlorophyll. The first 6 variables are measured in units of nitrogen (mmol NO 3 m −3 ) and fulfil mass-conservation equations (they take part in the nitrogen cycle); the last variable is not coupled to the others. Chlorophyll is explicitly calculated in the model from the phytoplankton biomass by conversion into chlorophyll units (mg chlorophyll m −3 ; Fennel et al. 2006), which implies multiplication with the ratio of chlorophyll to phytoplankton biomass (which is not linear and depends on the photo-acclimation of the phytoplankton cells) and the assumption that only a fraction of phytoplankton growth is used for the synthesis of chlorophyll.
The temporal and spatial evolution of the biological tracers is modelled by means of advection-diffusion-reaction ordinary differential equations (ODEs); a detailed description of each ODE for the Fennel model can be found in Fennel et al. (2006). The model simulates processes such as phytoplankton growth as a function of photosynthetically available radiation (PAR), temperature and nutrient availability (NO 3 and NH 4 ), mortality by zooplankton grazing and other linear processes. In this model, dead phytoplankton and zooplankton are transferred to the detritus pool, which is remineralized to NH 4 . The nitrification of NH 4 to NO 3 is also considered as well as the settling of phytoplankton and small and large detritus, which is responsible for the accumulation of organic matter at the sediment−water interface. The equations that describe the biological processes contain multiple parameters on which the final solution is strongly dependent. The values of the parameters were selected for simulating the shelf ecosystem of the Middle Atlantic Bight and are the default values in the ROMS distribution (see Table 1 in Fennel et al. 2006). In this study, we modified some of the parameters to values suitable for representing our ecosystem at the beginning of the spring bloom (Table 2); namely: (1) The half-saturation concentration for the uptake of NO 3 (K NO3 ) was lowered with respect to the default value to reproduce the conditions at the beginning of the spring bloom in NW Iberia. At this time, the diatom Chaetoceros socialis is dominant in the area (Bode & Varela 1998a,b). In a field cruise in the ría de Vigo at the end of June, Seeyave et al. (2013) obtained a value of K NO3 = 0.37 mmol NO 3 m -3 in a mixed diatom bloom. Taking into account that the spring bloom typically starts at the beginning of April, faster consumption of NO 3 would be expected compared to June. For this reason, we lowered the value of K NO3 to 0.25, assuming that the nutrient uptake (function L NO3 in Eq. 3 of Fennel et al. 2006) would increase for the same concentration of NO 3 . This value is within the ranges reported in the literature: 0.007−1.5 in Fennel et al. (2006), 0.02−10.2 for diatoms in Sarthou et al. (2005) or 0.1−0.7 for cultured oceanic phytoplankton in Eppley et al. (1969).
(2) The value of the initial slope of the P−I curve (α) was lowered with respect to the default. This param-  (Beckmann & Haidvogel 1993) 0.2 Haney number rx1 (Haney 1991) Maximum of 5  eter depends on the phytoplankton species composition of the bloom and varies with time and space, even for the same species. It reaches the highest values in early spring, which, according to a literature review in the area, range from approximately 0.02−0.07 mol C g chl −1 (W m −2 ) −1 d −1 in the first 50 m of the water column (Bode & Varela 1998a,b). A value in the middle of this range (α = 0.05) was selected for this study.
(3) The maximum grazing rate for zooplankton (g max ) was increased with respect to the default in the model (0.6 d −1 ); the selected value was an approximate average (1 d −1 ) of the values published for large and small zooplankton (0.96 and 1.2 d −1 respectively) in the N 2 P 2 Z 2 D 2 biological model configuration of Koné et al. (2005) for simulating the southern Benguela upwelling ecosystem.
The ODEs that describe the evolution of the state variables of the biogeochemical model are solved by using a Euler implicit scheme, with the exception of settling, which is treated explicitly. Implicit methods guarantee unconditional stability, which is a desired property for complex biological problems. Implicit schemes also restrict biological losses to positive concentrations. In our configuration, we ensured that all consumption terms were treated implicitly, including phytoplankton mortality and zooplankton excretion, which were treated explicitly in the default ROMS code. We switched off processes (i.e. fluxes were set to zero) when concentrations of species were below a reasonable threshold (which was already done for zooplankton and phytoplankton). This treatment reduced negative values in the biogeochemical fields and allowed us to avoid the use of a positive definite algorithm such as the Multidimensional Positive Definite Advection Transport Algorithm (MPDATA) for tracer advection.
Details of how the system of ODEs is solved numerically and some modifications carried out to the ROMS subroutine 'fennel.h' that contributed to guaranteeing positive definite solutions, mass conservation and numerical stability are provided in Text S1 in the Supplement at www.int-res.com/articles/suppl/ m708p045_supp.pdf.

Forcing and boundary conditions
For the hydrodynamic model ROMS, realistic atmospheric forcing was obtained from 2 different limited area operational models: either MM5/WRF from the regional agency Meteogalicia (www.meteogalicia.es), with a time resolution of 1 h and a spatial resolution of 30 km, or HIRLAM from the national agency AEMET (www.aemet.es), with a temporal resolution of 6 h and a spatial resolution 0.2°. High temporal resolution atmospheric forcing allows us to resolve the time scale of wind events that exert a strong influence on the dynamics in the area.
To properly reproduce the river plume dynamics and their effect on the biological components, the input of 26 rivers was accounted for (see Fig. 1) We highlight the Douro and Miño for being the most relevant in the formation of the Western Iberia Buoyant Plume, together with other minor Portuguese and Galician rivers : the Nalón, Sella, Deva, Pas, Nervión and Adour, which are the main contributors to the Cantabrian Sea river plumes (González-Nuevo & Nogueira 2014, Fernández-Nóvoa et al. 2019) and the Loire and Gironde, forming river plumes at the French continental shelf (Lazure & Jegou 1998). The methodology for constructing the time series of river runoff was published in Otero et al. (2010), although several improvements and updates have been carried out since the publication of that paper, mainly due to the availability of new data. Gaps in real data measured at the gauge station closest to the river mouth are filled with climatological values or estimations based on contiguous rivers. The sources of river runoff data were the Portuguese Agency for the Environment (www.apambiente.pt) for Portuguese rivers, the Confederación Hidrografica Miño-Sil for Miño river, the Galician regional agency for water (https:// augasdegalicia.xunta.gal/) for Galician rivers, the Confederación Hidrográfica del Norte (now Confederación Hidrográfica del Cantábrico; www. chcantabrico.es/) for Cantabrian rivers and Eau France (http://hydro.eaufrance.fr) for French rivers. A constant temperature of 12°C and constant salinity of 5 was imposed for each river.
Hydrodynamic open boundaries and initial conditions were obtained from a large-scale model, MER-CATOR PSY2V2, which was the configuration of the MERCATOR Ocean operational system in the North Atlantic for the years of our simulation period: 2006 and 2007 (a description of PSY2V1 and V2 can be found in Lellouche et al. 2006 and of later versions PSY2V3 and V4 in Lellouche et al. 2013). PSY2V2 is an operational forecasting system for the North Atlantic that assimilates sea level anomalies, sea surface temperature (SST) and temperature and salinity from in situ profiles. Horizontal velocity, salinity, temperature and sea surface height fields are provided at a horizontal resolution of 1/6° and a temporal resolution of 1 d, although the model resolution is 1/12°. These fields were imposed at the model boundaries by following the method described in Marchesiello et al. (2001), which combines radiation and nudging boundary conditions. For the barotropic mode, the radiation of waves that are generated inside the domain is allowed, with a weak nudging to the boundary values prescribed if the flow is directed outwards (TNUDG = 360 d). For the incoming flow, strong nudging is prescribed, and a relaxation time scale of 0.4 d (OBCFAC = 900) is selected. Since tides are considered (we use the TPX06.2 global inverse tide model), Flather conditions are also used for the barotropic mode. A sponge layer (20 grid cells) of increased horizontal viscosity (300 m 2 s −1 ) is set close to the boundaries and is zero elsewhere in the domain.
Open boundaries and initial conditions for the biogeochemical model were prescribed considering that NO 3 is the limiting nutrient in the area (Álvarez-Salgado et al. 2002, Howarth & Marino 2006. There are 2 main sources of nutrients coming into our system: upwelling and rivers. Nutrient content of upwelled waters is related to the nutrient content of the Eastern North Atlantic Central Water (ENACW; Pérez et al. 1995, Pollard et al. 1996, Bode et al. 2002, Álvarez-Salgado et al. 2002. ENACW in the area has 2 components: subpolar, driven by winter convective mixing in the NE Atlantic (north of our area) and subtropical, coming from the south (Bode et al. 2002). This water mass experiences interannual variability, and in the case of subpolar ENACW, its nutrient content mostly depends on winter convective mixing (Somavilla et al. 2009, Hartman et al. 2010. Being able to introduce this variability into the model would allow us to reproduce the interannual variability in primary production. This would not be possible if we used a NO 3 climatology like the monthly World Ocean Atlas (WOA) as in Dabrowski et al. (2014). Therefore, we applied polynomial fitting techniques to in situ data to obtain NO 3 concentrations from temperature, similarly to Fennel et al. (2006). With this relationship, we can 'translate' the MERCATOR PSY2V2 temperature fields, which capture variability in winter convection and advection, into NO 3 fields. Therefore, the description of the spatial and temporal variability of nutrient supply at the boundaries is substantially improved with respect to climatology.
The temperature−NO 3 relationship was computed using the IEO Radiales Profundas database (IEO-RADPROF). RADPROF is an IEO core project aimed at studying climate variability in N and NW Iberia (http://radprof.ieo.csic.es/). RADPROF cruises have been carried out every year since 2003 (Ruiz-Villarreal et al. 2006, Prieto et al. 2013. From 2003 to 2010, 3 hydrological standard sections were occupied twice a year (around winter and summer): Cape Finisterre (43°N), Cape Ortegal (8°W) and Santander (3°47' W). Measurements in RADPROF cruises include CTD, O 2 , nutrients and currents from vesselmounted ADCP and LADCP at regularly distributed stations located at a short distance over the slope and 15 and 20 miles over the oceanic planes. Additionally, measurements of chlorophyll a (chl a), phytoplankton and zooplankton are performed. For this paper, RADPROF data from the 2004 (February and September), 2005 (January and August), 2006 (July), 2007 (January) and 2008 (February and September) cruises were used. Data from the National Oceanographic Data Center (NODC) WOA2009 was also tested, but it was found that the RADPROF data set was a better fit to the ENACW properties as described by Ál varez-Salgado et al. (2002), especially at low temperatures; additionally, data are available below 500 m depths. This is illustrated in Fig. 2, which shows the WOA2009 data (blue dots) and the corresponding polynomial fitting (yellow line), together with the IEO-RADPROF data (red dots) and its associated polynomial fitting (black line). The linear regression analysis of NO 3 and temperature published in Álvarez-Salgado et al. (2002) is included in Fig. 2 (cyan line), and it is based on the characterization of NO 3 content in the Iberian margin made by the authors (which ranges from 2.6 ± 1.3 μmol kg −1 NO 3 in the branch of ENACW of subtropical origin [around 14°C] to 9.4 ± 1.2 μmol kg −1 in the subpolar ENACW water [approximately 12°C]). The better fitting of the IEO-RADPROF to the classical description of ENACW in Álvarez-Salgado et al. (2002) is clear in the figure. The way in which the above-described approach allowed us to consider the interannual variability in nutrient supply through the model boundaries for 2006 and 2007 is described in Text S2 and illustrated in Fig. S1.
Nutrient load by rivers is another source of nutrients to the shelf, and it must be included in the model to account for its role in fuelling production. We performed a review of available databases and literature reports of nutrient content in rivers. As an initial approach, we decided to impose a constant value of NO 3 concentration to all the rivers of 10 mmol m −3 , which is within the lower range of concentrations we found during spring. We reviewed NO 3 values for the Galician rivers in the yearly reports on the oceanographic conditions for Galicia from INTECMAR monitoring (www.intecmar.gal/). For the Douro in Portugal, our value is higher than that reported by Magalhaes et al. (2008), who measured a maximum of 0.5 mg l −1 in April (8 mmol m −3 ), and slightly lower than the values in the data set from the Portuguese National Water Resources Information System (https://snirh.apambiente.pt/), where values of around 15.5 mmol m −3 were obtained by averaging the measurements of several stations located at the lower basin of the estuary. Nutrient load for the rest of the state variables was considered to be negligible.

Simulation set
The model was run for the years 2006 and 2007. The time step for hydrodynamic calculations was 180 s for the baroclinic mode and 4.5 s for the barotropic mode. No spin-up period was considered to be necessary for the following reasons: (1) the model was initiated in winter when the water column is well mixed. The initial and boundary conditions were extracted from the same large-scale modelwhich is a realistic configuration with data assimilation, meaning that the initial condition was close to the model dynamical equilibrium (e.g. in contrast to a climatology) -whose original horizontal resolution and physics are similar to our model. This meant that there were no initial instabilities in the results as shown by the domain-averaged eddy kinetic energy (see Fig. S2); (2) For the biogeochemical model, we took into account that primary production in the area is largely driven by the wind regime, which favours the upwelling of subsurface water (ENACW in our area; Álvarez-Salgado et al. 2002), as is typical in eastern boundary upwelling ecosystems (Aristegui et al. 2006). We also considered the additional source of nutrients from river plumes, which is an important contribution of phytoplankton production in the area during the spring bloom (Bode et al. 2019). Both processes are physically driven, so we assume that if the model is in dynamic equilibrium and the nutrient inputs from the ENACW and rivers are well reproduced, then no long spin up is necessary for the biogeochemistry component either. We recall that NO 3 was obtained from the temperature fields used for the initial and boundary conditions of the hydrodynamic model. This implies that no inconsistencies between the initial and boundary conditions existed and that the resolution of the boundary forcing was similar to that of our model, which minimized the instabilities. In fact, the model initialization in winter implied starting conditions with (almost) zero plankton, which is realistic for the area in this season. No spin up also implies that the potential sources of nutrients from in situ organic matter regeneration were considered to be of less relevance for the spring bloom, which is consistent with the fact that the contribution of remineralization to the nutrient content in shelf waters is expected to be maximum in autumn as a consequence of the decomposition of organic matter produced in spring . To further investigate the effect of the spin up, additional tests were carried out, modifying the ini-

Remote sensing data
Daily images of SST and chl a obtained by the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument on NASA's Aqua satellite were used in the present study covering the area of the model domain (see Fig. 1). All images were processed at a spatial resolution of 500 m by the Remote Sensing Group of the Plymouth Marine Laboratory in the UK. The chl a images were obtained by means of the OC3M algorithm (O'Reilly et al. 2000).
Each pixel of these images was assigned a quality flag (L2-flags; http://oceancolor.gsfc.nasa.gov/). In SST images, we masked out pixels with SST warning or failure flags (28 and 29, respectively). In the case of the chl a images, pixels flagged as chlorophyllwarning, failure or turbid were masked out (flags 22, 16 and 12, respectively).

Pelacus cruises
The Pelacus cruises are core surveys of the IEO which have been carried out every spring -the spawning season of most pelagic fish in the areasince 1986. The main objective of this survey, which covers the continental shelf between 30 and 200 m depth from the Portuguese border to the French border along parallel transects, is the acoustic determination of biomass of pelagic fish stocks and the distribution of pelagic fisheries. During these pelagic ecosystem surveys, samples are obtained from different compartments: hydrography, plankton, ichthyo plankton, seabird and mammals presence (using trained observers), collection of information related to human pressures (i.e. presence, type and abundance of fishing vessels, marine debris) and, more recently, marine litter sampling (Bode et al. 2002, Bernal et al. 2007, Santos et al. 2013, Gago et al. 2015. For this paper, we used data from the Pelacus surveys in 2006 and 2007; in particular, temperature and salinity (see Fig. 5) and nutrients, chlorophyll and plankton (see Fig. 8). It is relevant to mention that the plankton data correspond to depth-integrated concentrations (up to 100 m) since they are obtained using 20 μm mesh Bongo nets and vertical hauls from 100 m depth to the surface, which needs to be considered for a homogeneous comparison with model results. Common units are also necessary to compare the model with the observations. In the particular case of zooplankton, the modelled concentrations are given in mmol NO 3 m −3 ; these were transformed into mg C m −3 (the observation units) by applying a C:N relationship of 8, which is slightly higher than the one reported by Bode et al. (2020).

Seasonal variability in 2006 and 2007
The weekly domain-averaged concentration of satellite and modelled surface chlorophyll for 2006 and 2007 is shown in Fig. 3 for the whole model area (Fig. 3A), the Atlantic oceanic part (computational area to the west of 8°W in Fig. 1; results in Fig. 3B) and the Cantabrian Sea (computational area to the east of 8° W in Fig. 1; results in Fig. 3C). The satellite domain average was obtained using daily MODIS chl a images, from which the weekly averages were calculated. The equivalent model domain average was calculated using the outputs closest in time to the available satellite images and applying the same flags. Results from 2 model runs with atmospheric forcing from 2 different meteorological models are plotted in the figure tic for the open ocean in the study area (see, for instance, Bode et al. 2011). Comparing the Atlantic Ocean part (Fig. 3B) and the Cantabrian Sea ( Fig. 3C), we can see that the concentration of chlorophyll is higher in the latter. The model does succeed in reproducing the spring peak, although the value is overestimated. In the MM5 forced run, the beginning of the bloom has a delay of around 2 wk in 2006 ( Fig. 3), whereas in 2007 the bloom timing is in reasonable accordance with the observations. The autumn bloom, which is a weak signal in the observations, is overestimated by the model. Moreover For 2006, if we used HIRLAM (red line in Fig. 3), the beginning of the spring bloom was better reproduced, and the spring chlorophyll peak, although still overestimated, was closer to the satellite ob servations. Indeed, although the correlation of both models with observations was similar in the common period (except for the Cantabrian Sea, which was lower for MM5), the RMSE was consistently lower when using HIRLAM (0.36 vs. 0.51 for the whole domain; 0.42 vs. 0.45 in the Atlantic Ocean; 0.59 vs. 0.92 in the Cantabrian Sea). To better understand the misfit between the model results and the observations in spring 2006, we plotted the model short-wave radiation against the observations of an AEMET radiometer (pyranometer sensor) located inland in the city of A Coruña (43. 367°N, 8.417°W). Fig. 4 shows the daily averaged global radiation from the radiometer (global radiation = direct radiation + diffuse radiation) and the net short-wave radiation from the MM5 and HIRLAM atmospheric models at the same location. Although albedo was not subtracted from the radio meter, it does not alter the conclusions. which had been strong in the area since mid-March, with a short period of relaxation at the beginning of April. In contrast, the Artabro Gulf (Fig. 5A, green area) experienced sustained upwelling conditions for most of the first half of April, with a small relaxation around 8 April (Fig. 6C), just when sampling was taking place in this area. Thus, the upwelling of deeper, colder waters explains the lower temperatures measured in this area (Fig. 5B). Sampling in the western Cantabrian Sea (Fig. 5A, orange area) comprised the period between 12 and 16 April 2006, which was characterized by weak winds (Fig. 6B). From 17 April on, the eastern Cantabrian Sea was surveyed (Fig. 5A, red area); this period was also characterized by rather calm wind conditions (Fig. 6A). Higher surface temperatures are observed as we advance to the east, which is usual during this period of the year in the southern Bay of Biscay. The most interesting feature of the surface salinity along Pelacus 2006 (Fig. 5F) is the freshwater plume around 43°N associated with the strong river discharge that This means that the model has some ability to replicate important features of the dynamics in the area (upwelling, river plumes or the IPC), although temperature was overestimated during upwelling conditions at the Artabro Gulf (around 1°C higher) as was the salinity of the freshwater plume around 43°N (more than 2 salinity units).
In contrast with Pelacus 2006, the Pelacus 2007 cruise started under strong upwelling conditions that persisted with variable intensity (Fig. 7A−E) from the beginning of the cruise in the western Iberia area (Fig. 5J, blue area) through the sampling in the Arta -bro Gulf (Fig. 5J, green area) and the western Cantabrian Sea (Fig. 5J, orange area). Hence, upwelling is responsible for the observed low surface temperatures in the 3 areas that are also reproduced by the model (Fig. 5K,L). Wind forcing in the eastern Can tabrian Sea (Fig. 5J, red area) was variable between 13 April and the end of the cruise; westerly winds were dominant the first 2 d and easterlies from 18 April on, with a period of relaxation in between. Surface temperatures increased to the east both in the observations and the model results. Surface salinity did not reflect the existence of a freshwater plume at the time of sampling in the western Iberian shelf (Fig. 5J, blue area), which is related to the fact that the river discharge was rather weak during this period (Fig. 7F). Nevertheless, some days later, freshwater was measured at the western and eastern Cantabrian Sea (Fig. 5J, orange and red areas), probably related to the river plume expansion under the dominant wind regime. The model was able to reproduce these low-salinity signals, although they were slightly overestimated (around 1 salinity unit; Fig. 5P).  2101  3101  2104  1101  1002  0203  1203  2203  1104  0101  2002  0104  0105   2101  3101  2104  1101  1002  0203  1203  2203  1104  0101  2002  0104  0105   2101  3101  2104  1101  1002  0203  1203  2203  1104  0101  2002  0104  0105   2101  3101  2104  1101  1002  0203  1203  2203  1104  0101  2002  0104  0105   2101  3101  2104  1101  1002  0203  1203  2203  1104  0101 2002 0104 0105  Fig. S4, right panel (Text S4). The model's performance is reasonably good; the model results were better at the surface than at 100 m depth.
We have proven that the model's ability to reproduce the hydrographic observations during the Pelacus cruises is fairly good. Short-term variability at the event scale (3−14 d) characterizes the system; therefore, conditions may vary strongly among hydrographical sections and can only be taken into ac count with a numerical model with realistic meteo rological forcing.  (Fig. 8B,F). To provide a hint on the state of the bloom, the percentage of diatoms obtained from the phytoplankton counts is included (Fig. 8C for 2006; Fig. 8D for 2007). Note that pico plankton is not included in phytoplankton counts although it can be a significant contribution to total phytoplankton biomass in some areas, as shown by Calvo-Diaz et al. (2004) in an analysis of Pelacus 2002 data. Finally, the concentration of surface NO 3 in mmol m −3 is also shown (Fig. 8G), although only for 2007 since this variable was not sampled in 2006.

Biological observations in
At the beginning of the Pelacus 2006 cruise, the concentration of chlorophyll measured on the western Iberian shelf (Fig. 5A, blue area) was low, with the exception of the section at 42° N (Fig. 8A), which was sampled on 4 April under strong downwelling conditions and high river discharge (see Section 3.2.1). Downwelling conditions cause the confinement of river plumes at the coast and their transport northwards (recall Fig. 5F,G). In this sense, retention, confinement and nutrient enrichment from rivers in the plume could explain the higher values of chlorophyll at the coastal stations. The high percentage of diatoms at 41°N, very close to the Douro mouth (Fig. 8C), can be an additional confirmation of this bloom at the plume. Zooplankton biomass was generally high at most of the sections in this area (Fig. 8F), which suggests that we are observing a late stage of a previous bloom. This hypothesis is supported by the phytoplankton counts corresponding to most of the stations in this area, which show a dominance of flagellates (low percentages of diatoms in Fig. 8C), except the section at 41°N as mentioned above. The concentration of surface chlorophyll was also low in the Artabro Gulf (Fig. 5A, green area) despite the occurrence of a weak upwelling pulse before sampling. The percentage of diatoms from plankton net hauls (depth-integrated in 100 m) was around 30%. The zooplankton biomass was lower than on the western Iberian shelf. Higher chlorophyll concentrations than in the previous areas were measured in the western Cantabrian Sea (Fig. 5A, orange area), with increasing values offshore. Dinoflagellates were the dominant phytoplanktonic species (> 90%), with the exception of the westernmost section (around 7°W), where diatoms dominated at the innermost station, probably influenced by the upwelling conditions of the previous days. Zooplankton biomass was also lower than on the western Iberian shelf. In the eastern Cantabrian Sea (Fig. 5A, red area), the concentration of chlorophyll decreased offshore but in creased inshore. Diatoms were the dominant phytoplanktonic species, and zooplankton concentration was low.
The situation was completely different during the Pelacus 2007 cruise, which started under upwelling conditions that had persisted since mid-March (see Section 3.2.1). This was reflected in the high concentration of chlorophyll, especially inshore (Fig. 8B), in the western Iberian shelf and Artabro Gulf (Fig. 5J, blue and green areas) and in the concentration of NO 3 in the same areas (Fig. 8G), which showed that the surface layers were being effectively fertilized. Apparently, we were observing an initial stage of the bloom because nutrients were not exhausted, chlorophyll was high and zooplankton biomass was still low. Moreover, phytoplankton counts in these areas showed that diatoms were the most abundant species; flagellates were almost negligible (Fig. 8D). The lowest concentrations of nutrients, lower concentrations of chlorophyll and higher zooplankton biomasses were found in the easternmost part of the Cantabrian Sea, especially close to the coast. Diatom relative abundance with respect to flagellates decreased to the east. These observations suggest that this area was sampled at a later stage of the bloom than the western part of the domain.

Variability in circulation and ecosystem
response in spring: model and satellite data The hydrographical and biological observations of the Pelacus cruises provide spatio-temporal information on the ecosystem; however, they are snapshots of a variable system and hence cannot completely characterize the spring season in the area. The upwelling index and knowledge of the dynamics of the area can help in evaluating conditions during the Pelacus sampling. Satellite images (when available) combined with biophysical models are the right tool to interpret the observations and especially to characterize the spring season and the response of the ecosystem. Fig. 9 shows the sequence of weekly average maps of MODIS surface chlorophyll and the corresponding model averages (considering the closest dates and  (Fig. 9D). Another interesting feature that is visible in the satellite images (Fig. 9B) and model results (Fig. 9E) during the first week of the cruise is the strip of low chl a along the coast and especially in the Artabro Gulf and western Cantabrian Sea. This low-chlorophyll band is consistent with the low observed in situ values at this time in the same areas (Fig. 8A). We also recall from Section 3.2.2 that the high chlorophyll concentration measured in one section of the western Iberian shelf was related to enhanced productivity at the river 60 Fig. 9. Weekly averages of surface chl a obtained from MODIS and the model for the period ranging between the week before the beginning of the Pelacus 2006 cruise and the week after its end. Units are mg m −3 plumes ( Fig. 8A) which, according to the satellite images ( Fig. 9B,C,G,H,I), was a sustained feature all along the cruise in this area. The model underestimates production at the river plume (Fig. 9E,F,J,K,L). From the second week of the cruise on the satellite and model images (Fig. 9C,F), there is a decline in the concentration of oceanic chlorophyll, especially off Portugal and in the eastern Cantabrian Sea. The following 3 weeks (Fig. 9G,H,I for the satellite and Fig. 9J,K,L for the model) seem to suggest a northwestward displacement of the chlorophyll maximum; the concentration in the remaining areas is low, mostly inshore, although information about most of the Cantabrian Sea is missing due to cloud coverage. This pattern is well reproduced by the model. The statistics of the model performance against satellite observations during the period of the Pelacus 2006 cruise are provided in Fig. S5 (Text S4). The average values are similar in both cases, although the model overestimates the chlorophyll in the eastern Canta brian Sea and underestimates it at the river plumes on the western Iberian shelf, which is clearly re flected by higher positive and negative percentage biases. Correlation is generally positive off Portugal and the Galician rias and in the central Cantabrian Sea but is poor in the western and eastern oceanic parts of the Cantabrian Sea.
In addition to surface chlorophyll, the model can provide information on the variability around the spring bloom, which is not possible to obtain from other data sources. To illustrate this ability, we include Hovmöller plots for 2 sections on the west coast (41.5°N and 42.5°N; Fig. 10) and 2 sections on the North coast (6°W and 3°W; Fig. 11), ranging from 1 January until 30 April for salinity, SST, mixed layer depth (MLD, in meters) calculated from the density profiles considering a threshold of 0.03 kg m −3 , chlorophyll (in mg m −3 ), NO 3 (in mmol m −3 ) and zooplankton (in mg C m −3 ) obtained from the model on a daily basis (12:00 h). Fig. 10D,J provides insight into the onset of the oceanic spring bloom along the western Iberian shelf in 2006, which started after 12 March, some weeks before the Pelacus 2006 cruise. In the Cantabrian Sea, the bloom began around 20 March (Fig. 11D,J). It is interesting to note that the thermal stratification was very weak in West Iberia at the beginning of the bloom (see MLD in Fig. 10C,I) and that a complete MLD shoaling did not occur in the Western Cantabrian Sea either (Fig. 11C). Only at the eastern Cantabrian Sea (Fig. 11I) was proper thermal stratification present at the bloom onset. The model can also be used to better understand the influence of river plumes on production. For instance, comparing the Hovmöller plots for 41.5 and 42.5°N in Fig. 10, some differences related to the degree of river in fluence arise. Since the river influence is higher at 42.5°N (Miño and the Galician rías) than at 41.5°N (north of the Douro river), the extension of the river plumes is wider in the first case (Fig. 10A,G). The low-salinity signal is correlated with low chlorophyll (Fig. 10J) and high zooplankton (Fig. 10L), a situation mostly sustained at the river plumes since the end of January (the same is observed at 41.5°N). Just after 22 March, an exceptional river discharge under neap tides (see the subplots corresponding to the Douro discharge and the tides at A Coruña in Fig. 6) led to very low salinity close to the coast and provided new nutrients (Fig. 10K). The 31 isohaline is included in the plots (magenta line) corresponding to this latitude to emphasize this fact. Confinement, stratification and enrichment led to the increase in the concentration of chlorophyll that is observed before and during the cruise (Fig. 10J) and the corresponding increase in zooplankton biomass (Fig. 10L). This is a plausible mechanism to explain the higher chlorophyll values at 42°N in Fig. 8A as well as the high zooplankton values (Fig. 8E). On the western Can tabrian shelf, the Hovmöller plots for 6°W show that chlorophyll production associated with the river plumes has existed since the end of January (Fig. 11D). Two upwelling pulses occurring at the end of February and around 12 March (Fig. 6B) caused the northward expansion of the plumes and the shoaling of the MLD in the river plume area (Fig. 11C), which promoted production (Fig. 11D). It is interesting to note that the nutrients are consumed in the plume in the course of the northward expansion (Fig. 11E) during the 12 March pulse and that zooplankton grow during the 2 pulses (Fig. 11F). Upwelling-favourable winds are dominant during most of April 2006 in this area (Fig. 6BK), explaining the river plume expansion that is observed during the cruise. It seems that the river input of nutrients during this period is enough to sustain some production but at lower levels than offshore, where the conditions seem to be more favourable (Fig. 11D,E). This is just the opposite of the situation observed during the winter months.
The weekly-averaged satellite images and model results from the week before the Pelacus 2007 cruise until the week after its end are shown in Fig. 12. High concentrations of chlorophyll were observed the week before the cruise in the satellite images (Fig.  12A) and model results (Fig. 12D) in the oceanic and coastal areas of Portugal and Galicia, but not in the Cantabrian Sea. The satellite and model maps for the first week of the cruise are mostly cloud covered in this area (Fig. 12B,E); however, in the second week, it is clear that the bloom has extended to the whole Cantabrian Sea (Fig. 12C,F). During the first and second weeks of the cruise, sampling was taking place mostly on the western Iberian shelf and Artabro Gulf under upwelling-favourable winds, which is reflected in the high chlorophyll concentrations on the coast in these areas observed both in the satellite (Fig. 12B,C) and model images (Fig. 12E,F) as well as in the in situ observations (Fig. 8B). The second week of the cruise also covered the sampling in the western Cantabrian Sea, where the concentrations of chlorophyll were high, not only on the coast but also in the open ocean (see Fig. 12C for the satellite and Fig. 12F for the model results, which over estimate the chlorophyll concentrations). This is in agreement with the in situ observations as shown in Fig. 8B. The remaining weeks (Fig. 12G,H,I for the satellite and Fig. 12J,K,L for the model) seem to show the decline of the oceanic bloom, which is reflected by a northwestwards displacement of the maximum chlorophyll, leaving low concentrations along the coast with the exception of the Artabro Gulf and the river plumes on the western Iberia shelf, which are seen by the satellite but not well reproduced by the model. Fig. S6 (Text S4) summarizes the statistics of the comparison of the model with the satellite product in the Pelacus 2007 period. The average values are similar off Portugal and the Galician rias (low biases), but they differ especially in the Cantabrian Sea, with a strong overestimation in the easternmost part and underestimation in the westernmost part, which are reflected as positive and negative biases, respectively. The correlation is generally positive, with the highest values occurring in the eastern Cantabrian Sea. Correlation is negative mostly in the western Cantabrian Sea and the coastal areas of the Artabro Gulf.
As with 2006, we used Hovmöller plots obtained from the model results to complete the picture of variability around the 2007 spring transition in the area. Fig. 10P,V shows that high concentrations of chlorophyll are calculated in the open ocean some days after 12 March in close relation with the stratification imposed by the offshore expansion of the river plumes (Fig. 10O,U) and the increase in solar radiation around this date (see Fig. 4C for March 2007). A second and more intense bloom occurred after 11 April, when stratification of the water column was stronger (Fig. 10O,U). In the Cantabrian Sea, the oceanic bloom started around 11 April at 6°W (Fig. 11P), whereas in the easternmost part, at 3°W (Fig. 11V), it had started at least 10 d earlier according to the model results. However, we recall that this is the area that presents higher biases with respect to observations (Fig. S6).
Focusing on the upwelling taking place on the western Iberian shelf, we see that some days before the cruise, from approximately 20 March onward, the dominant winds (see Fig. 7) caused the expansion of the river plumes far offshore and even their separation from the coast (Fig. 10M,S). Moreover, a decrease in river discharge some days before the cruise (Fig. 7F) led to a contraction of freshwater onshore (Fig. 10M,S) and a decrease in stratification, especially observed at 42.5°N (Fig. 10U). The uplift of nutrients at the coast as a consequence of upwelling is clearly seen in Fig. 10Q,W, leading to the increase in the concentration of chlorophyll that is observed in Fig. 10P,V. Notice that the peak in chlorophyll occurs some days after the cruise sampling in this area, confirming that we are at an early stage of the bloom, as we suspected from the analysis of the observations (see Section 3.2.2). This fact is also supported by the low zooplankton concentrations at this time, which peak just after the cruise (Fig. 10R,X).
Production associated with river plume variability can also be seen in the Cantabrian Sea Hovmöller plots. In Fig. 11P, we see that the highest chlorophyll concentrations in the western Cantabrian Sea were registered between 12 March and the beginning of April, associated with the northward expansion of river plumes and the shoaling of the MLD in the river plume (Fig. 11O). Prevailing upwelling conditions during April (see Fig. 7B) introduced new nutrients and thus sustained production. The zooplankton biomass started to increase from the beginning of April onwards (Fig. 11R). Notice that the concentration is higher during sampling (around 11 April) than that measured on the western Iberian shelf (at the beginning of April; Fig. 11R,X), which is in agreement with the observations at the coastal stations (recall Fig. 8F). In the eastern Cantabrian Sea, Fig. 11V,W,X confirms that we are at a later stage of the chlorophyll bloom during sampling, even in the open ocean, characterized by low concentrations of chlorophyll, low concentration of nutrients and high zooplankton biomass, which helps to contextualize the observations in this area (Fig. 8).

Stratification and timing of the spring bloom
Phytoplankton growth has traditionally been linked to the development of the seasonal thermocline at the beginning of spring. This is formally known as the critical depth hypothesis (Sverdrup 1953), which predicts that blooms occur when the MLD shoals to become less than a critical value at which gross primary production and respiration are balanced. In addition to this hypothesis, other mechanisms have been proposed to explain the onset of the spring bloom, the most relevant of which are the critical turbulence hypothesis (Huisman et al. 1999), which states that a phytoplankton bloom can occur independent of the MLD as long as the local turbulence is low enough to keep the cells in the illuminated surface layer, and the disturbance−recovery hypothesis (Behrenfeld 2010, Behrenfeld & Boss 2014, according to which the phytoplankton blooms start in winter when the ocean surface is cooling and mixing is strong. These conditions would be favourable for the accumulation of phytoplankton due to reduced encounter rates with grazers as a consequence of dilution. The resulting bloom increases prey−predator interactions, and predation ultimately consumes the bloom in the recovery phase. The mechanisms that trigger the onset of the spring bloom remain an object of research and scientific debate (see, for instance, the review by Chiswell et al. 2015 andthe comments from Marra 2016). This discussion has been attributed, in part, to a lack of data on phytoplankton stocks in winter and early spring, which the relatively recent availability of new observational platforms with high-resolution sampling capabilities (such as gliders) can help to solve (Rumyantseva et al. 2019).
We have seen that our model presented a reasonable ability to reproduce the timing of the spring bloom; hence, we can use it, together with meteorological information, to provide insight into the variability of mechanisms that trigger the onset of the spring bloom at different spatial and temporal scales. Our model results have shown that the onset of the bloom occurred under different hydrodynamic and atmospheric conditions within the same year in different geographical areas of our computational domain and that a combination of different processes can trigger the bloom at the same location. For instance, the onset of the bloom in NW Iberia in mid-March 2006 was not linked to stratification of the water column (Fig. 10C,D). Rather, it seems to be consistent with the critical turbulence hypothesis since it occurred in a period of wind relaxation (Fig. 6D,E), and wind is the main generating force of turbulence in the water column. Blooms under similar conditions have been reported in the North Atlantic (Townsend et al. 1992, Ellertsen 1993, the North Sea (Backhaus et al. 1999, Dale et al. 1999) and, more locally, in the north and northwest shelf waters of the Iberian Peninsula (Fernández & Bode 1991) and in the Central Cantabrian Sea (Álvarez et al. 2009). However, as pointed out by Álvarez et al. (2009), it is not as simple as having a physical mechanism that keeps the cells in the euphotic layer of the water column. A physiological activation phase, driven by increased solar radiation, is necessary to explain these blooms. Indeed, we see that solar radiation remained at high levels between 11 and 17 March 2006 (Fig. 4A), which could have provided these initial favourable conditions. On the other hand, the beginning of the bloom the same year in the eastern Cantabrian Sea is most likely explained by the classical critical depth hypothesis (Fig. 11I,J).
The processes that trigger the onset of the bloom and their variability can be further illustrated with our model results for the western Cantabrian Sea. In Fig. 13, we see that in 2006, phytoplankton started to grow and accumulate around 21 March, a period that is characterized by low turbulence (Akv is a ROMS model output variable that is associated with the vertical mixing in the momentum equations: low values indicate low vertical turbulence). In the previous days, short-wave radiation was significantly high, meaning that the light conditions might have been favourable for phytoplankton growth.  Fig. 7). The ability of the model to properly reproduce the timing of the spring bloom implies that it can be used as a tool to explore further aspects of the ecosystem. For instance, the timing and intensity of the spring bloom has been associated with the survival of haddock larvae, and hence recruitment, off the eastern continental shelf of Nova Scotia, Canada (Platt et al. 2003, Platt & Sathyendranath 2008, or with the sequestration of atmospheric CO 2 below the permanent thermocline in the sub-polar North Atlantic due to the sinking of particulate organic carbon (Martin et al. 2011).

River plumes and fronts
The dynamics and variability of river plumes in W and N Iberia have been the object of previous detailed studies ( Ruiz-Villarreal et al. 2006, Mendes et al. 2017. The response of river plumes to wind events is well described in our model, which reproduces the expansion (contraction) of the river plumes under upwelling-(downwelling)-favourable winds (Otero et al. , 2009). See, for instance, the expansion of the river plume just before 20 February 2006 (Fig. 10A) under upwelling conditions (Fig. 6). Additionally, the impact of cooling−warming on the river plume is well simulated. In our results, we clearly reproduce its thermal signal, with lower values than the adjacent waters (Ruiz-Villarreal et al. 2006, Otero et al. 2009). In addition to the re sponse of the plume to meteorological events, our model realistically accounts for the freshwater content on the shelf and its variability since it is mainly forced by daily river runoff. Prevailing downwelling conditions in autumn−winter in the area, although interrupted by some upwelling pulses, result in the ac cumulation of freshwater on the shelf , which is exported offshore in spring with the change to upwelling conditions (Ruiz-Villarreal et al. 2006). Our model results also provide insight into the effect of river plumes on plankton productivity on the W and N Iberian coasts (Figs. 10 & 11), in line with previous studies in the area (Santos et al. 2001, Álvarez-Salgado et al. 2002, Bode et al. 2002, Ribeiro et al. 2005) that have pointed out that river plumes provide favourable environments for enhanced primary production due to haline stratification (which promotes stability of the water column and concentration at the surface layers, where light is available) and nutrient enrichment. Interestingly, our model results show productivity associated with the river plumes as early as the end of January in 2006 (Fig. 10D,J) and even earlier in 2007 (Fig. 10P,V) in W and N Iberia. Winter blooms have been described for the Portuguese shelf (Santos et al. 2001, Ribeiro et al. 2005, in the Galician rias (Varela et al. 2008(Varela et al. , 2010 and in the central Cantabrian shelf (Álvarez et al. 2009) associated with unusually calm, sunny weather characteristic of high-pressure periods as well as in some river plumes in the Bay of Biscay (Loire Estuary; Gohin et al. 2003) and Gironde estuary (Labry et al. 2001). In the Portuguese and French shelf river plume examples, the combination of haline stratification and optimal light conditions was found to trigger the bloom. Calm conditions with sufficient light are more frequent on the west coast (see Figs. 6 & 7, where it is clear that downwelling is more intense on the north coast), and the accumulated freshwater on the western shelf is also much higher than on the northern shelf (Figs. 10 & 11). The influence of river plumes in the Cantabrian Sea has been described as reduced and limited to coastal zones due to the relatively low discharges associated with the reduced drainage basins (González-Nuevo & Nogueira 2014). Our results show that river plumes have a clear impact on the chlorophyll signal in the narrow Cantabrian Sea shelf (see Fig. 11D,J,P,V) and that the stronger runoff of French rivers clearly impacts the easternmost part. Therefore, river plumes in the Cantabrian Sea need to be considered for understanding primary and secondary production in winter and spring, and this has not been sufficiently recognised so far in the literature. Al though the summer season is outside the scope of this paper, we should mention that freshwater advection from French rivers during summer upwelling results in freshwater intrusions into the Cantabrian Sea shelf .
River plume frontal zones have also been related to high productivity due to convergence and strong mixing (i.e. Acha et al. 2004), thus constituting food and reproductive habitats for higher trophic levels such as zooplankton and fish (Morgan et al. 2005). The model presented here is a suitable tool to investigate further productivity at river plume frontal zones since it realistically reproduces factors affecting fronts like river plume variability and the variability of ambient circulation (Otero et al. 2009).
The interannual variability of the river plumes between 2006 and 2007 is remarkable in the model results, showing a wider extension in 2007 than in 2006 due to higher river runoff (see Fig. 10F vs. Fig. 10R and Fig. 6 vs. Fig. 7). A particularly strong expansion event took place in W Iberia between the middle and end of March 2007 (Fig. 10M,S) caused by high river discharge and sustained upwelling. The high primary and secondary productivity associated with the river plume extend significantly offshore and illustrate the potential of river plumes to transport material offshore. The interannual variability of river plumes, together with the variability in the dynamics of other relevant processes, like upwelling, result in strong variability in the oceanographic and biogeochemical conditions on the shelf, with implications that extend to the whole food web. In particular, such variability could be responsible for part of the variation in recruitment of several fish species for which the environmental con ditions at the river plumes have been found to be favourable for larvae survival (i.e. sardine and an chovy; Santos et al. 2004, Bernal et al. 2007, Rod riguez et al. 2009

Model simplicity and parameterization
Rather than considering the simplicity of an ecological model such as ours a disadvantage, it has been suggested that simple models can be as good as complex ones when adapted for the particularities of a specific site (Friedrichs et al. 2007). Indeed, complex models do not necessarily guarantee better results due to their uncertainties in the description of ecological processes, lack of data and difficulties involved with parameterizations. These concerns about the efficacy of increasing model complexity were analysed some years ago (i.e. Anderson 2005) and, more recently, questions have been addressed such as how complex a model should be to reliably represent the interplay between biology, physics and biogeochemistry (Kriest et al. 2010, Ward et al. 2013, Kwiatkowski et al. 2014, Xiao & Friedrichs 2014. In this context, we have selected one of the simplest ecological models within ROMS (Fennel et al. 2006), putting our effort into evaluating the model's skill against satellite and in situ observations and into the selection of parameters for proper characterization of the ecosystem -especially in the spring season we wanted to simulate. We recall that in our model, we only account for one group of phytoplankton and one group of zooplankton to which we have attributed certain constant parameters.
The model's ability to reproduce the bloom dynamics in the area was satisfactory in general, although the spatial statistics showed the challenges of pointto-point comparisons with satellite products, as discussed in Stow et al. (2009) andIOCCG (2000), and chlorophyll was overestimated during the spring and autumn blooms (Fig. 3). Indeed, in the first stages of a bloom, when nutrients have not yet been consumed, diatoms are the dominant species, being re placed by dinoflagellates when nutrients are de pleted (Tilstone et al. 2003). The parameters selected for the only group of phytoplankton in our model correspond to a rapid-growing species/functional group that would mimic the ones present at the beginning of the spring bloom. Therefore, we impose this rapid-growth behaviour onto our unique phytoplankton group with independence in the amount of available nutrients, thus not accounting for observed species succession as the bloom progresses. This fact may have contributed to the observed chlorophyll overestimation, although we imposed a control on phytoplankton growth by increasing g max with respect to the default reported by Fennel et al. (2006), which also contributes to avoiding fast depletion of NO 3 . In García-García et al. (2016), we showed that this value of grazing resulted in some qualitative skill in comparison to in situ measurements of zooplankton bio mass during the Pelacus cruises in spring 2006 and 2007. Tweaking model parameters to adjust default values (generally defined for a different area) to local conditions is a common procedure in biogeochemical modelling. Table 3 shows some examples of how different ROMS-based biogeochemical models that include part of our study area had their parameters adjusted (focusing only on the ones we modified) with respect to the reference model. In some cases, the changes were done in the same direction as we did, although not simultaneously; i.e. K NO3 was decreased with respect to the default value (Marta-Almeida et al. 2012, Reboreda et al. 2014b) and g max was increased (Dabrowski et al. 2014), but the model parameters are all different for the different models. In all cases, the selected parameters led to an improvement of the model's performance with respect to the default, although most of the models, including ours, overestimated chlorophyll.
Biogeochemical model parameterizations are indeed relevant for reliable modelling of plankton 70 Reference Biogeochemical model type K NO3 α g max Model performance and Reference (chl a) This paper N2PZD2 (Fennel et al. 2006) 0.25 (0.5) 0.05 (0.125) 1 (0.6) Overestimation spring bloom. Good timing Dabrowski et al. (2014) N2PZD2 (Fennel et al. 2006) 0.5 (0.5) Not specified 1.3 (0.6) Overestimation spring bloom. Good timing Marta-Almeida et al. (2012) NPZD (Koné et al. 2005) 1.5 (2) 1 (0.004)* 0.9 (0.9) Overestimation. Variable timing Reboreda et al. (2014b) NPZD (Koné et al. 2005) 1.5 (2) 1 (0.004)* 0.9 (0.9) Overestimation spring bloom. Anticipation Reboreda et al. (2014a) N2ChlPZD2 (Gruber et al. 2006) 0.9 (0.75) 1 (1)* 0.6 (0.6) Overestimation spring bloom. Anticipation Table 3. Parameter values of different studies that cover part of the study area with respect to the reference biogeochemical model. The focus is on the parameters modified in this study (see Table 2 for the meaning of symbols and units). Default values are inside brackets. The values with an asterisk have different units than the ones above (see Koné et al. 2005) dynamics, but the role of the physical model is crucial and, in many cases, better explains model misrepresentations (Doney et al. 2004, Skogen & Moll 2005, Najjar et al. 2007, Sinha et al. 2010, Popova et al. 2012). The biogeochemical model results are also dependent on meteorological forcing, either indirectly through the impact of ocean circulation  or directly through radiation, as we have illustrated in this study. With respect to the indirect effect through changes in ocean circulation,  remarked that, in general, the time scales of wind events are similar in all meteorological model products (although some events are not captured by the models), which results in a similar qualitative response of circulation to wind because the general conditions of upwelling− downwelling are represented by all the models. This means that we should expect both model configurations (model forced with MM5/WRF or HIRLAM) to reproduce the event variability in a similar way. It is the differences in radiation into the ocean between models which seems to produce the first-order effects, and we have shown that if the radiation is underestimated around spring, then a delay in the spring bloom onset is likely to occur in the model (see Fig. 3 for MM5/WRF in 2006 and Fig. 4). Finally, another important feature in biogeochemical modelling is the choice of the numerical scheme. For example, Lévy et al. (2001) found that the new production estimates could vary as much as 30% by comparing different advection schemes in a simplified channel simulation, and Kriest & Oschlies (2011) also highlighted the importance of the numerical schemes considered in the calculation of organic matter sedimentation and remineralization in biogeochemical ocean models. We attribute part of the model skill against observations to the revised handling of terms in the ecological model, which reduced negative values in biogeochemical fields and allowed us to avoid the use of a positive definite algorithm, such as MPDATA for tracer (including salinity and temperature) advection. Other implementations of NPZD-like models in ROMS use the MPDATA for horizontal (and vertical) advection of all tracers, including temperature and salinity (Dabrowski et al. 2014, Solé et al. 2016, Wang et al. 2020. However, MPDATA is a non-monotonicity preserving, total-variation-diminishing scheme, meaning that over-shoots and under-shoots may appear in the numerical solution. In our coupled hydrodynamic− biogeochemical model configuration, using MPDATA resulted in a worse representation of the temperature and salinity fields than when we used the U3H scheme for the hydrodynamic-only simulations, although it guaranteed the positive definition of the biological tracers. Examples of NPZD-like models in ROMS that do not use MPDATA can be found in the literature. For example, Fennel et al. (2011) used a third-order upwind advection of momentum and a fourth-order horizontal advection of tracers in a study of the biochemistry of the shelf affected by the Mississippi River plume. A possible improvement of our model performance could be the splitting of the phytoplankton compartment into 2 functional groups, representing diatoms and dinoflagellates, respectively. Koné et al. (2005) showed that a model with 2 phytoplankton groups (although size-based) provided more realistic chlorophyll distribution in the southern Benguela upwelling system. Another possibility for better representation of plankton dynamics is the treatment of the model parameters as time−space-varying rather than constant. This approach was used by Mattern et al. (2012), who used statistical techniques that minimized the difference between satellite and modelled chlorophyll to estimate time-dependent values of 2 parameters in a 3D biogeochemical model for the coastal region of the northwest Atlantic, and they obtained an improvement in the chlorophyll comparison to observations.

Nutrient input by rivers
We have observed that the modelled concentration of chlorophyll at the Galician and Portuguese river plumes was underestimated with respect to the satellite values in 2006 and 2007, especially after the bloom decline (Figs. 9 & 12). The estimation of chlorophyll from ocean colour in coastal waters is challenging because the optical response is affected not only by phytoplankton but also by other particles and dissolved material. In this sense, although we explicitly removed the pixels flagged as chlorophyll-warning, failure or turbid, the OC3M algorithm used to process the satellite product in our study is likely to have resulted in an overestimation of the satellite chlorophyll on the shelf in Figs. 9 & 12. Note that chlorophyll overestimation might affect waters further away from the coast on the eastern Bay of Biscay since the French shelf is broader than the N and W Iberian shelf. The way chlorophyll is calculated in the model from the phytoplankton biomass considering the Geider et al. (1996Geider et al. ( , 1997) models (further details in Fennel et al. 2006) is not without uncertainty since the models employed are simplifications of reality.
On the other hand, these differences can also be related to the effect of nutrient inputs by the rivers in the model. In Section 2.1.2, we pointed out that the selected constant concentration of NO 3 input by the rivers was in the lower range of the values that we found in a review of several literature and data sources. However, the values for the NO 3 input from rivers, although scarce and scattered, are variable. For instance, according to Azevedo et al. (2006), NO 3 input by river Douro ranges between 57.1 ± 8.9 mmol m −3 close to the river mouth and 103.9 ± 8.9 mmol m −3 at an upper location, with both samples taken during flood tide over an annual cycle (December 2002− December 2003. In Couto et al. (2018), averaged concentrations of NO 3 at the Douro lower basin between 76 and 171 mg l −1 (1226−2752 mmol m −3 ) were measured at low tide in spring (May 2013). These high values were attributed to intense precipitation in the sampling period associated with intense soil lixiviation. Also, the biogeochemical characterization of river Miño in the zone of tidal influence that was carried out as one of the tasks of the national project ZOTRACOS (2003ZOTRACOS ( −2006 obtained NO 3 concentrations of around 50 mmol m −3 at low tide and 20 mmol m −3 at high tide, which is twice as high as the concentration imposed for river NO 3 in our model configuration. Similar averaged NO 3 concentrations were measured at some stations of the Confede ración Hidrográfica Miño-Sil network (https://www. chminosil.es/) in the lower basin of the Miño river.
This discrepancy among the different literature and data sources raises the question of whether the underlying reason for the chl a underestimation in our model at the Galician and Portuguese river plumes after the bloom decline could be related to an underestimation of NO 3 input by rivers Miño and Douro, the largest rivers in that area. To evaluate this hypothesis, the period ranging from 1 January 2006 until 31 May 2006 was re-run, considering a concentration of 100 mmol m −3 of NO 3 for river Douro and 50 mmol m −3 for river Miño. Similar values were used in Marta-Almeida et al. (2012) (88 mmol m −3 of NO 3 for river Douro and 34 mmol m −3 for river Miño). Note that the nitrate units in their Table 2 should be mmol m −3 instead of μg l −1 . Fig. 14A,B,C shows the MODIS daily satellite images for some time after the bloom decline when the river plumes at the Galician rías are not cloud-covered. In the same figure, the results of the reference model configuration (10 mmol m −3 of NO 3 input for every river) regarding surface chl a for the time step closer to the date of the satellite image are included (Fig. 14D,E,F). Finally, the last row of Fig. 14 (Fig. 14G,H,I) corresponds to the model results regarding chl a when the NO 3 concentration of rivers Miño and Douro has been increased.
It can be seen that increasing the concentration of NO 3 in rivers Douro and Miño has an impact on the concentration of chlorophyll at the river plumes, which also increases, becoming closer to the observations (although still underestimated). Therefore, imposing realistic time-varying fluxes of NO 3 through the rivers, similar to the realistic river flows we are including in our model, will result in an improvement in the model results. Indeed, an enhancement of the land discharge forcing into models has been claimed as one of the priorities for improving (operational) models (see, for instance, Capet et al. 2020 and references therein). Therefore, the effort to improve land discharge forcing should include constructing time series of nutrient fluxes from rivers.
The above comments on NO 3 inputs are especially relevant in the case of models that assume primary productivity is nitrogen-limited, as in our case. This assumption, which is common in marine systems (Eppley et al. 1973), implies that if we increase the concentration of NO 3 it would immediately be consumed under favourable light conditions, leading to an increase in primary productivity. However, in other coastal areas under the influence of river plumes, phosphate has been shown in modelling exercises to limit primary productivity; for example, in the northern Gulf of Mexico affected by the Mississippi and Atchafa laya plume (Laurent et al. 2012) or in upwellingdominated systems such as the Oregon shelf (Spitz 2010). Indeed, the role of phosphate limitation on phytoplankton growth is increasingly being recognised, especially in coastal waters under the influence of freshwater discharges (e.g. Labry et al. 2002). For instance, Laurent et al. (2012) proved, by means of numerical modelling, that phosphate limitation would slow down NO 3 uptake in the Mississippi River plume, resulting in a wider spreading of this nutrient -and hence, primary productivity -over a larger area.
According to Azevedo et al. (2006), the Douro estuary is generally phosphate-limited, except when river discharge is low. Therefore, it would not be surprising if the river plumes on the western Iberian shelf were phosphate-limited as well. In this sense, an increase in the realism of our configuration would be achieved not only by properly characterizing variability in NO 3 fluxes, but also by increasing the complexity of our biogeochemical model by accounting for phosphate. The effect of including phosphate dynamics is likely to have an impact in our area, and primary productivity would probably spread over a larger area as Laurent et al. (2012) found with the Mississippi plume, thus better fitting the satellite chl a images and improving the representation of the ecosystem dynamics.

SUMMARY AND CONCLUSIONS
A 3D biophysical model for western and northern Iberia has been implemented. The physical model performs well in reproducing the main oceanographic conditions in the area (upwelling, river plumes, fronts, stratification), and the biogeochemical model accounts for the sources of nutrients (nutrient content of ENACW that uplifts at the coast under upwelling-favourable winds and river input). Comparing the model results with in situ and satellite observations, we see that a simple N 2 PZD 2 model, adapted to the particularities of the area, provides valuable information on phytoplankton dynamics at different temporal and spatial scales. At a seasonal time scale, the onset of the phytoplankton blooms, together with their duration and interannual variability, are reasonably simulated compared with a satellite ocean colour product, although the chlorophyll concentration is slightly overestimated. At shorter time scales, we compared the model with in situ observations (including plankton) and put into relevance the capacity of the model to be used as a tool to provide information on the data gaps of the non-synoptic sampling spring cruises that the IEO carries out every year, providing at the same time the geographical differences among areas as well as a contextualization of the observations. In this sense, we have demonstrated how to take advantage of monitoring data -specifically, plankton data -to validate and complement a biogeochemical model. Compared with observations, our realistic numerical model configuration shows that productivity on the shelf is the result of the interplay of wind events, river plumes and light intensity, all of which vary at interannual, seasonal and event scales. This variability is propagated through the food web via trophic interactions, thus affecting the whole ecosystem. Other specific findings of our study are summarised in the following: (1) We have shown that our biophysical model is able to reproduce seasonal variability of plankton biomass, showing the typical spring and autumn blooms in the area. It also captures event variability and therefore constitutes a support tool not only to interpret the observations but also to investigate the causes of the spring bloom. We have demonstrated that the spring bloom is linked to the development of the seasonal thermocline, as expected, but radiation in the meteorological forcing used for the biophysical model is essential to obtain accurate timing of the spring bloom. We have shown that if atmospheric model radiation is underestimated around spring, then a delay in the spring bloom onset is likely to occur in the biochemical model.
(2) Our model results indicate that part of the interannual variability of production is driven by variability of the NO 3 content of intermediate waters (ENACW). In our model configuration, we can account for this interannual variability, driven by convective mixing, using temperature outputs from the parent model at the open boundary and translating temperature values into NO 3 through empirical temperature−NO 3 relationships obtained from in situ data of ENACW.
(3) Although upwelling is a main source of nutrients on the NW Iberian shelf, proper characterization of the nutrient input and its temporal variability by rivers is necessary to better explain the observed levels of productivity in the coastal areas, especially when upwelling conditions relax. River plumes promote and sustain productivity in the coastal areas of W and N Iberia by generating the appropriate conditions for phytoplankton growth: water column stability and nutrient enrichment. Our model results highlight that river plumes in the Cantabrian Sea need to be considered for understanding primary and secondary production in winter and spring -a fact that has not been clearly recognised so far in the literature. The model results, in agreement with observations described in the literature, show that blooms occur at the river plumes even in winter, associated with calm, sunny weather conditions that favour phytoplankton growth. These calm conditions with ade-quate light are more frequent on the west coast, where the accumulated freshwater on the shelf is also much higher than on the northern coast. This results in higher productivity in autumn−winter in the western area. We have learnt with our model simulations that considering the time variability of nutrient input through rivers would be necessary to reproduce the observed levels of chlorophyll on the shelf.
(4) From a modelling point of view, the biogeochemical model parameterizations are important in order to adjust the timing and intensity of the blooms in simple models such as the one used here, but a physical model that includes realistic forcing (wind, radiation, rivers, tides, etc.) is key to reproducing the oceanographic conditions that trigger the plankton blooms and determine their distribution. Using 2 model configurations with different meteorological forcing, we have shown the influence of physical forcing (in particular, short-wave radiation) on the timing and intensity of the spring bloom. The model seems to represent the main features of the bloom dynamics with no need to increase the model's complexity. However, part of the differences in chlorophyll between the model and observations we report (overestimation in spring or different spatial distribution) could stem from the fact that we have chosen a broadly used but still too simplified ecosystem model approach. It is likely that a slight increase in the complexity of the biogeochemical model, in terms of increasing the number of phytoplankton groups or considering phosphate limitation, will result in an im proved representation of species succession along the bloom and the plankton distribution on the shelf affected by river plumes. The selected numerical schemes also play a role in model performance: we attribute part of our model's skill against observations to the revised numerical handling of terms in the ecological model, which reduced negative values in biogeochemical fields and allowed us to avoid the use of a positive definite algorithm such as MPDATA for tracer (including salinity and temperature) advection.
(5) Through this study, we have shown the strengths and weaknesses of the model and have provided details on where it needs improvement. In other words, we have demonstrated that the model is not perfect. However, the observations we have around the spring bloom are not perfect or complete either. By combining the available observations and the model, through its construction and results, we have provided an example of how both data sets together can help to 'disclose the truth' (Skogen et al. 2021) about the relevant processes around the spring transition in NW Iberia.