⁎Corresponding author.
E-mail address:antti.westerlund@fmi.fi(A. Westerlund).
Contents lists available atScienceDirect
Journal of Marine Systems
j o u r n a l h o m e p a g e :w w w . e l s e v i e r . c o m / l o c a t e / j m a r s y s
and on the river runoff. The sill in the Åland Sea towards the Baltic Prop-er is shallow (c. 70 m), thus only the watProp-er above the halocline (with sa-linity of c. 6.5–7 PSU) penetrates to the Bothnian Sea forming the bottom water of the basin. The halocline in the Bothnian Sea is at the depth of 60–80 m. The surface salinity varies between 4–6 PSU decreas-ing towards the north. Overall, the stratification in the Bothnian Sea is rather weak and mixing is able to penetrate also to deeper layers. A more detailed description of the Baltic Sea and its hydrography can be found e.g. inLeppäranta and Myrberg (2009), and of the Gulf of Bothnia specifically e.g. inHåkansson et al. (1996).
In many ways, the Gulf of Bothnia is one of the least investigated sub-basins of the Baltic Sea, as noted by e.g. Omstedt and Axell (2003). It has had good ecological status and extended monitoring has not been a key priority unlike in other, more eutrophicated areas of the Baltic Sea. There has been a renewed interest in the Gulf of Bothnia, as ecosystem health indicators and eutrophication-related parameters have shown a decline in the eutrophication status of the area (HELCOM, 2014; Fleming-Lehtinen et al., 2015; Lundberg et al., 2009).
Warming trends in surface air and sea temperatures seem stronger in the Gulf of Bothnia than elsewhere in the Baltic Sea (BACC II Author Team, 2015; Lehmann et al., 2011), but relatively large uncertainties re-main in terms of the impact of climate change on the biogeochemistry of the area (e.g.Meier et al., 2012). In addition to monitoring, evalua-tions based on model studies are also needed. Unfortunately, compari-sons of physical-biogeochemical models have found that typical model performance in the area is poor. There are problems in the repre-sentation of both biogeochemistry and physics (Eilola et al., 2011).
One notable factor that has previously limited process and modelling studies of the northernmost Baltic Sea has been the limited amount of observational data. In Baltic Sea off-shore areas, temperature and
salin-models, observations need to have sufficient spatial and temporal cov-erage. New measurement techniques that are already frequently used in the oceans, such as Argofloats and gliders, give the opportunity to collect new data for model validation and development. The Finnish Meteorological Institute (FMI) together with Aalto University have been developing Argofloats that can also be operated in shallow seas (Purokoski et al., 2013). Thefirst tests were made in 2011 and since then Argofloats have been operatedfirst in the Bothnian Sea and later also in the Baltic Proper. This new dataset is well suited for evaluating the capability of hydrodynamic models to produce the vertical structure of temperature and salinity. It provides a time series of profiles from the area of interest with good temporal resolution, showing the dynamics of temperature and salinity in the water column throughout the summer.
State-of-the-art 3D ocean models, such as NEMO (Madec and the NEMO team, 2008), provide a good basis for implementing and further developing model applications for coastal and shelf seas. However, sim-ilar to other models originally developed for oceans, several adaptations are required, including adjusting bottom friction and tuning the turbu-lence schemes. A specific configuration of NEMO, named NEMO Nordic, has been made for the Baltic Sea studies byHordoir et al. (2013a,b, 2015). This configuration has been previously used to study e.g. air– sea coupling (Gröger et al., 2015) and degradation of dissolved organic matter (Fransner et al., 2015). By studying the performance of this model configuration in the Bothnian Sea, we can gain further under-standing of its capabilities and also of any improvements that may be re-quired. We focus our studies on the vertical structure of temperature in the Bothnian Sea, from which new observations from autonomous Bal-tic Sea Argofloats are available. With this new data we can evaluate the ability of 3D models to produce vertical temperature profiles and mixed layer depths in the area with unprecedented temporal detail. We also Fig. 1.Model bathymetry (in metres) in the study area. BoB: Bay of Bothnia, BS: Bothnian Sea, BP: Baltic Proper, GoF: Gulf of Finland, AS: Archipelago Sea, Å: Åland Sea. Locations of the Argo profiles are marked with red circles (2012) and green squares (2013). The magenta hexagons show from north to south the locations of the Bay of Bothnia wave buoy, the Finngrundet wave buoy in the Bothnian Sea, and the Northern Baltic Proper wave buoy. The orange triangles show from north to south the locations of Kemi I and Märket weather stations. The yellow diamond shows the location of the SR5 monitoring station. The inset in the upper corner of the map shows the whole model domain and bathymetry. (For interpretation of the references to colour in thisfigure legend, the reader is referred to the web version of this article.)
A. Westerlund, L. Tuomi / Journal of Marine Systems 158 (2016) 34–44 35
2. Materials and methods
2.1. NEMO ocean model
NEMO 3D ocean model version 3.6 has been set up at FMI for the Bal-tic Sea, based on the NEMO Nordic configuration byHordoir et al.
(2013a,b, 2015). The model has been discretised on a Baltic Sea–North Sea grid with two nautical mile horizontal resolution and 56 z-coordinate vertical layers. The topmost vertical layer is 3 m, and the ver-tical resolution gradually increases with depth, being about 13 m at 150 m depth. The deepest point in the domain is in the Norwegian Trench in the North Sea at around 720 m depth and 22 m layer thick-ness. Our study area, the Bothnian Sea, is mostly less than 150 m deep.
The time step of the model was 360 s. Model domain and bathymetry are shown inFig. 1.
The NEMO Nordic configuration uses the TVD advection scheme (Leclair and Madec, 2009) for both tracers and momentum. Lateral mixing in NEMO Nordic uses Laplacian isopycnal diffusion. Isopycnal vis-cosity was set to 50 m2s−1above the depth of 30 m, and to 10−3m2s−1 below. Isopycnal diffusivity was set to 10% of viscosity. The model was configured with the LIM3 sea ice model (Vancoppenolle et al., 2009).
The model was run for the years 2012 and 2013. Temperature data from the model was saved as 3-h averages for surfacefields and daily averages for full profiles. Initial conditions for both runs were taken from the FMI operational ocean model HBM-FMI (Berg and Poulsen, 2012), which provides a daily model output for roughly the same do-main as the NEMO Nordic configuration.
2.1.1. Surfacefluxes, forcing and lateral boundary conditions
Forecasts from FMI's numerical weather prediction (NWP) system HIRLAM (HIRLAM-B, 2015) were used as atmospheric forcing. Its do-main covers the European region with a horizontal resolution of 0.15 degrees (V73 and earlier; before 6 March 2012) or 0.068 degrees (V74; after 6 March 2012). Vertically the domain is divided into 60 (V73) or 65 (V74) terrain-following hybrid levels, the lowest level being about 12 m above the sea surface. A forecast is run four times a day (00, 06, 12, and 18 UTC) using boundary conditions from the Boundary Condition Optional Project of the ECMWF. Each day of forcing was extracted from the 00 forecast cycles with the highest available temporal resolution in the model archive, varying from 1 to 6 h.
Forcing taken from HIRLAM includes the 2-metre air temperature, total cloud cover, mean sea-level pressure and 10-metre winds, and ei-ther the 2-metre dew point temperature or relative humidity, depend-ing on the availability in the model archive. Forcdepend-ing was read into the NEMO run with the CORE bulk formulae (Large and Yeager, 2004). Pre-cipitation and river runoffs were climatological (as described byHordoir et al., 2013a, 2015). The two open boundaries in the model, in the En-glish Channel and between Scotland and Norway, were configured as inHordoir et al. (2015), but with only tidal boundary surface height contribution and no storm surge model.
2.1.2. Vertical turbulence schemes in NEMO
In order tofind the optimal vertical turbulence scheme for simulat-ing the vertical temperature structures in the area, different schemes and settings available in the NEMO vertical mixing package were com-pared. The following model test runs are referenced later in this article:
•GLS/k-ε, with CB parameterisation (“k-e”)
• GLS/k-ω, with CB parameterisation (“k-w”)
• GLS/k-ε, no CB parameterisation (“k-e, CB off”)
Here GLS/k-εand GLS/k-ωrefer to the Generic Length Scale (GLS) turbulence model (Umlauf and Burchard, 2003, 2005) implemented in NEMO, which reduces tok-ε(Rodi, 1987) andk-ω(Wilcox, 1988)
2014) and thek-εmodel in particular has a long history in Baltic Sea modelling community (e.g.Svensson, 1979). The CB parameterisation refers to theCraig and Banner (1994)parameterisation of wave-breaking induced turbulence, which was used to study the effect of wave-induced mixing. We usedαCB= 100 for the wave-age related pa-rameter, as suggested byCraig and Banner (1994). For Charnock's con-stant we used the default NEMO Nordic value ofβ= 7.0⋅104.
More information of the turbulence parameterisations in NEMO can be found inReffray et al. (2015), and a review of some issues encoun-tered in numerical modelling of coastal ocean turbulence inBurchard et al. (2008).
2.2. Measurements 2.2.1. Baltic Sea Argofloats
Observations from autonomous Argofloats operating in the Baltic Sea were used to study the accuracy of the model results. Data were col-lected during two separate missions in 2012 and 2013. Thefirst mission lasted approximately six months from 17 May 2012 to 5 Dec. 2012, dur-ing which time around 200 vertical profiles of temperature and salinity were collected. The second mission lasted approximately four months from 13 Jun. 2012 to 2 Oct. 2013 with over 100 acquired profiles. Argo tracks and locations of the profiles are shown inFig. 1.
2.2.2. Moored buoys
Surface temperature measured at FMI's Directional Waveriders (loca-tions shown inFig. 1) was used to validate the model. Waveriders have a temperature sensor in the mooring eye at the bottom of the buoy and the measurement depth is c. 0.4 m below the sea surface. Temperature mea-surements are provided every 30 min, excluding the ice season, when the buoys do not perform measurements. Temperature data from FMI's wave buoy in the Bothnian Sea could not be used as the temperature sen-sor in the buoy did not measure within the specifications during the study period and was eventually replaced in 2014. Instead, data from Swedish Meteorological and Hydrological Institute's (SMHI) Finngrundet wave buoy were used for the Bothnian Sea. This buoy mea-sures SST at 0.5 m depth (http://www.smhi.se/kunskapsbanken/
oceanografi/en-vagboj-finngrundet-1.22076) and the data is provided at 1-h intervals.
2.2.3. Monitoring data
Temperature and salinity profiles from monitoring cruises to the Bothnian Sea in 2012 and 2013 were collected from the ICES (Interna-tional Council for the Exploration of the Sea) Dataset on Ocean Hydrog-raphy. A total of 55 profiles were obtained for 2012 and 82 for 2013. All these profiles originated from theR/V Aranda.
2.3. Statistical methods
The model performance was evaluated using statistical analysis based on the following equations. The mean of a dataset (observations or model) was defined as:
x¼1 N
XN
i¼1
xi ð1Þ
whereNis the number of data points andxiare the data points.
The bias of a dataset was defined as
B¼y−x ð2Þ
36 A. Westerlund, L. Tuomi / Journal of Marine Systems 158 (2016) 34–44
The root-mean-square error was defined as
and the correlation coefficient as
R¼
where the standard deviationσis defined as
σx¼
The model configuration used in the study has been previously val-idated byHordoir et al. (2013a, 2013b, 2015), who showed that gener-ally the configuration performs in a satisfactory manner. The emphasis of these previous validations has been on salinity dynamics, while validation performed by the present authors concentrated on the sea-sonal temperature variations in the model. Thek-εvertical mixing scheme with the CB parameterisation was chosen here for closer exam-ination, since it is the parameterisation used in the NEMO Nordic configuration.
This study focused on the summers 2012 and 2013, periods for which there was Argo data available for the study area. Both summers had roughly average air temperatures in the area based on long term data from the area (Pirinen et al., 2012). May 2012 was some degrees warmer compared to the 30-year average (1981–2010), while May 2013 was colder than average. This is in line with data from FMI's ice charts, according to which the ice season ended in the Baltic Sea approx-imately two weeks earlier in 2012 (15 May) than in 2013 (30 May).
There were also a number of interesting wind-induced mixing events during the period in question. Overall, these two summers proved inter-esting test cases for the model.
3.1.1. Meteorological forcing
As the performance of any 3D hydrodynamic model is highly depen-dent on the accuracy of the atmospheric forcing, the FMI-HIRLAM forc-ing was evaluated against data from coastal weather stations Märket and Kemi I (locations shown inFig. 1). Both stations can be considered generally representative of open sea conditions. The comparison showed that the forcing wind speed was fairly well represented by HIRLAM. The air temperature was forecast with good accuracy except for the beginning of June 2013. During this period there was a technical problem in the SST assimilation of the NWP system, which caused a cold bias in air temperature over the Baltic Sea lasting at least four days (Fig. 2).
3.1.2. Sea surface temperature
The statistical analysis of the modelled sea surface temperature (SST) in the northernmost sub-basins of the Baltic Sea (Table 1) shows that the SST was generally well reproduced by NEMO. Bias and RMS error are typically less than 1 °C in all the areas. The largest RMS error, of 3.31 °C with a bias of 2.95 °C, is for summer 2013 in the Bay of Bothnia. This resulted from incorrectly timed melting of ice cover in
wave buoy in that area on 11 June 2013 means that the number of avail-able data points for the summer period was slightly lower than in the other stations.
In addition to the statistics, the time series of modelled and measured SST from Finngrundet (located in the Bothnian Sea, location shown in Fig. 1) is shown inFig. 3. In both years, the model followed the seasonal cycle of SST well. In 2013 the model slightly underestimated the SST throughout the summer as was already shown in the statistical analysis inTable 1. In both years the cooling of the surface layer in autumn was slower in the model than in the measurements. The above-mentioned problem in meteorological forcing in June 2013 can be seen as a drop in modelled SST. Furthermore, the measured data shows higher tempo-ral variation than the modelled data, due to modelled SST data being Fig. 2.Two metre air temperature measured at the Märket weather station in summer 2013 compared to FMI-HIRLAM forcing used for the NEMO model. The drop in modelled air temperature in early June 2013 was due to a technical problem in the NWP system, seeSection 3.1.1.
Table 1
Modelled daily SST averages from NEMO compared against wave buoy observations in 2012 and 2013 for the study area and adjacent basins. Model RMS errors, biases (°C) and correlation coefficients R shown for spring, summer, autumn and all three seasons combined. N/A means that observation data was not available for that time, or there were too few observation data points for statistical comparison.
Baltic Proper Bothnian Sea Bay of Bothnia
RMSE Bias R RMSE Bias R RMSE Bias R
2012
Mar–May 0.60 b0.01 0.98 0.74 −0.18 0.87 N/A N/A N/A
Jun–Aug 0.70 0.23 0.98 0.78 0.33 0.98 1.07 −0.09 0.97
Sep–Nov 0.39 0.075 0.99 0.63 −0.46 0.99 0.64 −0.39 0.99 Mar–Nov 0.57 0.11 0.99 0.71 −0.08 0.99 0.89 −0.23 0.98 2013
Mar–May 0.94 0.68 0.97 N/A N/A N/A N/A N/A N/A
Jun–Aug 0.97 0.65 0.97 1.48 1.09 0.96 3.31 2.95 0.92
A. Westerlund, L. Tuomi / Journal of Marine Systems 158 (2016) 34–44 37
saved as temporal averages (as 3-h averages) and observations being in-stantaneous data.
3.1.3. Temperature and salinity from ship observations
We validated the model results against CTD monitoring data from the Bothnian Sea. The dataset was relatively sparse both in space and time. Many stations were sampled only once during the study period.
Most frequently visited monitoring stations were sampled three or four times a year.
Comparison of temperature profiles showed that the model could reproduce seasonal variations as expected in the study area. Compari-son of salinity showed that the model results were sensible for the pur-poses of this study for the Bothnian Sea. Due to its sparseness, only limited conclusions about salinity or temperature dynamics could be drawn from this data. An example of typical results is shown inFig. 4, where all salinity and temperature profiles taken at the SR5 monitoring station (location shown inFig. 1) in 2012 are shown. In the southern Baltic Sea outside our study area, the model showed some overestima-tion of salinity below the permanent halocline.
3.2. Profiles from Argofloats
Figs. 5 and 6 show measured and modelled (k-ε with CB parameterisation) temperature profiles during Argo measurement campaigns in the Bothnian Sea in 2012 and 2013. In both years, the sea-sonal thermocline developed in the spring, as expected, and the surface mixed layer reached its maximum depth and temperature in August.
NEMO was able to reproduce the vertical structure of temperature near the surface, when the atmospheric forcing was sufficiently accu-rate. The deepening of the thermocline as well as the temperature gra-dient was well represented by the model. The response of the surface layer to forcing was similar in measurements and model data. However, the warming of the surface layer in spring was slower in the modelled than in the measured profiles. Moreover, the layers below thermocline had a warm bias, and the dicothermal layer or the old winter water layer was not as pronounced as in the measurements. The most proba-ble reason behind this is the combination of initial conditions, limited vertical resolution, and over-mixing in the deeper layers. Overall, tem-perature gradients in the model were gentler than measured and
calmer periods, their structure in the model was smoother than in reality.
Argo salinity profiles from these two campaigns (not shown) did not reveal notable features with respect to the subject of this study, present-ing a picture similar to the CTD profiles inFig. 4.
3.2.1. Near-surface temperature and thermocline depth
Next, we further investigated measured and modelled temperatures near the surface along the Argofloat route and the depth of the mixed layer for the three tested parameterisations listed inSection 2.1.2. We defined near-surface temperature as the model data point that has been sampled at the depth of the topmost data point in the Argo data (typically around 4 m, depending on the profile). In our modelled data this was in most cases very close to the sea surface temperature (mean difference around 0.1 °C), but at single locations the difference was larger (up to 3 °C). Overall, when compared to the sea surface tem-perature of the model at the same points, the near-surface temtem-perature time series was very similar but slightly smoother. Thermocline depth was defined as the location of the maximum temperature gradient with depth.
Figs. 7 and 8show measured and modelled near-surface tempera-tures, along with estimated thermocline depths. Modelled near-surface temperature captured the overall seasonal heating and cooling.
However, in spring 2013 heating was slower than in the measurements with most mixing schemes, with the exception being the run with no wave-breaking parameterisation. All schemes were able to produce the near-surface temperatures in late July/August 2013.
Thermocline depths were fairly well reproduced by the model when thek-εscheme was used together with the wave-breaking
Thermocline depths were fairly well reproduced by the model when thek-εscheme was used together with the wave-breaking