1 Introduction

The mixed layer, where all tracers are almost homogeneous, owes its own homogeneity to mixing processes caused by turbulence injected from the atmosphere. It defines the primary extent of the exchange of energy and mass between the atmosphere and the underlying ocean and the penetration of turbulence. Changes within this layer regulate the ocean–atmosphere as a coupled system through the heat storage preserved on it and its influence on surface currents (McCreary et al. 2001; Seager et al. 2002; Montegut et al. 2004; Willis et al. 2004). Besides that, changes within this layer are a determinant of variability in physical and chemical properties that control the biological productivity in the ocean. Most mixed layer depth (MLD) variability studies have traditionally been focused in short-term time scales (i.e. diurnal, intra-seasonal and seasonal), and these frequencies dominate MLD variability. However, the lack of long time series of high-resolution vertical profiles in the ocean has made long-term studies of MLD variability to be quite infrequent.

The combination of long time series of hydrographic data sets and oceanic models is a quite common feature in mixed layer studies. The complex interactions among physics, chemical and biological processes have motivated the development and use of a number of models which try to give insight into these relationships. In this work, we will combine the free downloadable one-dimensional water column model General Ocean Turbulence Model (GOTM; available at http://www.gotm.net) with a long-term oceanic time series running since the early 1990s in the southern Bay of Biscay as part of an interdisciplinary monitoring program.

A number of simulations will be carried out with the aim of obtaining a better understanding of the local patterns of variability in the upper ocean in relationship with atmospheric forcing. At the same time, these exercises allow to evaluate the performance of GOTM under different model configurations and diverse forcing data and products against the actual series of profiles. After getting confidence in the results yielded by the model, we will attempt to reconstruct the historical evolution of the upper ocean, including the winter MLD variability, since the 1950s.

The structure of this paper is as follows: In Sections 2, 3 and 4, the GOTM model, data sets and study cases to be carried out are described. The results and discussion are presented in Sections 5 and 6, beginning with the simulation of the upper level variability in oceanic waters of the Santander standard section (Sections 5.1 and 5.2) and then focusing on the analysis of the ability of the climatological profiles and the GOTM model to reproduce the upper layer variability (Section 5.3). After that, preliminary results obtained from the application of the GOTM model to reproduce MLD variability in the last 60 years of the Bay of Biscay are analysed (Section 6). Finally, in Section 7, the main findings of the study are summarized.

2 The GOTM model

The study of the variability in mixed layer depth and upper ocean properties in oceanic waters of the southern Bay of Biscay will be addressed by means of simulations using the one-dimensional water column model GOTM (Burchard 1999, available at www.gotm.net). GOTM simulates vertical mixing and advection processes in the ocean forced by momentum and buoyancy fluxes at the surface. A brief and extended review of GOTM model equations can be found in Burchard et al. (2006) and Umlauf et al. (2008). The equations for temperature (θ) and salinity (S) vertical distribution are written as follows:

$${\partial_t \theta} - \nu' {\partial_{zz} \theta} + {\partial_z \langle w' \theta' \rangle} = \frac{\partial_z I}{\rho_0 {\kern 1pt} C_p {\kern 1pt}} + {\tau_r}^{-1} (\theta_{\rm osb} - \theta) $$
(1)
$${\partial_t S} - \nu'' {\partial_{zz} S} + {\partial_z \langle w' S \rangle} = S (P-E) + {\tau_r}^{-1} (S_{\rm osb} - S) $$
(2)

where ν′ denotes the molecular diffusivity of heat, C p denotes the heat capacity and ν′′ denotes the molecular diffusivity of salt. The only source term in these temperature (Eq. 1) and salinity (Eq. 2) balance equations results from the vertical divergence of shortwave radiation, I, and net precipitation (P)–evaporation (E) balance, respectively. The second term on the right-hand side describes a relaxation of θ and S with time scale τ r towards prescribed profiles θ obs and S obs. This relaxation term is used to account for lateral unresolved effects. Vertical velocities are set to zero in this application; consequently, no vertical advection is considered.

GOTM allows a choice from the different turbulence schemes, especially second-order turbulence closures, for simulating vertical turbulent diffusion in geophysical flows. Recent works from Burchard and Bolding (2001) and Umlauf and Burchard (2005) described the soundness and reliability of these second moment turbulent closure schemes in representing the development of mixed layer through entrainment in stratified shear flows. We will use the k − ϵ model with second-moment closure by Cheng et al. (2002) which has been shown to realistically simulate vertical mixing (Umlauf and Burchard 2005). A value of steady-state Richardson number of Ri st  = 0.25 and length scale limitation were used. Parameterizations of mixing by internal waves and double diffusion in stable stratification following Large et al. (1994) were introduced.

3 Data set

3.1 The Bay of Biscay

The present work deals with upper water variability in the southern Bay of Biscay, so we begin with a short description of the area. The Bay of Biscay is an adjacent sea to the Eastern North Atlantic basin, and the circulation patterns are weak in comparison with the main North Atlantic currents west of the Iberian Peninsula, the North Atlantic Current and the Azores Current (for reviews on the area see Lavín et al. 2006). The main circulation is the residual of the dominant southwards currents at the west, and central waters (Eastern North Atlantic Central Water, ENACW) get quite stagnated following roughly a weak anticyclonic gyre with velocities around 1 cm s − 1 (Pingree 1993; van Aken 2001, 2002; Colas 2003; see Fig. 1). This area is influenced by the extension of the Iberian Poleward Current (IPC), a geostrophically balanced poleward flow, trapped within 50 km of the shelf edge and due to a cross-slope density gradient, which is established in wintertime at the slope of the Iberian margin. The IPC carries warm and salty waters towards the southern waters of the Bay of Biscay (Pingree and Cann 1990; Lavín et al. 1998; García-Soto et al. 2002; Cabanas et al. 2003). Thus, in general, in oceanic waters of the southern Bay of Biscay, advective effects are weak supporting the use of a one-dimensional model.

Fig. 1
figure 1

Eastern North Atlantic basin and location of Santander standard section (black dots) and AGL buoy (black asterisk) in the southern Bay of Biscay. A composition of ENACW circulation (black arrows) in the Bay of Biscay taken from several sources (modified from González-Pola et al. 2005) and IPC (grey arrow) are also included

Full size image

3.2 The Santander standard section

In the southeastern corner of the Bay of Biscay, northwards of Santander, a oceanographic section is running since 1991 as part of the Radiales program from IEO (Spanish Institute of Oceanography; see Fig. 1). This program occupies monthly oceanographic sections around the Iberian Peninsula, mainly located over the continental shelf and sampling hydrographic and biogeochemical parameters (Valdés et al. 2002). Due to the proximity of the shelf break to the coast near Santander, the more external stations of this section are located at oceanic waters. However, meteorological conditions in the Cantabrian Sea make it sometimes difficult to reach outer stations being better covered those closer to the coast.

Station 6 is the best sampled oceanic station with 183 CTD profiles from May 1992 to December 2008, and we will refer to this station hereinafter as IEOS6. IEOS6 is only 6 km apart from the 200-m isobath, and so in some aspects it could be considered as a shelf break station although it is at 850 m depth. Its position at the inner slope is affected by advective processes such as IPC intrusions (Lavín et al. 1998). The next station number 7 (hereinafter IEOS7) is at a depth of 2,400 m. It began in August 1994, existing 130 profiles from May 1992 to December 2008. Thus, because it is our objective to investigate MLD variability in oceanic waters of the southern Bay of Biscay, along this work we will make use of both stations, trying to combine the much better coverage at IEOS6 with the weak influence of shelf-break processes at IEOS7.

On June 2007, an oceanic–meteorological buoy (AGL buoy) was set up in the Santander Section (43°50.67′ N; 3°46.20′ W, 2,850 m water depth and 40 km offshore, www.boya_agl.st.ieo.es). Coinciding with its position, a new station number 9 has been added to the monthly sampling (hereinafter IEOS9). This station is, however, very close to IEOS7. The buoy is equipped with a large set of sensors that register hourly the main atmospheric, hydrological and biogeochemical parameters at the air–sea interface, as well as wave data sensors and an Acoustic Doppler Current Profiler for the measurement of the currents in the first 100 m of the water column. The hourly data and heat fluxes estimation from AGL buoy, together with temperature and salinity profiles measured at IEOS9, provide an accurate local source of forcing fields for comparisons with large-scale meteorological products.

3.3 A climatological seasonality of the upper ocean variability

The quite large and continuous set of CTD profiles series from IEOS6 allowed the development of a method for characterizing the upper ocean variability, which was applied to construct a local climatology of the upper layer vertical structure evolution in oceanic waters of the Santander standard section (González-Pola et al. 2007). The new algorithm developed in that work provides, besides an improved estimation of MLD when compared with other previous algorithms, an accurate representation of the seasonal and permanent thermocline beneath the mixed layer.

This set of climatological profiles will be also used to simulate the upper layer variability in oceanic waters of the Santander standard section to investigate the combined capability of the climatological profiles and true forcing running under the GOTM model to properly describe the upper layer temporal evolution. Positive results will enable us to use them to reproduce past evolution and to speculate on future variability of the upper layer vertical structure under different climatic scenarios in this area.

3.4 Atmospheric forcing

The GOTM model is forced at the sea surface by the heat and momentum fluxes involved in the temperature and U and V momentum equations, respectively, as well as by the precipitation–evaporation balance participating in the salinity equation (see Eqs. 1 and 2). The atmospheric forcing represented by these fluxes is obtained from National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) Reanalysis data set (Kalnay et al. 1996) as six hourly data and with a spatial resolution of 1.9° × 1.9°.

Analyses carried out in order to evaluate the accuracy of the NCEP estimates for turbulent heat fluxes, as well as for the local main meteorological and oceanic variables, against those computed from the AGL buoy and other meteorological data from the Spanish Agency of Meteorology (AEMET) have shown a good agreement (Somavilla et al. 2009; Lavín et al. 2009; Somavilla 2010). The best matching with AGL buoy data is obtained for the closest grid point from NCEP (44.76° N and 2.75° W) to the buoy position than for any multi-grid point average around the buoy. So, NCEP heat and momentum fluxes and precipitation time series used in the different simulations will be extracted at this location.

In order to assess the impact of the accuracy of atmospheric forcing on the results, we will also make use of hourly sensible, latent and momentum fluxes estimated from the AGL buoy measurements following Fairall et al. (1996) formulation (although AGL buoy is currently equipped with long and short wave radiative sensors, there are not available data for the time period of this study and those data have been taken from NCEP/NCAR Reanalysis data set). We will also consider the European Centre of Weather Meteorological Forecast (ECWMF) reanalysis for comparison purposes.

4 Numerical experiments to be carried out

We will outline in this section the sequence of simulation runs that we are going to perform explaining the purpose of each. The specific settings of all runs are summarized in Table 1. The first experiments will take advantage of the relatively extensive available information in recent periods to explore the model performance and gain insight in the dominant processes affecting the upper layer temporal evolution. The comparison of these results with those obtained using the upper layer vertical structure climatology to simulate the upper layer variability will assess us in the way that we must set the model in order to simulate periods where information is scarce.

Table 1 Resume of the study cases and simulations carried out indicating the time period covered and the atmospheric forcing and relaxation used
Full size table

The differences among model runs will be the period of simulation, the sources of air–sea forcings and the background hydrographical field used to relax the water column properties. Relaxing to hydrographical profiles is necessary to account for lateral advection, which affects density distribution and is not resolved in a one-dimensional model. There is lateral advection at time scales of the event variability of meteorological forcing in the area, particularly at IEOS6, that is not resolved by the monthly sampling (Somavilla 2010). More important for the scales of interest in the present study is the meridional advection of central waters that shape thermocline water properties. Equatorwards transport of ENACW formed in wintertime northwards of our area of interest continuously renovates the properties of the permanent thermocline below the mixed layer (van Aken 2001). Therefore, a relaxation of the observed profiles is necessary to represent the effect of advection in thermocline water mass properties in this one-dimensional model.

For the first set of simulations will use the monthly frequency time series of profiles, and they will help us to investigate and isolate the prevailing factors affecting upper waters variability. Three different study cases will be carried out forcing by NCEP/NCAR fluxes and relaxing towards temperature and salinity profiles at IEOS6, the best sampled oceanic station. These study cases, presented in Section 5.1, have been named as 5.1.A, 5.1.B and 5.1.C to systematize the nomenclature. In simulation 5.1.A, the complete surface forcing fields—heat and momentum fluxes and precipitation–evaporation balance—are included. Simulation 5.1.B includes the forcing from heat fluxes and precipitation–evaporation balance but discards the effects of wind stress in the momentum equation and in the shear production of TKE. This simulation aims to isolate the effect of pure convection in mixed layer development. Simulation 5.1.C in contrast includes the forcing from the momentum flux and evaporation–precipitation balance but disregards the contribution of heat fluxes in the temperature equation and their effects in the buoyancy production of TKE. This simulation attempts to isolate the effects of wind shear in mixed layer development. In this respect, we should note that convection and wind stress-driven turbulence occur simultaneously in nature and interact non-linearly, so any ‘ad hoc’ removal of one of the terms may cause inconsistencies (see Rhines and Hakkinen 2003, for a rationale on the problem of surgical removing of terms in a model). However, this exercise can illustrate roughly what are the dominant processes governing MLD variability.

A second set of runs will be executed to investigate possible bias in the simulation results due to uncertainty in atmospheric forcing fields and their temporal resolution. Three simulations will be performed: the first one forced by six hourly data from NCEP/NCAR (simulation 5.2.A) and the second and third by six hourly and hourly air–sea fluxes computed from AGL buoy (simulations 5.2.B and 5.2.C). In both cases, the relaxation is performed towards IEOS9.

The next step will be to evaluate the skill of the model in reproducing the observed interannual variability when using true atmospheric forcings, but relaxing to the climatological water column structure instead to an actual series of profiles. As these climatological profiles cannot account for short-term advective processes like the IPC, the relaxation time scales should increase and this third exercise will also assess on the influence of larger or shorter relaxing time scales. This third set of simulations 5.3.A and 5.3.B share forcings with 5.1.A but will be relaxed respectively towards the seasonal climatology and IEOS6 with a time scale of 6 months.

The last step will be the analysis of mixed layer depth variability in the last 60 years in the Bay of Biscay. Two runs, 6.A and 6.B, will be performed towards the upper waters climatology used in 5.3.A and that climatology modified to incorporate the slow decadal warming in the Bay of Biscay. In general, all simulations are executed in a vertical domain of 500 m with 1 m vertical resolution and a time interval of 1 h.

5 Upper ocean variability and GOTM performance from the monthly oceanic series (1995–2008)

5.1 Seasonal to interannual evolution. Winter MLD variability

As stated in Section 4, three model simulations 5.1.A, 5.1.B and 5.1.C are performed in order to analyse the prevailing factors driving the temporal evolution of the upper layer vertical structure. The origin of simulations is set on January 1, 1995, when a winter profile at the beginning of deep sampling at the Santander section is available. Figure 2 shows the upper layer temporal evolution of temperature and salinity from observations and from model output. A fair agreement between observations and the complete simulation 5.1.A is apparent (Fig. 2a, b, e, f). The seasonal cycle of the mixed layer evolution, consisting of a deepening of the mixed layer between 150 and 200 m depth during the winter and the development of seasonal thermocline at around 30 m depth during the summer, is well reproduced by the GOTM model. The interannual variability is also captured in the model results.

Fig. 2
figure 2

Upper layer evolution from observations and model runs in the southern Bay of Biscay. a, e Temperature and salinity from IEOS6 series of profiles. Black dots represent MLD estimation following the González-Pola et al. (2007) algorithm applied to IEOS6 and IEOS7 temperature profiles. b, f Temperature and salinity from simulation 5.1.A. Black line indicates MLD estimated from GOTM model as the depth of the turbulence surface layer (from simulated turbulent kinetic energy TKE). c Same as b for 5.1.B. d Same as b for 5.1.C

Full size image

The area is affected by seasonal advective patterns that influence the yearly cycle. The main one is the intermittent entrance in wintertime of slope currents associated with the IPC, carrying warm and salty waters to the southeastern Bay of Biscay. These events disrupt and blur the water column structure specially at IEOS6 and that is the cause of including also the MLDs calculated for IEOS7 in Fig. 2 (less affected by IPC intrusions and therefore better representative of oceanic conditions. We expect the simulations to be representative of the oceanic area). Another advective pattern is related to the arrival in summer months of surface low-salinity waters of continental origin (Lavín et al. 1998). These advective sources are accounted for in the model through the relaxation term, so they depend on the existence of a high-frequency time series of profiles. These advective patterns are especially important for the salinity, which is not subjected to a much greater large-scale seasonal cycle as it is the case of temperature. Both in the data and the model output, it is observed the presence of saltier water during winter months in surface waters. The origin of this pattern is the combination of the advective effects and a convective source associated with the mixing of surface waters with saltier deeper waters in wintertime.

Simulations 5.1.B and 5.1.C (temperature fields in Fig. 2c, d, respectively) are used to assess which terms (processes) are dominating in the MLD development. Winter MLDs achieved by 5.1.A and 5.1.B are quite similar among them while 3.1.C underestimates the winter MLD. This exercise just reflects the well-known feature that winter deepening of the mixed layer on mid-latitude areas is mainly driven by convective processes (Kantha and Clayson 2000; Alexander et al. 2000). In short, the heat fluxes in the temperature equation and the buoyancy production of TKE are the principal responsible for the depth reached by winter mixed layer. During winter time, short and long wave radiation heat fluxes almost cancel in the study area so sensible and latent heat fluxes are the only contributors to the highly negative net heat flux in autumn–winter months. On the other hand, wind stress takes on more significance during spring–summer period, causing mixing at shorter time scales that can be observed in simulations 5.1.A and 5.1.C (Fig. 2a, c). During this period, the net heat flux is positive, and so no convective contribution to mixing processes is present until early autumn being generated by wind stress-induced turbulence.

The most outstanding event in these series is the extreme MLDs reached in winters 2005 and 2006, which caused the surfacing of the 12°C isotherm. This event, described in the work of Somavilla et al. (2009), had important consequences for the upper ocean hydrographic structure in the Bay of Biscay. Although described in general for the winter deepening of the mixed layer, the role of convection becomes especially noticeable during winter 2005. Thus, when only convective effects (simulation 5.1.B, Fig. 2c) are taken into account the extreme MLD that occurred in winter 2005 is reproduced, results not achieved when only wind stress effects (simulation 5.1.C, Fig. 2d) are considered, even when the wind stress is increased by a factor of 2 and 3 (results not shown). No significantly different results in MLD are obtained when this increased wind stress is considered together with heat fluxes. This conclusion coincides with that pointed out in Somavilla et al. (2009), where the exceptional cold conditions that caused huge buoyancy losses (convection) and that only had precedents during the 1960s were considered responsible for the unprecedented depths of winter mixed layer. The following winter 2006 reproduced a very deep mixed layer due to what is known as the re-emergence mechanism (Timlin et al. 2002), a process by which the low-stratification anomaly consolidated during a hard winter gets stored year-round beneath the summer surface thermocline, thus facilitating the advance of the newly formed mixed layer in the following winter. The ultimate cause behind such a deep mixed layer in year 2005 seems to be linked to a couple of specific episodes of very strong buoyancy loss at the end of January and February. The second one triggered a sudden increase in MLD from ~150 m on February 25 to more than 330 m on March 8. This sequence of events is fairly well reproduced by the daily model runs of GOTM model.

5.2 Daily to seasonal surface variability. Reliability of atmospheric forcing fields

Figure 3 shows a comparison of sea surface temperature (SST) and sea surface salinity (SSS) of simulation 5.1.A and IEOS6 data. Surface temperature variability is reasonably well represented, being the mean difference \(\lesssim \)0.5°C. Larger differences occur during spring and summer, peaking above 2.5°C in some sporadic events and reaching even 4.5°C in spring 1998. These differences disappear quickly at depth, pointing to uncertainty in atmospheric forcing fields used to run the model as their cause. During summer months, an eventual overestimation of net heat flux would act on a very shallow mixed layer, leading to an excessively large SST increase. Even so, the main anomalies along the series are well reproduced by modelled SST and SSS. We can highlight the cold summers of 2002 and 2007, the abnormally warm summers of 2003 and 2006 (Luterbacher et al. 2004, 2007), the extremely cold winter of 2005 (Somavilla et al. 2009), the warm autumn–winter periods of 1998 (Levitus et al. 2000) and 2006–2007 (Luterbacher et al. 2007) etc.

Fig. 3
figure 3

a Sea surface temperature and c salinity from observations at IEOS6 (dot solid grey line) and from simulation 5.1.A (solid black line). b, d Difference between modelled and observed SST and SSS, respectively, at IEOS6 time series sampling moments (dot in dot solid grey line)

Full size image

As pointed in Section 4, we take advantage of the availability of data from the ODAS buoy AGL to investigate the impact of changing the atmospheric forcing fields and their temporal resolution. Simulations 5.2.A, 5.2.B and 5.2.C are designed to discriminate among NCEP vs. AGL forcing fields (see Table 1). Note that the depth of the water column is set to 300 m, the maximum depth at IEOS9 sampling. In Fig. 4 MLD, SST and SSS for these runs are shown. Simulations show good agreement among them. Differences between simulation 5.2.A (NCEP heat fluxes six hourly) and simulations 5.2.B and 5.2.C (heat fluxes from AGL) are apparent in some events in winter, when 5.2.A simulates a deeper MLD. Local heat fluxes computed using local SST (see Section 3.4) like those in 5.2.B and 5.2.C result in shallower MLD, although the temporal variability of deepening events is similar. Wind differences between simulations (NCEP vs. local, hourly vs. six hourly) do not result in clear MLD variations, although some differences in SST are apparent in spring and early summer 2008, when wind stress (typically affected by a higher frequency temporal and spatial variability) is the main contributor to mixing processes. The increase of the temporal resolution of atmospheric forcing from six hourly to hourly data does not result in remarkable differences between both AGL runs (5.2.B and 5.2.C in Fig. 4).

Fig. 4
figure 4

a MLD from simulation 5.2.A (dash grey line), 5.2.B (solid grey line) and 5.2.C (solid green line) from July 2007 to December 2009. b Sea surface temperature and d salinity from observations at IEOS9 (dot solid red line) and AGL buoy daily mean values (black line) and from simulation 5.2.A (dash grey line), 5.2.B (solid grey line) and 5.2.C (solid green line). c, e Difference between modelled and AGL buoy daily means of SST and SSS, respectively, for simulation 5.2.A (dash grey line) and 5.2.B (solid grey line). The difference found in Section 5.1 between IEOS6 time series and simulation 5.1.A is also shown (dot solid grey line). The dot solid red line in c and e represents the difference between IEOS9 and AGL buoy daily mean of SST and SSS at IEOS9 sampling moments

Full size image

In terms of SST, the mean difference between modelled and AGL buoy SST daily mean time series \(\lesssim \)0.5°C (Fig. 4c) is similar to the difference found previously between simulation 5.1.A and IEOS6. There is also a good agreement between SST data from IEOS9 at the sampling snapshots and AGL or modelled SST (dots in dot solid red line in Fig. 4c). However, we should notice that the direct comparison against IEOS9 daily interpolated values (as is mandatory whenever only monthly data is available) would yield much larger discrepancies (Fig. 4b). Thus, the inefficiency of IEOS9 to resolve higher frequencies of variability is a source of error in SST.

In relation with this fact, the difference found between simulation 5.1.A and IEOS6 is also shown in Fig. 4c, e. There are not large SST difference excursions within these periods July 2007 to December 2008, so it is not possible to attribute such peaks to errors in the atmospheric forcing fields. Extra runs have been performed to see whether finer vertical resolution (0.5 m) may improve the performance of the GOTM model, although this yields occasionally temperatures up to 1.5°C colder than the 1-m resolution standard runs, and the norm is to observe a similar overestimation of SST.

As expected, divergence between modelled and observed SSS is greater than for SST. Mesoscale structures associated with the advection of low salinity waters during summer months (Lavín et al. 1998) have a quite high-frequency variability, so they cannot be accounted for properly through the relaxing term. During that time together with the presence of superficial low salinity waters, a shallow halocline is observed at IEOS9 profiles that hinders mixing processes and the entrainment of cold waters from below (Somavilla 2010). The entire effect of this advection is poorly represented in the one-dimensional simulations causing partly the subestimation of SST in AGL buoy runs (5.2.B and 5.2.C) observed in Fig. 4 during summer months. In general, the one-dimensional modelling of vertical mixing processes can only obtain an approximation of the SSS seasonal cycle.

5.3 Reproduction of the interannual variability from the climatological seasonality

The main goal of this section is to check whether the upper ocean variability from air–sea forcings (either measured or derived from a reanalysis) can be reproduced by the GOTM model relaxing towards a climatological (average) seasonal evolution of the vertical ocean profiles. Figure 5 shows the temporal evolution of temperature and salinity in the upper layers obtained from simulations 5.3.A and 5.3.B (6-month relaxation towards climatological series of profiles and to IEOS6, respectively). Thus, simulation 5.3.B differs from 5.1.A in the relaxation time scale (6 and 1 month, respectively). It allows, as a previous step, the evaluation of the role of the relaxation time scale in the model when using an actual and climatological series of profiles. In the same figure, it is also displayed the temperature and salinity at 100 and 200 dbar from IEOS6, IEOS7, model runs 5.1.A, 5.3.A, 5.3.B and also an extra model run without relaxation (Fig. 5d–f). The effect of the relaxation term in the simulation is evident from this no-relaxation run. In such case, the only information on the actual stratification of the water column is the initial profile of the simulation (January 1995). Every winter, the development of the mixed layer erodes the initial profile, and the thermocline gets cooler and deeper every year and MLD reaches deeper values progressively. Thus, the relaxation in the model is needed to preserve realistic thermocline properties (hence the stratification) that, in the Bay of Biscay, is continuously modulated by northern large-scale advection.

Fig. 5
figure 5

Upper layer temporal evolution of temperature from simulations 5.3.A (a) and 5.3.B (b) and salinity (c). The black line indicates the MLD. d Temperature at 100 dbar from model simulation 5.1.A (solid black line), 5.3.A (solid brown line), 5.3.B (solid red line), no relaxation (magenta line), IEOS6 (dot solid blue line) and IEOS7 time series (dot solid green line). e, f Idem for temperature and salinity at 200 dbar

Full size image

Simulation 5.3.B (Fig. 5b) yields roughly a similar pattern of seasonal and interannual variability of MLD than simulation 5.1.A (Fig. 6a, c), but some interesting differences arise. As time relaxation increases, the strength of the advective component that maintains the permanent thermocline diminishes. When a 6-month relaxation scale is used, waters tend to be fresher and colder (Fig. 5d–f). This feature is related to the reduction of the effect of the advection of warmer and saltier waters from the IPC pulses induced by the increase in time relaxation. Simulation 5.3.B do show the interannual variability determined by IPC-transported waters and also the intense IPC years are evident (Fig. 5b, c), but events are more apparent in 5.1.A (Fig. 2b, e) when w e relax tightly. Mixed layer depths are also slightly affected, being for instance deeper those obtained in 1998 and 2003 by simulation 5.3.B (Fig. 6a, c). The biggest difference in winter MLD among both exercises is found in 2006. Simulation 5.3.B provides a much deeper winter 2006 MLD than 5.1.A does, even exceeding that of 2005 (Fig. 6a, c). The very low stratification achieved in winter 2005 below the seasonal thermocline is over preserved year-round with a longer relaxation term, thus easing the advance of the next year MLD (re-emergence). In agreement, temperature time series at 100 and 200 dbar shows colder values for simulation 5.3.B during winter 2006 (Fig. 5b, d, e). Actually, we cannot discard the presence of deeper MLD in 2006 than in 2005 in the region. Although we did not capture such signal in our standard section, this may have been missed from a monthly series.

Fig. 6
figure 6

a Mixed layer development computed from model simulation 5.1.A (solid black line) and González-Pola et al. (2007) algorithm applied to IEOS6 and IEOS7 time series (red dots). b Idem for MLD computed from model simulation 5.3.A and c MLD computed from model simulation 5.3.B

Full size image

Simulation 5.3.A (Fig. 5a, c), relaxed towards climatological temperature profiles and salinity series of observations, provides essentially a similar pattern to the initial simulation 5.1.A (similar years of deep and shallow MLD, see Fig. 6a, b). This result is of paramount importance. It confirms the pre-eminence of local air–sea interaction as primary driver of winter MLD variability, indicating that a climatological treatment of the large-scale lateral advection allows the GOTM model to achieve a reliable water column evolution. In particular, for the extreme winter mixing of 2005, both simulations diagnose a similar MLD. As expected, differences arise in relationship with the effect of IPC intrusions in MLD detailed in the previous paragraph, an effect that is in no way accounted for when relaxing to climatological profiles. Good examples are the winters 1998, 1999 and 2007 that present a deeper MLD in 5.3.A (Fig. 6b respect to Fig. 6a). The absence of warmer and saltier waters advected by the IPC is observed through the colder and fresher water column in simulation 5.3.A than for simulation 5.1.A in those years (Fig. 5d–f).

The re-emergence event in 2006 supposes a challenge for our climatological approximation. Advected waters reaching the point of interest in 2006 were actually affected by the same large-scale atmospheric and mixed layer temperature anomaly in 2005 (see Figs. 6 and 7 in the work of Somavilla et al. 2009) and, as we know, the low-stratification anomaly still persists on its arrival. By relaxing to climatological profiles, we are incorrectly forcing the model towards a climatological state. However, simulation 5.3.A does show a re-emergence signature replicating a deep mixed layer in 2006. It seems that setting a 6-month time scale we are relaxing slow enough to keep locally year-round part of the anomaly and reproduce (partially) the re-emergence event. It is observed how temperatures at depth drifts just after winter 2005 in simulation 5.3.A (Fig. 5d, e), but the event in 2006 is not completely lost, and the coldest temperatures at 100 and 200 dbar are actually reached during the winter of 2005 and 2006.

In summary, simulation 5.3.A is able to reproduce the seasonal cycle and the most conspicuous events observed at IEOS6 time series, which are direct consequence of atmospheric forcing. At the same time, short-term advective effects on the MLD like the IPC entrance are neglected. In this sense, if we peer into isolated signatures of warm and salty water (e.g. 1998, 2002 and 2007), we observe that simulations 5.3.A and 5.3.B are more similar to IEOS7 data than to IEOS6 (Fig. 5d–f). As said before, IEOS6 at the shelf-break is much more affected by shelf-slope advective processes than IEOS7 (unfortunately the later station is much poorly sampled and several-month gaps complicate its use for relaxation purposes). Therefore, if we are interested in studying the mixed layer along an oceanic large-scale region, such avoidance of shelf-break advective anomalies is a plus for the use of a climatological series for relaxing purposes. Finally, large-scale advection of anomalies like the 2006 re-emergence gets partially simulated with climatological series through the retention of the anomaly locally. The reproduction of such re-emergence phenomena depends of course on an appropriate election of relaxation time scale and requires the anomaly generated locally to be equivalent to that generated upstream and then advected.

In the next section, we will attempt to study past conditions of the upper layer vertical structure using the GOTM model. We have seen that we need to relax to a set of temperature and salinity profiles in order to account for the effect of large-scale lateral advection in thermocline water properties and stratification. We must be aware that historical observations are sparse with long gaps and information about large-scale ocean advection variability is scarce. Any interpolated series constructed from such scattered data, which even skips whole winters without having a single profile, would add an important bias to the model runs (early occurrence of mixed layer development and seasonal thermocline and important variations in the final MLD reached at each time step). It can result in important errors and severe consequences when attempts are made to couple physical model results with some other modules such as biogeochemical ones (Friedrichs et al. 2006). Within this section, we have demonstrated that the relaxation to a climatological series of profiles, carefully constructed from a large data set, results in a fair reconstruction of the upper water variability in our region. We will therefore combine the use of this climatology, information about past states and the GOTM model with the aim of reconstructing the historical water structure evolution in our region.

6 Upper ocean variability in the long term (1949–2008)

6.1 Results of the long-term run. Constrains and reliability

The availability of NCEP/NCAR reanalysis since 1948, together with the success achieved in the use of climatological profiles for estimating water column structure, allows us to carry out a simulation of the upper layer temporal evolution for the last 60 years in the Bay of Biscay. The simulation was initiated on 1 January 1948 (beginning of NCEP/NCAR Reanalysis data set) and finished on 1 October 2008 with a time relaxation of 6 months towards the climatological profiles. As the remaining simulations throughout this work, it was executed for a vertical domain of 500 m with 1 m vertical resolution and a time step of 1 h.

Figure 7a presents the resulting MLD time series from 1948 to 2008. This simulation yields two remarkable features. Firstly, during the 1970s and 1980s, MLDs were strikingly and unexpectedly shallower than MLDs from 1995 onwards. Secondly, during the 1950s and 1960, some years of very deep MLD are apparent, namely 1956, 1963 and 1965, but none of them as deep as that in 2005. Nevertheless, accepting this deepening of MLD in recent times—including a record-breaking event—in an area with a strong upper ocean warming trend (González-Pola et al. 2005; Castro et al. 2009) is, at least, quite counter-intuitive. Therefore, before proceeding to analyze these results, we will discuss further whether the setting and forcings of the model may be distorting its overall behaviour.

Fig. 7
figure 7

a Mixed layer development computed from GOTM model simulation 6.A for the last 60 years in the Bay of Biscay. b Winter mean net heat loss and momentum flux (c) from 1948 to 2008

Full size image

The response in MLD seems to be strongly linked to net heat flux (Fig. 7b). Winter mean net heat losses shifted from ~90 W m − 2 experienced during 1970s and 1980s to ~105 W m − 2 after 1995, and this feature causes such apparent mixed layer deepening. Oceans are heat reservoirs, acting like a flywheel in matters related to weather on time scales of weeks and longer. They gain heat during the spring and summer losing it slowly during fall and winter. Winter net heat loss increase observed since mid-1990s is ultimately a combination of variability of storminess activity and air temperature during winter months acting combined with excess of heat gained in summertime. Therefore, the first question is about the reliability of the forcing fields. The results of Section 5.2 indicate that MLD show limited sensitivity to variation of atmospheric forcing products for a modern period. However, historical reanalysis data may have some bias, and for this reason, NCEP/NCAR fluxes were compared with those obtained from ECWMF Reanalysis data set. No significant differences between both were found, so a potential bias in NCEP/NCAR atmospheric forcing would be shared by ECWMF. On the other hand, both reanalysis are affected by a similar coarseness that may result in limitations when small-scale physical processes are important. Nevertheless, for studies focused in interannual variability, the air–sea forcings have a large-scale character, thus generating similar effects in the upper layers in broad areas.

Apart from atmospheric forcing, another factor that could result in unrealistic behaviour of the model is the water column structure imposed through relaxation to account for unresolved large-scale lateral advection. The climatology used so far was generated with CTD profiles obtained during the 1990s and onwards at IEOS6. It represents a warmer water column than that in existence during the previous decades in this area, and stratification may have also varied in an unknown way. In principle, one could suspect that extracting the ‘true’ NCEP heat fluxes from this extra-heated water column may result in an underestimation of MLDs in previous decades. To overcome this issue, the climatology was modified to incorporate the slow decadal warming in the Bay of Biscay. The procedure was to calculate the winter SST anomaly between the climatology and the NOAA_ER SST (Smith et al. 2008) and then using a low-passed version of such anomaly to modulate the whole water column. We are aware of the caveats inherent to this approach, but we should not forget that the climatology is only intended for accounting for missing the large-scale lateral advection with a time scale of several months, and the interannual variability should be determined by the short-scale resolving model run. The purpose of this experiment is to avoid a systematic bias towards a too warm water column for decades.

This new model run is presented together with the previous one in Fig. 8a, b. Although some small differences between both simulations are observed in terms of absolute MLD, the same long-term variability is obtained, i.e. (a) large MLDs in 1963 and 1965 are still shallower than in 2005 and (b) MLDs in the 1970s and 1980s are in general shallower than those from 1995 onwards. This confirms once again the dominance of net heat fluxes in the determination of the MLD.

Fig. 8
figure 8

a, b Upper layer temperature evolution from 1948 to 2008 from simulation 6.A and 6.B, respectively. The black line indicates the MLD estimated from GOTM model simulation. c Winter minimum SST from NOAA_ER SST (Smith et al. 2008) (dot solid black line) and satellite data from NOAA_OI_SST_V2 (Reynolds et al. 2002) (dot solid green line) and from GOTM model simulation 6.A (dotted red line) and 6.B (solid red line). d 200-300-m-depth temperature time series from available historical data in the area 43–48° N; 008–003° W (dot solid black line), IEOS7 (green dots), simulation 6.A (dotted red line) and simulation 6.B (solid red line). e NAO index from 1950 to 2009 (black line). f AEMET winter mean air temperature from 1961 to 2007 at Santander station. The red line in e and f represents the 5 years running mean of the time series

Full size image

This new run can be compared with the available data sets at the area, which are at the surface the aforementioned NOAA_ER SST used to modulate the climatology and the satellite SST data from weekly optimum interpolated NOAA_OI_SST_V2 (Reynolds et al. 2002) and at depth the historical dataset of hydrographical profiles at the Bay of Biscay (sources are WODB2005 global database (Boyer et al. 2006) and SISMER French data set www.ifremer.fr/sismer/ among others) and IEOS7 data. Although historical profiles are scarce and scattered, we can construct a reliable temperature series focusing on the fairly stable level at the upper permanent thermocline that is eventually hit by the mixed layer advance in late winter (~200–300 m). Modelled and observed evolution of winter minimum SST and temperature at 200–300 m depth is shown in Fig. 8c, d. There is a tight agreement among trends from both the data and the model run at depth, indicating a reliable reproduction of water properties at the upper thermocline. Model run 6.B almost corrects the systematical deviation towards warm winter SST values in the 1970s and 1980s in simulation 6.A compared with what is observed in NOAA_ER SST, but some deviation still persists.

This discrepancy raises questions on the overall consistency of the long-term model run. As the model works fine within the time domain for which we have the series of profiles, it is difficult to assess the cause of these warm mixed layers in the 1970s. There are, however, good reasons for thinking that MLD has deepen in recent times.

Similar warming trends are observed both at the surface and at depth in Fig. 8c, d (time series at both depths evolve in parallel) in observations and modelled data (run 6.B). However, this pattern changes in 1995 and a ‘stabilization’ of winter minimum SST is evident both in observed and modelled data while at depth the warming trend continues until 2005 (extreme winter mixing event) at a rate of 0.234°C · decade − 1 in IEOS7 data and 0.218°C · decade − 1 in historical hydrographical data. It supports the counter-intuitive outcome of shallower winter mixed layers during some previous periods colder than today.

A simplified view of the process is that the deeper level is warming faster than surface waters which implies higher temperature differences between both levels during the 1970s and 1980s. In that case, equal net heat losses would produce deeper MLDs at the end of the time series. Besides if a winter net heat loss increase exists, as it is the case here, this tendency would be reinforced. Considering that the net heat flux across the air–sea interface, Q 0, was perfectly balanced with the heat storage in the water column, Q strg, that is:

$$ \rho {\kern 1pt} C_{p{\kern 1pt} } {\kern 1pt} h{\kern 1pt} \frac{{\partial T}}{{\partial t}}{\kern 1pt} = {\kern 1pt} Q_0 {\kern 1pt} $$
(3)

where ρ is the sea water density and C p is the heat capacity of sea water, we could estimate the mean depth reached by the winter mixed layer, h. Taking into account the observed temperature differences between the surface and 200–300 m depth (observed winter mean SST − \({\rm Tem}_{200\text{--}300\,{\rm m}}\)) before and after 1995 (~0.9°C and ~1°C, respectively) and the winter mean net heat loss for both periods (see above), a mean winter deepening of the mixed layer of ~190 and ~237 m would be obtained for each period, respectively. Although smoothed, the differences in the MLD obtained between both periods would be maintained based on observations.

Finally, further experiments with model SST restoring to SST observations result in deeper MLDs in the whole period and in a smaller difference between MLD in the 1970s and 1980s and in recent years. However, the observed decadal difference in MLD, although smaller, is also still present in those experiments.

6.2 Characteristics of extreme MLD years

The long-term run in simulation 6.A,B indicates the existence of large winter MLDs in the 1950s and 1960s. Actually, Somavilla et al. (2009) suggested, based on the similarity in atmospheric pattern variability in the North Atlantic, that similar or even deeper MLDs to those found in 2005 in the Bay of Biscay may have taken place during the winters of 1963 and 1965. Unfortunately, there are no winter profiles in the area to confirm or rule out such hypothesis. Nevertheless, the present simulations provide a tool to investigate further these events. In principle, simulated MLDs in 1963 and 1965 are shallower than those observed in 2005, and contrary to the suspicious period of the 1970s, modelled SST does not disagree much with respect to NOAA_ER data.

Figure 9 compares the succession of mixing events in terms of heat loss and MLD evolutions during the winters of 1963, 1965 and 2005. The different sequences of mixing events show themselves as a remarkable difference in these years. In 1965, frequent mixing events are observed evenly distributed throughout the winter while 1963 and 2005 are characterized by episodes of great heat loss. In particular, 1963 is highly conditioned by an extreme heat loss mixing event on January 13 (761 W m − 2) that cause the MLD to shift abruptly from 70 to more than 150 m. After that, mixed layer development continues until mid-February, alternating stabilization and moderate mixing episodes to achieve a maximum depth of 250 m. In 2005, two episodes of extreme heat loss occurred at the end of January and February. The second one was as large as that observed in 1963 (764 W m − 2) causing the mixed layer to progress from a typical winter mean value (150 m) to the 330-m depth that were luckily observed at Santander section temperature profiles (IEOS7).

Fig. 9
figure 9

a MLD (black line) from simulation 6.B and net heat loss during winter 1963 (solid grey line). b Idem for winter 1965 and c winter 2005. Red dots in d represent MLD estimation following the González-Pola et al. (2007) algorithm applied to IEOS6 and IEOS7 temperature profiles

Full size image

Although winter mixing and convection processes are mostly governed by the accumulated effect of heat loss throughout the winter, extraordinary convection episodes appear more likely to occur when heat loss is concentrated in intense mixing events, rather than distributed evenly but more weakly throughout the winter. This conclusion has also been pointed in other works (Marshall and Schott 1999; Vage et al. 2008). It is also worthy of mention, as previously discussed in Sections 5.1 and 5.2, that wind stress-induced turbulence does not play a fundamental role in the generation of these extreme mixing events. Thus, the highest wind stress was registered in winter 1965 followed by 1963 and 2005. On the contrary, it is the existence of strong air-sea temperature and relative humidity gradients involved in sensible and latent heat flux equations respectively what determines strong heat loss to occur during winter months.

6.3 Low-frequency variability in MLD and large-scale atmospheric patterns

Low-frequency and interannual variability in the extent of the winter MLD influence a number of key oceanographic processes, from modal water masses formation to primary production processes. Determining such variability from field data requires extensive data sets with enough coverage at the critical end-winter period, so there are actually few works in this field. Prominent examples of description of long-term trends in MLD are Polovina et al. (1995) and Freeland et al. (1997) for the North Pacific and Michaels and Knap (1996), Nilsen and Falck (2006) and Carton et al. (2008) for North Atlantic. The ultimate drivers of these low-frequency changes in the MLD are large-scale atmospheric forcing patterns.

Although our exercise has been performed with a one-dimensional model, subjected to uncertainties in the representation of the structure of the water column, the model runs and the relationship among trends in SST and in the upper permanent thermocline temperature suggest a pattern of shallower MLDs during the 1970s and 1980s and deeper ones since mid-1990s for the southern Bay of Biscay. Taking into account the MLD cycles described for the North Atlantic and recent reviews of the effects of the North Atlantic Oscillation (NAO) on ocean and air–sea fields (Marshall et al. 2001; Visbeck et al. 2003; Hurrell and Deser 2009), we will discuss in this section the low-frequency MLD oscillation inferred in oceanic waters of the Santander standard section and its differences and similarities with other regions in the North Atlantic.

The patterns of air–sea interaction associated with the NAO have a robust spatial pattern, and so it has been shown in several works for air temperature, SST, heat fluxes, wind stress or precipitation (Marshall et al. 2001; Visbeck et al. 2003, and references therein). Net heat flux anomalies associated with high and low NAO index years show a four-lobe structure which is also observed in maps of correlation between winter SST and NAO (Fig. 10a, b). The feature is commonly known as ‘NAO SST tripole’ (Visbeck et al. 2003, and references therein) because one anomaly is centered in the western North Atlantic between 30◦ N and 45◦ N and another two anomalies of opposite sign are centered in the subpolar gyre and in the northern tropical Atlantic. This tripole pattern explains the counterphase oscillations in mixed layer depth described at Hydrostation S in the Sargasso Sea (Michaels and Knap 1996; Paiva and Chassignet 2002) and the subpolar gyre (Carton et al. 2008; Henson et al. 2009, and references therein). A fourth lobe, located at the eastern Atlantic margin between 40° N and 60° N, seems to present an analogous behaviour to that observed in the western part (i.e. being in phase with the temperate western north Atlantic lobe and so opposite to that found in the subpolar gyre). This fourth centre of anomaly is of special interest in our case as it extends to the Bay of Biscay.

Fig. 10
figure 10

Net heat flux anomaly (watts per square metre) during high NAO positive (a) and negative (b) index years

Full size image

The NAO index underwent a strongly positive period in the 1970s and 1980s to reverse from mid-1990s onwards (see Fig. 8e). Ascribed to positive NAO index periods, it is expected in the subpolar gyre an increase of storm frequency concurrent with a positive anomaly in winter net heat loss (higher than the mean value, Fig. 10a) that has been related to the MLD deepening trend and SST decrease observed during the 1970s and 1980s in this area (Carton et al. 2008; Henson et al. 2009, and references therein). On the contrary, at the east and west lobes (Fig. 10a), it is expected a decrease in the passage of winter storms and a negative anomaly in winter net heat loss (lower than the mean value). Therefore, a shallowing of MLDs and a winter SST increase would be expected to occur in these areas. At the southern Bay of Biscay, winter net heat loss provided by NCEP/NCAR does indeed show these low anomaly of heat-released (Fig. 7b) that in turn caused the GOTM model to develop warm and shallow MLDs (although as previously discussed the effect appears to be exacerbated; Fig. 8b). At the eastern lobe, similar results (shallower MLD) have been reported at Hydrostation S in the Sargasso Sea for that period (Michaels and Knap 1996; Paiva and Chassignet 2002).

As NAO enters the opposite phase from 1995 onwards, the tendencies all around the North Atlantic reverse, i.e. positive [negative] anomaly in winter net heat loss at mid-latitude lobes [subpolar and subtropical lobes] (Fig. 10b) causes increase [decrease] of MLD. It also coincides with that reported in the works of Carton et al. (2008) and Henson et al. (2009) for the subpolar gyre; Michaels and Knap (1996) and Paiva and Chassignet (2002) for the western Atlantic margin; and found in this simulation for the Bay of Biscay.

Thus, air–sea forcing fields in the North Atlantic and its relation with NAO index decadal variability explain long-term variability in MLD time series all around the North Atlantic and including the Bay of Biscay. Due to the fact that the entire North Atlantic basin is affected by this low-frequency variability which also has effect in North Atlantic circulation (Curry and McCartney 2001), a change in both the thermohaline properties and the velocity of the advected waters maintaining the main thermocline in the study area cannot be rejected before and after 1995 (see Fig. 12 in Curry and McCartney (2001) showing lower and higher mass transport across the North Atlantic for positive and negative NAO index years). Such variation (increase) may be considered as a possible contributor for the warm MLDs obtained during the 1970s and 1980s in GOTM model runs. Moreover, this interannual variability is nowadays superimposed on a generalized global warming that influences temperature at all levels in a way not totally understood. As long as global warming does not result in geographically and seasonally homogeneous warming, its effect on this oscillatory pattern of variability may differ among regions.

Other low-frequency atmospheric modes of variability different from NAO have also a role in determining the winter MLD. The Eastern Atlantic (EA) pattern has been already suggested as an important driver in this area (Somavilla et al. 2009). In its negative phase, a high sea level pressure anomaly is located west of the British Isles, inducing strong cold easterly winds which enhances the heat loss over the Bay of Biscay. This atmospheric feature was the ultimate cause of extreme mixed layer depths in winter 2005 and 1965. A high negative winter EA is also observed in years 1956, 1971, 1976 and 1981, when deeper MLDs than those observed in their neighbouring years are found, mainly in 1956 and 1971, see Fig. 7a. In 1976 and 1981 an even distribution of net heat losses throughout the winter is observed, and as has been explained previously (Section 6.2), such uniform pattern implies a lower capacity to induce extraordinary convection episodes than other configurations exhibiting intense mixing events. Interestingly, a strengthening and a northeastward displacement of the Iceland low and Azores high by 10–20° in latitude and 30–40° in longitude has been described in a global warming scenario (Hu and Wu 2004). The resulting high level pressure close to a position west of the British Isles resembles in some way the negative phase of EA, which could have important consequences in the future configuration of upper layer vertical structure in this area.

It is interesting to note that, due to the similarity in the large-scale atmospheric forcing, re-emergence phenomena could have taken place in the years following the winters of 1963, 1965, 1956 and 1971, as occurred in winter 2006 after the extreme cold winter of 2005. Simulation results could be partially underestimating MLD in 1964, 1966, 1957 and 1972.

7 Conclusions

The combination of a long-term hydrographic time series running since 1991 in the southern Bay of Biscay, a climatology of upper layer temporal evolution obtained from this time series and the GOTM model to study the upper layer temporal evolution in the southern Bay of Biscay has allowed an extensive study of the prevailing forces governing the seasonal, interannual and long-term MLD variability. Different simulations have been carried out obtaining a high level of detail of MLD variability and extreme phenomena observed at Santander standard section and good agreement with observations.

As expected, the different simulations have shown that on seasonal time scales, winter mixed layer development is mostly conducted by convection processes while wind stress is responsible for mixing events during the spring–summer season. The role of convection processes is proved to be especially noticeable during the exceptional winter mixing occurred in winter 2005, coinciding with that noted in Somavilla et al. (2009). The sequence of events is also well simulated within a single year MLD development, and the re-emergence phenomenon in winter 2006 responsible for the easy reproduction of a new deep and cold mixed layer is reproduced with different levels of success by the GOTM model in the different simulations.

On the other hand, the comparison of the results of the different simulations and the observed variability has shown that, although atmospheric forcing is mainly responsible for MLD variability and the extraordinary phenomena observed in oceanic waters of the Santander standard section, advective processes such IPC whose intensity diminishes towards oceanic waters can determine slight differences in the final depth reached by winter mixed layer.

Tests carried out to investigate climatological profile skill together with true forcing to reproduce the upper layer temporal evolution in the southern Bay of Biscay have demonstrated its ability to provide a fair reconstruction of the upper water variability in our region. They allow to account for the effect of large-scale lateral advection in thermocline water properties and stratification while avoiding the inclusion of shelf-break advective anomalies which is a plus for their use for relaxing purposes in studying the mixed layer along an oceanic large-scale region.

It has permitted a first trial of the reproduction of the last 60 years of MLD variability in the study area. Remarkable results have been obtained. Although in certain grade overestimated by simulation results, an unexpected period of shallower MLDs seems to have occurred during the 1970s and 1980s which would have been concluded from mid-1990s with a period of deepening trend in MLD. The reproduction of sea surface and 200 m depth temperature time series and the warming trend at both levels supports the counter-intuitive outcome of shallower winter mixed layers concurrent with generalized upper water warming trends reported on several occasions for the area.

As found in other recent studies, long-term trends in MLD in the southern Bay of Biscay seem to be related to decadal variability in the NAO, being in phase and opposition to other cycles of deepening and shallowing trends in MLD found from subtropical-to-subpolar areas in the North Atlantic. Favourable sequences of mixing events result in intense convection processes explaining interannual differences and the extraordinary deepening of the winter mixed layer as occurred in the years 2005, 1963 and 1965. The negative phase of the Eastern Atlantic pattern seems to determine important interannual variability through intense episodes of cooling and mixing as occurred in the winter of 2005 in the Bay of Biscay.