Introduction

In the search for sustainable, low-carbon, power production, many coastal regions are turning to offshore wind farming to meet rising energy demands. Through continuing technological development, offshore wind farms (OWFs) are now being built, or planned, in waters that are deep enough to be seasonally stratified. This is particularly true in the coastal regions of the German Bight sector of the North Sea, where a significant fraction of the total area could eventually be covered with wind farms (Fig 1). Such large-scale development raises a number of questions regarding the impacts of these structures on the coastal ocean environment. However, we focus on the possible large-scale impact that is caused by increased turbulent mixing of the seasonal stratification as the tides drive currents back and forth through the arrays of OWF structures (see diagram in Fig 2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Map of wind farm developments in the German exclusive economic zone of the North Sea.

Red denotes wind farms that are either operational or in construction, while blue denotes areas in the planning stage. The existing Bard and Global Tech (GT) farms are labelled, as well as the North Sea Buoy 3 (NSB3) measurement station. Contours are of mean water depth in m. Data obtained from the Bundesamt für Seeschifffahrt und Hydrographie (BSH).

https://doi.org/10.1371/journal.pone.0160830.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Basic sketch of the idealised setup considered.

A typical density profile is illustrated where stratification (i.e., change in ρ(z)) is confined to a pycnocline layer with thickness, b, at height, h, from the sea bed. Only a single foundation structure is shown in the sketch.

https://doi.org/10.1371/journal.pone.0160830.g002

Stratification forms every summer in the North Sea when increased heat from solar radiation and increased air temperatures warm the upper layer of the water column, often creating temperature differences of 5 to 10°C from the cooler bottom layer. This stratification is only able to form in waters that are sufficiently deep, since both wind stress acting at the water surface, and bottom friction at the sea bed, are the primary sources of turbulence that act to destroy the stratification [1, 2]. The seasonal development of stratification, as well as its evolution in time through turbulent mixing, is known to have a large influence on the productivity and food web of the North Sea [3–5], with implications for carbon fixation [6], dissolved oxygen concentrations and hypoxia [7].

The literature on the impacts of OWFs is large, but has focussed to a great degree on marine life, e.g., birds [8, 9], sea mammals [10] and construction noise [11–13], and biofouling [14, 15]. The impacts span both negative and positive effects, such as the harmful impact of construction noise on mammalian behaviour [11], or the hypothesis that a single foundation structure could provide 2.5 times the amount of habitat that it destroys [16]. However, the hydrodynamic effects and related impacts on turbulence and mixing have not received nearly as much attention [17–19]. A closely related study was recently carried out by Rennau et al. [18] with respect to wind farm development in the Baltic Sea, where they focussed on impacts on the dense inflows to the Baltic from the North Sea. They found that a realistic OWF construction scenario, according to plans approved in 2010, could cause mixing in Baltic inflows and decrease bottom water salinity by 0.1 PSU. In an extreme scenario, with turbines filling the Danish Sounds, they calculated a 0.3 PSU decrease in bottom salinity, and a 2-metre shoaling of Baltic inflows. In that scenario wind farm generated mixing could cause the highest density inflows to disappear. It is important to note that the flow through the OWFs considered by Rennau et al. [18] is not caused by tidal currents, as we consider in this study.

The approach that we take in this study is highly idealised, and designed to understand whether or not it is possible for wind farms to generate enough mixing to significantly affect the stratification of the German sector of the North Sea. There are many uncertainties in our ability to parameterise structure-induced mixing, and the developed models are intended only as a series of order-of-magnitude estimates. Our principle goal is to bound and estimate two different time scales: (i) the mixing time scale, τ_mix, that quantifies the residence time for stratification given a mixing rate (i.e., removal) that is characteristic of the OWF foundation structures, and (ii) the advective time scale, τ_adv, that quantifies the residence time of a water parcel within a wind farm given the mean residual currents. From these two time scales we quantify the impact of the OWF structure mixing on the large-scale stratification.

The paper is organised as follows: idealised models for both turbulence production and the turbulent mixing of stratification are first developed, and then estimates of the mixing time scale are made with the use of a numerical model and different field measurements. A numerical model for residual currents is then used to estimate the advective time scale. Discussion and conclusions are presented in the final sections.

Turbulence Production

We begin by developing a parameterisation for the power consumption of the wind farm foundations. Two wind farms currently built in deep (40 m) North Sea waters (Fig 1) that are often seasonally stratified are Bard 1 and Global Tech 1 (see [5], and later discussion). Both have foundation structures that consist of a number of cylindrical sections that are mounted to the sea floor. The Bard 1 “tripile” foundation consists of three cylinders mounted vertically, however, the Global Tech 1 “tripod” foundation is more complex and described in detail in the supplementary information. These foundation structures are modelled by assuming the idealised situation of a group of cylinders in a tidally oscillating cross flow (see the diagram in Fig 2). The following formula gives the drag force that a cylinder exerts on the passing (tidal) flow, (1) where ρ₀ is the density of the fluid, C_D is the drag coefficient, A is the frontal area of the cylinder that is exposed to the free stream, with A = 3DH for the Bard tripile foundation where D is the diameter, and H the total water depth, and is the velocity of the free stream.

Since we are interested in time scales longer than the dominant tidal periods (semidiurnal periods M₂ and S₂ are dominant in the German Bight [20]), we assume that the water depth, H, is constant. In addition, despite the major role that differences in water density will play in this paper, in the above formula we shall take the water density to be constant, represented by the naught subscript on ρ₀. This is justified because the changes in water density caused by stratification, Δρ, are much less than the mean density, i.e., Δρ/ρ₀ ≪ 1, and so have little affect on the cylinder drag formula.

The value of C_D is known to depend on a number of parameters, and is one of the most important sources of uncertainty in our parameterisation of foundation structure turbulence. C_D has been found to depend on the Reynold’s number , the relative roughness κ/D, with ν the kinematic viscosity of the fluid and κ a measure of the roughness length scale, as well as on the amount of turbulence in the free-stream approach flow [21]. To retain this uncertainty in our analysis, we shall use two estimates of C_D in this study that span the expected range, C_D = 0.35, 1.0, which we refer to as the low- and high-drag cases, respectively. Note that Rennau et al. [18] use a C_D = 0.63 in their study of OWF induced mixing in the Baltic Sea. Estimates of all these parameters are listed in Table 1, with details of our estimates of OWF structure parameters outlined in S1 File in the supporting information.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Summary of representative values for wind farm and stratification parameters.

https://doi.org/10.1371/journal.pone.0160830.t001

Because the force is parallel to the velocity, and the power is given by , we can write the power removed from the flow by the cylinder as (2) However, it is more convenient to express the power consumption per unit area, so our formula will now apply within a wind farm, providing a mean value, and is given by (3) where P_str, in W m⁻², is the power removed by the structures from the flow per unit area of wind farm, with ℓ the distance between equally spaced wind turbines in a much larger farm, and represents the mean cubed current velocity averaged over a period much longer than the dominant tidal periods. This estimate of the power removed by the wind farm foundation structures is assumed to be equal to the power put into the turbulence production in order to calculate mixing rates, as described in the following section.

Using the wind farm parameters listed in Table 1, all that is needed to calculate P_str are mean cubed tidal current velocities. These were obtained from using a two-dimensional version of the TRIM-NP hydrodynamic model [22, 23]. Barotropic currents were calculated in the North Sea and part of the north eastern Atlantic Ocean using a regular grid of 12.8 km, with a finer nested grid of 1.6 km resolution in the German Bight region. The tidal currents and surface elevations are prescribed at the open lateral boundaries using the FES dataset [24], and the atmospheric forcing is provided from the regional COSMO-CLM historical simulation [25]. Simulations for the period of Jan. 1948 to Aug. 2015 are freely available (doi: 10.1594/WDCC/coastDat-2_TRIM-NP-2d). The results of the TRIM model for the calculation of were found to be within ±5% of the values using the Copernicus Marine Environmental Monitoring Service data set (http://marine.copernicus.eu) for the Northwest European Shelf Seas region for the Bard and Global Tech OWFs. Additional validation of the model is discussed in Weisse et al. [26].

This estimate of allows us to map P_str for the North Sea, as shown in Fig 3(c) and 3(d). The quantity is also plotted in Fig 3(a) and 3(b), giving an estimate of the average power that could potentially be delivered to the turbulence by the OWF foundations in any region of the German Bight per unit mass of water, with units W kg⁻¹. Two extreme cases are shown in Fig 3, and these correspond to (a, c) the low-turbulence case where we choose the low drag coefficient and the Bard foundation structures, and (b, d) the high-turbulence case where we choose the high drag coefficient and the Global Tech foundation structures. The difference between the low- and high-turbulence cases is a factor of 4.6, demonstrating that there are significant uncertainties in the amount of OWF structure induced turbulence. The maps of Fig 3 show that exhibit large regional variation within the German Bight, and generally decreases northwards with distance from the coast, which follows from the fact that the tidal velocities are strongest in these shallower regions. The large changes in shown in each of the low-turbulence panels of Fig 3(a) and 3(c) are, however, likely to be an overestimate. This is due to the dependence of C_D on the current velocity. Our low-drag value of C_D = 0.35 is chosen to be representative of circular cylinders that have surpassed the drag crisis, i.e., a sudden drop in C_D that occurs at . As the tidal currents reduce away from the coasts, we should therefore see a higher likelihood of the high-drag scenario. Despite this effect, we can generally expect lower in areas further offshore.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Map of estimated potential wind farm induced turbulence production in the German Bight region of the North Sea.

Contours and colours show in units of W kg⁻¹ (a, b) and P_str in mW m⁻² (c, d). We have shown both low- and high-turbulence cases for the Bard 1 parameters with C_D = 0.35 in (a, c), and Global Tech 1 parameters with C_D = 1.0 in (b, d). The spatial pattern in (a, b) reflects the modelled distribution of , and the panels differ only by a factor of 4.6 (the difference in wind farm parameters for the high- and low-turbulence cases). The location of both the Bard 1 and Global Tech 1 (GT) farms are shown, along with the North Sea Buoy 3 (NSB3) measurement station. Grey areas are not included because they are on average less than 10 m deep, or comprise tidal flats or estuaries.

https://doi.org/10.1371/journal.pone.0160830.g003

Values for P_str in the Global Tech and Bard 1 farms for both the high- and low-drag cases are shown in Table 1. We find an approximate range that is bounded by 3.3 and 15 mW m⁻² for the average rate of turbulence production for these two farms. These two extreme values will be used throughout the rest of the paper, and shall be referred to as the high- and low-turbulence cases.

Note that since P_str depends on the cube of the tidal velocities, it is relatively sensitive to errors in . We have compared the estimates of P_str from the TRIM model with measurements using an Acoustic Doppler Current Profiler (ADCP) at the North Sea Buoy 3 (NSB3) station, taken over a 45 day period from mid-July to the end of August in 2014 (to be discussed in the observations section). The comparison shows only a 22% difference in the value of P_str calculated from the ADCP measurements compared to the model results (with model results yielding larger values). This discrepancy is quite small considering the other uncertainties, e.g., C_D and foundation structure type, that lead to a factor of 4.6 between the high- and low-turbulence estimates.

To put the calculated turbulence levels in perspective, it is interesting to compare the turbulence production by the turbine structures, P_str, with that occurring in the boundary layer at the sea bed. The latter can be quantified by (4) where the bottom drag coefficient [27]. Taking the ratio leads to (5) and when we substitute typical values we find that despite the much larger OWF structure drag coefficient , the equivalent area of the structures relative to the total sea floor area A/ℓ² ≈ 0.5 − 1.0 × 10⁻³ is much smaller, and results in a relative power consumption of P_str/P_bot ≈ 0.037 − 0.21. This implies that the total power that is extracted from the tides by the turbine structures is approximately 4 − 20% of the bottom boundary layer extraction. As farms are built in deeper waters this percentage will increase linearly with the water depth, H. This suggests that the OWF-induced turbulence could be significant compared to the natural sources. In addition, the OWF turbulence will be distributed throughout the water column, whereas the bottom boundary layer turbulence is mainly confined to the near-bed region [27].

Mixing of Stratification: Idealised Modelling

Using the estimate of the power removal by the wind farm foundations that was developed above, we now examine the potential impact on the North Sea density stratification. The basic assumption is that the power removed from the free-stream tidal currents is converted into small-scale turbulent motions. These motions then dissipate their energy through both internal friction (i.e., viscous dissipation) and by a mixing of the water column stratification. In this section we provide a series of idealised estimates of this mixing in the North Sea by different methods under the assumption that no other turbulent mixing processes are acting simultaneously. Our main goal is to estimate characteristic residence time scales for the stratification, which can be compared with the time of seasonal stratification build up, and an advective time scale quantifying the residence time of water within an OWF.

The general approach we take to estimating the turbulent mixing of the wind farm foundations is to satisfy the conservation of energy for the turbulent flow. As a first approximation we shall neglect horizontal advection by the large-scale residual currents (i.e., averaged over periods much larger than the dominant tidal periods). The effect of residual currents is discussed in a later section. Neglecting this residual advection for the moment, we can write an evolution equation in time, t, for the mean density profile, ρ(z, t), within the wind farm as (6) where g is gravitational acceleration, and B is the vertical buoyancy flux due to turbulent mixing acting on the density field. The sign convention of B < 0 for the mixing (i.e., destruction) of the stratification is used. In turbulent mixing studies the buoyancy flux is often written explicitly as , where the primes denote values associated with turbulent eddies, and the overbar represents an ensemble average over many turbulent events. We can derive an equation for a measure of the potential energy of the water column, ϕ(t), defined by (7) where ρ_mix is the density if the entire water column were completely mixed. Note that the presence of ρ_mix only constitutes a choice of datum with respect to which the potential energy is measured, and that ϕ really represents the amount of energy required to mix the water column. Throughout the rest of this paper we refer to ϕ as the “stratification”. A mixing of the water column, caused by B < 0, corresponds to a decrease in ϕ over time. This can be seen by forming an evolution equation for ϕ, obtained directly from Eq (6) to give (8) where we have assumed that there are no fluxes of buoyancy (and potential energy) across the boundaries at z = 0, H.

In the often-used eddy diffusivity formulation of turbulent mixing, B is parameterised by the diffusive law (9) with K(z, t) the turbulent diffusivity. This parameterisation results in a diffusion equation for ρ with variable diffusivity, i.e., (10)

In order to connect the stratification and the turbulence we use a conservation equation for the kinetic energy of the turbulence that consists of a balance between three terms: (i) the production of turbulent kinetic energy, , (ii) the rate at which this energy performs work on the buoyancy field (i.e., mixing stratificaton), B, and (iii) the rate at which the energy is lost to viscous dissipation (i.e., frictional losses), ε. This three-way balance is written (11) and is often referred to as the local equilibrium hypothesis. In general, the production term may have a number of different contributions, one of which is through the turbulence generation of the turbine structures (Fig 3). In this case, we have , which is independent of time, since we are dealing with time scales longer than the dominant tidal periods, and assumed to be independent of z.

Using the basic equations above, we now describe three idealised mixing models (i.e., parameterisations) of increasing complexity in order to estimate a characteristic residence time for the stratification, τ_mix.

Model 2: A time-dependent pycnocline model

In this model we attempt to account for a time-dependent pycnocline thickness. In contrast to the last section, we choose a well-defined continuous error function (erf) representation given by (16) where h is the height of the pycnocline, Δρ the density difference across the pycnocline, and b(t) is a measure of its thickness, which we now assume is increasing in time due to mixing. An increasing pycnocline thickness with time has been observed in the laboratory experiments of Whitehead [29], where a cylindrical rod was used to produce turbulence by stirring water with a two-layer salt stratification in a tank. Whitehead [29] found that the thickening pycnocline resulted when the stratification was relatively weak (as measured by an appropriately defined Richardson number). The error function profile used above is a well known solution to the diffusion equation and also resembles the solutions found in the Whitehead [29] experiments.

Using this continuous profile, b is defined more precisely as (17) The vertical coordinate has the range 0 < z < H, with z = 0 the sea bed, and z = H the water surface. At the moment we will ignore the fact that this profile does not satisfy the appropriate no-flux boundary conditions at z = 0, H, and limit ourselves to values of b that do not exceed the depth of the pycnocline, i.e., times for which b < min{h, H − h}. We therefore cannot use this method as such to precisely calculate the time to complete mixing. By treating Δρ as a constant, but allowing b to vary in time, we can however, determine the rate of change of ϕ in time due to a thickening pycnocline as (18) Performing the derivative, the change in ϕ can be written as (19) where (20) Evaluation of this integral can be shown to result in the relation f(b) = b³/2π.

We can now use this to equate the rate of change of ϕ in Eq (19) with the amount of power that is put into mixing from Eq (13), and rearrange to get the following simple relationship for the (constant) rate of increase in interface thickness (21)

Again using the typical values from Table 1, gives a rate of interface thickening of 0.25 and 1.1 m day⁻¹ for the low- and high-turbulence cases, respectively. Note that as a quick estimate of the mixing time we can calculate how long it takes for the interface to thicken to the water depth, H/(db/dt) = 160 and 35 days. This analysis highlights how important it is to account for the thickness of the pycnocline in our evaluation of the total power that is put into turbulent mixing.

Accounting for boundaries.

We shall now extend this model to account for the appropriate no-flux boundary conditions to model the evolution of ρ(z, t) until the time of complete mixing. This model is essentially equivalent to a solution of the diffusion Eq (10), with a diffusivity, K(t), that is a function only of t.

We begin with the following initial condition of the erf-type (22) with b₀, h, H, ρ₀ and Δρ all fixed by the initial condition, and given by the values in Table 1. However, in order to satisfy no-flux boundary conditions of (23) we must add an infinite number of additional images [30], identical to our profile, but whose pycnocline lies outside the domain, i.e., we let (24) In practice, however, we have found that 6 images are sufficient to satisfy the boundary conditions and conserve mass. The evolution of ρ(z, t) proceeds while ensuring that Eq (13) is satisfied at all times.

The solution is shown in Fig 4(a) for both low- and high-turbulence cases, which show a collapse when plotted using the dimensionless time defined by t_* ≡ tP_str/(gΔρH²). This scaling of the time variable is suggested by Eq (21), and the time that it takes the interface to thicken to the full water depth, i.e., H/(db/dt). The time to complete mixing can be defined by this dimensionless time variable, and through inspection of Fig 4(b) we choose τ_mix ≡ gΔρH²/P_str, which corresponds to t_* = 1.0. This choice leads to mixing time scales of 170, and 38 days (the approximate duration of the plots shown in Fig 4a), for the low- and high-turbulence cases. Also plotted in Fig 4(b) are the constant pycnocline predictions for the fixed values of b = 6 m and b = H = 40 m (red lines).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Solution for the time-dependent pycnocline model 2.

The results of both low- and high-turbulence simulations are identical but span different time periods, since the solutions collapse if the dimensionless time t_* is used. The evolution of the dimensionless density profile 2(ρ − ρ₀)/Δρ is shown in (a), along with the decrease in the stratification over time in (b). Also shown in (b) are the constant pycnocline thickness model mixing rates given by Eq (13) with b = 6 m and b = H = 40 m for comparison (slopes of red lines). In the low- and high-turbulence cases of (a), density profiles are identical, and plotted every 18 and 4 days, respectively.

https://doi.org/10.1371/journal.pone.0160830.g004

Despite the support for a model of this type from the experiments of Whitehead [29], in which the pycnocline thickens in time due to a source of turbulence that is constant in depth, in reality the interface will eventually grow until it is influenced by the turbulence of the upper mixed layer, and the bottom boundary layer. After this time, the model should not be considered to be valid, since these sources will tend to sharpen the pycnocline and decrease b. Nonetheless, this mixing model shows that the mixing time scale is sensitive to the evolution of the pycnocline.

Model 3: A one-dimensional mixing model

In this section we shall formulate a simple one-dimensional (1D) mixing model to account for the vertical structure of turbulent mixing. The goal is, once again, to determine a representative mixing time for the turbine structures to destroy the summer stratification, assuming that no other mixing processes are occurring. We use the diffusive model for the density fluxes from Eq (10) to determine the density profile over time, where K = K(z, t) now. In addition, we use the local equilibrium form for the turbulent kinetic energy as in Eq (11) with the production term , and the buoyancy flux term B = −KN², i.e., (25) Once again is taken independent of t, and is assumed independent of z. It remains to provide parameterisations for both K and ε.

Following successful 1D turbulence modelling in the Irish Sea by Simpson et al. [27], we shall begin by using a similar approach. The dissipation of turbulent kinetic energy, ε (in m² s⁻³, or equivalently W kg⁻¹), shall be modelled by the relation (26) where C₁ is an empirical constant (see Table 2 for a summary of constants), l is the turbulent length scale which must be modelled, and we have chosen to use the turbulent velocity scale instead of the turbulent kinetic energy, k. In addition, the diffusivity shall be given by (27) where C₂ is another empirical constant.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Empirical constants used in the one-dimensional mixing model.

https://doi.org/10.1371/journal.pone.0160830.t002

In order to choose an appropriate form for the turbulent length scale, l(z), we account for the effects of both upper and lower boundaries, and density stratification. Following the study by Simpson et al. [27], the boundary effects were accounted for by choosing l = l_H(z) ≡ κz(1 − z/H), where κ = 0.41 is von Karman’s constant. The influence of stratification will be modelled by the methods presented in Umlauf [31] and Jackson et al. [33], which can be applied to shear free stratified mixing situations. In the limit of vanishing shear, and finite stratification, Umlauf [31] has let l → l_N ≡ C₃ k^1/2 N⁻¹, with C₃ an empirical constant (Table 2). This form can be understood in terms of the Ozmidov length scale, defined as L_O ≡ (ε/N³)^1/2, which is proportional to the largest vertical eddy size in turbulent stratified waters, at large Reynolds numbers. If we associate l_N ∼ L_O and take ε ∼ k^3/2/l in the above definition, we get the scaling l_N ∼ q/N. Therefore, we model the turbulent length scale in bounded, stratified waters by the form (28)

Our 1D mixing model consists of numerically solving Eq (25) to give q(z), once the forms in Eqs (26)–(28) are substituted. This is then used to calculate K(z), so that the diffusion equation can be stepped forward in time to give a new ρ(z, t). We use the Newton-Raphson method to solve the energy balance Eq (25), and a semi-implicit finite difference predictor-corrector scheme for the diffusion equation.

The results of the modelling are shown in Fig 5, where a number of important features can be seen. First, the evolution of the density profile shows that the pycnocline reaches an approximate equilibrium thickness close to the initial thickness of b = 6 m. We would, therefore, expect that the stratification, ϕ, should decrease at near constant rate as described in Eq (13). This is indeed what is seen in Fig 5(b), and reflected in the mixing efficiency, 〈R_f〉, in Fig 5(c) where after a short initial adjustment it remains at a relatively constant value up to the time of complete mixing at the end of the simulation. Note that as with the previous case 2 mixing model, there is generally a collapse of the data for the low- and high-turbulence simulations if we use the dimensionless time, t_*. However, slight differences in 〈R_f〉 were observed between these cases (Fig 5c).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Results of the one-dimensional mixing model 3.

(a) The evolution of the dimensionless density profile, 2(ρ − ρ₀)/Δρ, as a complete mixing of the stratification occurs. The results of the low- and high-turbulence simulations show very similar results when the dimensionless time, t_* is used. The time to complete mixing is found to be 469 and 107 days for the low- and high-turbulence cases. (b) The evolution of ϕ with t_* where the red lines indicate the equivalent constant mixing efficiency lines at fixed pycnocline thickness (b = 6 and 9 m are shown). (c) Bulk mixing efficiencies 〈R_f〉, and the corresponding effective pycnocline thicknesses, b. Slight differences in mixing efficiency for the low- and high-turbulence cases are seen.

https://doi.org/10.1371/journal.pone.0160830.g005

This behaviour of the pycnocline to preserve its thickness can be seen to be an asymptotic feature of our model with a strongly stratified water column and a constant-in-depth production of turbulence. If the boundary influence on the turbulent length scale is negligible compared to the stratification then l_H ≫ l_N so that l → l_N. The energy balance leads to the relationship that q ∝ N^−1/2. This relation then gives a diffusivity K ∝ N⁻², which is just the right vertical dependence so that the pycnocline does not grow in time. This can be understood as the diffusive spreading of the pycnocline being just balanced by the divergence of K. In other words, when the pycnocline is strong, and located away from the boundaries, the mixing model predicts that the pycnocline thickness remains relatively constant in time. This lends support for the case 1 mixing model with a constant thickness.

However, despite the relatively constant pycnocline thickness, we find that there is an enhanced rate of mixing compared to the case 1 model with b fixed at the initial condition. Assuming that the mixing efficiency, R_f = 0.17, an equivalent value of b required to produce the simulated decrease in ϕ is found to be approximately b = 9 m (Fig 5b, red lines).

Using the above results, we can calculate the appropriate stratification residence time as the time required for complete mixing. This gives a residence time of τ_mix = 469 and 107 days for the low- and high-turbulence cases, respectively. A summary table of our mixing time scale for all three idealised modelling cases is shown in Table 3.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Summary of the mixing time scale estimates, τ_mix.

https://doi.org/10.1371/journal.pone.0160830.t003

Model Predictions in the Context of Observations

In this section we use a number of available observations and recent measurements in the German Bight region of the North Sea to make inferences as to the potential impact of the wind farm foundation structures on turbulence and mixing. The data sources used are through both measurements made using ocean gliders and fixed moorings, as well as numerical model results.

Description of data sources

Measurements made from a number of different instruments, and at different times are each described here in turn.

Gliders.

Field measurements were collected using autonomous ocean gliders [34, 35] during the summers of 2012 and 2014. Gliders adjust their buoyancy by means of an internal pump in order to profile the water column vertically from near the surface to approximately 2 − 3 m above the bed, while using wings to provide a forward momentum. Each glider is equipped with a Seabird CTD to provide measurements of conductivity, temperature and pressure, at a frequency of 0.5 Hz that are used to calculate stratification parameters. The surfacing times of the gliders were approximately 3 hours apart, and allow for the determination of an exact position. These positions for the 2012 and 2014 campaigns are shown in Fig 6.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Glider positions for the 2012 and 2014 campaigns.

https://doi.org/10.1371/journal.pone.0160830.g006

Mooring.

A series of temperature and salinity measurements made at the NSB3 station by the BSH (Bundesamt für Seeschifffahrt und Hydrographie) were used to calculate stratification related parameters during the summers of 2004, 2005, 2009-2013. Temperature and conductivity were recorded every hour at the fixed depths of 4, 6, 10, 15, 20, 25, 30, 35 m for temperature, and 6 and 35 m for conductivity. Data was recorded in 40 m total water depth, at a sampling interval of 10 mins, and averaged to produce hourly values. There are significant gaps in the data for many of the sensor depths, and it was only possible to calculate ϕ(t) for the year of 2009. However, changes in temperature, salinity, and density across the pycnocline were calculated for the other years when available.

ADCP.

During the summer of 2014, an Acoustic Doppler Current Profiler (ADCP) was also moored on the sea bed close to the NSB3 station from July 16 to August 31. The ADCP was a RDI 600 kHz sampling every 10 minutes using an ensemble of 32 pings. This allowed for the measurement of both components of the horizontal velocity over the depth range of 5 − 38 m in a total water depth of 40 m.

In addition to these measurements, we also use results from the numerical simulation of North Sea stratification in the German Bight.

Stratification modelling.

Data from simulations carried out and discussed in van Leeuwen et al. [5] will be used to determine the duration of seasonal stratification, and the variability that is present in the German Bight region. These simulations were carried out for a 51 year time period within the North Sea domain (51 to 60°N and −5 to 9°E) using the General Estuarine Transport Model (GETM, see www.getm.eu and [36]). We shall only discuss results of the duration of seasonal stratification at the location 54.5°N, 7.0°E near to the NSB3 measurement station. For full details of the simulations, and validation see van Leeuwen et al. [5].

Observed stratification

As an example, we begin with an overview of the stratification for July and August 2014 that can be seen in the depth-time plot of water temperature, T, shown in Fig 7(a). At the beginning of the measurements the water column has developed a strong stratification with a well-defined pycnocline centred at around 12 m depth. With the arrival of storm Bertha on Aug. 9, persistent winds in excess of 10 m s⁻¹ eventually produced a complete mixing of the water column, and an end to the seasonal stratification. This build up of stratification and eventual destruction by fall and late-summer storms is typical of the North Sea seasonal cycle, and will be discussed further below.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Glider measurements of North Sea stratification, summer 2014.

The upper panel shows the conservative temperature in depth and time. The pycnocline region is indicated by the light coloured lines, corresponding to the b₁ definition, and the depth of the sea bed is indicated with dark grey fill. The bottom panel shows thermocline thickness using two different methods (b₁, b₂), and total temperature difference, ΔT. The approximate times of Storm Bertha are indicated by the dark bar at the top of the panel, and are defined by times of sustained wind speeds in excess of 10 m s⁻¹.

https://doi.org/10.1371/journal.pone.0160830.g007

A critical factor in determining the bulk mixing efficiency, 〈R_f〉, was the thermocline thickness, b. This is a quantity that can be readily calculated from the glider CTD measurements. However, exactly how b should be defined is not clear, and we are not aware of a generally accepted definition. We shall therefore, use two different definitions with the first, b₁, corresponding to the depth interval over which 80% of the total change in temperature occurs. In the calculation of b₁ we use a sorted temperature profile, T_*(z), obtained by rearranging the measured T values in the profile to be monotonically increasing with height. This is done to avoid unrealistic values of b₁ due to the presence of overturns. The second definition, b₂, is based on an integral representation that has previously been used in the study of stratified shear layers [37] given by (29) where we define the temperature difference across the thermocline as ΔT ≡ T_u − T_l with T_u and T_l the temperatures in the upper and lower layers, respectively. Note that the definition of b₂ has the property that b₂ = ΔT/(dT/dz)_max, in accord with Eq (17), when T_*(z) has the hyperbolic tangent form.

Time series of thermocline thicknesses b₁, b₂, as well as ΔT, measured during the campaign by the glider, are shown in Fig 7(b). Average values of b₁, b₂ over the stratified period, before the arrival of storm Bertha, correspond to 6.0 m and 5.8 m, respectively, and therefore a typical value of b = 6 m has been listed in Table 1. Note that if we had defined b₁ to be 90% of the total change in T then we find an average value of 8.2 m. We have chosen to use the 80% definition, since it corresponds very closely with the independent definition of b₂. However, this somewhat arbitrary definition of b is a limitation of the idealised modelling. It is also not clear what the spatial variability of b is in the German Bight, or how it depends on H, or other factors.

Stratification and the seasonal cycle

To understand the relevant time scales over which the enhanced OWF-induced turbulence can act to mix the stratification, we must first quantify the typical duration of the seasonal stratification (a second, advective time scale is estimated in the following section). Here we appeal to the North Sea-wide modelling study of van Leeuwen et al. [5], for a representative point close to the NSB3 measurement station. A histogram is plotted in Fig 8 that shows the frequency of the duration of longest continuous stratification for simulations carried out over the years 1958-2008. The figure shows a significant amount of inter-annual variability, with a mean duration of seasonal stratification of 85 days and a standard deviation of 34 days. It will be shown in the following section that this seasonal stratification time scale is not the limiting scale for the mixing of stratification, since advection of stratification through the German Bight is more rapid.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Histogram of the period of longest continuous stratification from the model of van Leeuwen et al. [5].

The data are presented for the location 54.5°N, 7.0°E near to the NSB3 measurement station. Mean and median values for the time of seasonal stratification are 85 and 80 days, respectively. Data supplied by S.M. van Leeuwen, Cefas (UK); contains public sector information licensed under the Open Government Licence v3.0, UK Crown Copyright.

https://doi.org/10.1371/journal.pone.0160830.g008

We now consider all dominant processes that are acting to build and destroy the stratification over the course of the seasonal stratification cycle. We begin with a simple depth-integrated model similar to that presented in Simpson and Sharples [1], where an energy balance for the potential energy of the North Sea water column may be written in the general form of (30) The terms on the right hand side (with all ) represent the power input due to atmospheric heating and net incoming solar radiation, and the possible turbulent mixing processes arising from turbine structures, bottom friction, and surface wind stress (including wind waves), respectively, with the asterisk indicating that each must be multiplied by an efficiency.

We evaluate the relative importance of the P_str term by plotting the observed ϕ(t) from measurements collected near NSB3. In order to calculate ϕ(t), we use Eq (7) together with the following form for the water density (31) where Δρ ≡ ρ_l − ρ_u, and the function f(z) ≡ [T(z) − T_u]/(T_l − T_u). This form assumes that the salinity profile has the same vertical structure as T(z). We make this assumption because the variation of salinity in the pycnocline is either not known (mooring data, with measurements at only 6 and 35 m), or unrealistic densities are calculated within the pycnocline from the mismatch and time response of the conductivity and temperature sensors (glider data). However, since in each case the upper and lower mixed layer salinity is accurately measured, this contribution to density is accounted for through the Δρ factor.

At these times, and in the measurement locations, we expect very little influence from OWFs that could produce a P_str contribution. We can therefore estimate the relative size of P_str by overlaying our previously determined estimate based on a fixed width thermocline. In other words, in Fig 9 we can compare the rate of growth of stratification (dϕ/dt) with the rate of removal by the turbine structures, i.e., (32) The left equation is given by the observations and the 〈R_f〉P_str term represented by the various lines depending on the b, C_D, and foundation type chosen, i.e., the different bulk mixing efficiencies 〈R_f〉 = R_f b/H and low- and high-turbulence cases. Fig 9 suggests that the mixing due to wind farm turbulence could represent a significant contribution to the overall potential energy budget of the stratification.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Measurements showing the build up of stratification, ϕ, over the summer months (data points) versus the rate of stratification removal by the turbine foundation structures (straight lines).

Measured stratification data are from a thermistor mooring at NSB3 in summer 2009 (black points), and glider data collected close to NSB3 from summer 2014 (green points), and from larger scale transects passing through NSB3 in summer 2012 (blue points). None of these stratifications are expected to be affected by OWFs. The lines plot the mixing power of the turbine foundations at two different thermocline thicknesses, b = 6, 9 m, based on results from models 1 and 3 using the slope R_f P_str b/H, or equivalently the slope 〈R_f〉P_str. Also, both the low- and high-turbulence cases are shown with the dashed and solid lines, respectively. All lines are shown starting from an arbitrary initial date of May 1. Mooring data from 2009 were obtained from BSH.

https://doi.org/10.1371/journal.pone.0160830.g009

As shown in Fig 9, ϕ gradually increases in the early summer months as the P_heat term increases relative to the other stratification removing terms. However, there is significant variability of the seasonal peak in stratification, as seen in the glider measurements from 2012. Based on the measurements in 2009, and 2014, we can take a typical value for the peak stratification as ϕ_max = 5.0 kJ m⁻²(Table 1). In order to provide a greater justification for taking this representative value based on only two seasons, we analyse other years from the mooring data at NSB3 (Fig 10). The data from these years have captured the total change in temperature and salinity (and therefore also density) across the pycnocline, but do not have the vertical resolution required for computing ϕ. Assuming the error function form for the density profile, as described in mixing model 2, as well as the other typical values for b and h listed in Table 1, this representative value of ϕ_max = 5.0 kJ m⁻² corresponds to a density difference of Δρ = 3.1 kg m⁻³. This value of Δρ is also found to be the average peak measured from the data in years 2009-2012 (Fig 10), supporting our choice of ϕ_max = 5.0 kJ m⁻². The measurements from 2004 and 2005, plotted in Fig 10, unfortunately do not resolve the peak in summer stratification, and have therefore not been used to calculate the average peak Δρ.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Measurements showing the difference in (a) temperature, (b) salinity, and (c) density across the pycnocline over the summer months at NSB3.

These differences were calculated between measurements in the surface layer at 4 or 6 m deep (depending on the availability of data) and at depths of 35 m in the lower layer. Data were obtained from BSH.

https://doi.org/10.1371/journal.pone.0160830.g010

It should be noted that the salinity contribution to Δρ is often large and highly variable (Fig 10b). This is expected to be the influence of the freshwater plume from the Elbe River that flows along the eastern coastline of Germany and Denmark. This freshwater is expected to have less of an influence in the German Bight further away from the coasts, and may lead to lower values of Δρ and ϕ there. This importance of freshwater influence was also noted in van Leeuwen et al. [5], where they found that conditions near NSB3 were classified as a region of freshwater influence 12% of the time.

Accounting for Advection and Wind Farm Size

Thus far our model analysis has assumed that the wind farms occupy the entire sea, so that there is effectively no advection of stratification through the wind farm area. However, it is known that there is a residual mean flow, superimposed on the tidal motions, that is responsible for an advection of the water column through the North Sea [38]. When the wind farm occupies a finite extent then the mixed water will eventually be replenished with fresh stratification from outside, and the time scale for this replenishment is denoted by τ_adv. This effect was ignored in the previous sections in order to estimate a mixing time scale, and time scale for the seasonal stratification, which are both independent of the wind farm size.

What is relevant for our analysis is the amount of time a water column is expected to spend within the elevated mixing region of the wind farm. This advective time scale is determined by both the length of the farm, and the strength of the mean circulation. Two scenarios will be examined: (i) OWFs on the scale of current construction (e.g., the Bard and Global Tech farms) with a typical length of L = 8 km, and (ii) a future scenario where OWFs fill the study area shown in Fig 11(a) with a length scale of L ≈ 100 km. These two scenarios will allow for impacts that span the full range of OWF development options in the German Bight.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Residual currents.

(a) Map of residual (depth-averaged) currents in the German Bight with contours indicating the magnitude, and arrows indicating direction and magnitude. The averaging period is taken over all years between 1958-2014. The dark outlined area represents the study region referred to in the text. (b) Histograms showing the distribution of times (τ_adv) that correspond to the minimum residence time within the study area of 70, 50, and 30% of the drifters. The mean values of each distribution are indicated by the lines at the top, and bin edges are at τ_adv = −5, 5, 15, 25, ….

https://doi.org/10.1371/journal.pone.0160830.g011

Residual currents are estimated from the TRIM model [22], as previously described in the tidal current estimates section. To get a picture of the mean (depth averaged) circulation through the German Bight region, we have used an averaging period running over the course of the entire simulation length, from years 1958 to Aug. 2015 (freely available from doi: 10.1594/WDCC/TRIM-NP-2d-PCA_ResCurr). The result is shown in Fig 11(a), and shows that a bulk advection that is towards the north-east through the proposed OWF development area.

To estimate a representative τ_adv for the relatively small-scale (L = 8 km) OWFs currently built, we make the assumption that the spatial gradients of the residual current do not vary appreciably over the length of the OWF. By inspection of Fig 11(a), this appears to be a good assumption, and allows us to write , where refers to the time-averaged residual current magnitude. Using the representative range of mean residual current speeds in the study region of cm s⁻¹, and L = 8 km, gives 1.5 < τ_adv < 9.3 days as a typical range in the German Bight. Specifically, for the Bard and Global Tech OWFs we find cm s⁻¹, and therefore, τ_adv = 4.0, 4.4 days, respectively.

When considering mean advection speeds through larger regions where significant spatial gradients exist, it is not possible to use a simple spatial average for the calculation of mean residence times. This is due to the Lagrangian nature of particle advection, since particles spend more time in regions of lower velocities than in regions of larger velocities. This makes the computation of τ_adv more complicated for the large-scale OWF development scenario. The procedure that we take to estimate τ_adv is as follows. First, the study region indicated in Fig 11(a) is seeded with 2,000 evenly distributed Lagrangian drifters. The position of each drifter is integrated over a time period of 90 days, and the total duration that the drifter spends in the study region is calculated. This is done with an integration time beginning on each separate day from May 1 until July 1, sufficient to cover the summer stratification season when the 90 day integration time is included, for the years 2001-2005, 2009, 2012, 2014 (largely chosen to match with years in which stratification data is available). Each 90 day integration is treated as an independent estimate of the distribution of time values from each drifter, which is the duration that the drifter spends in the study area, allowing for exit and reentry. The minimum times recorded by 70, 50, and 30% of the drifters for each distribution are then calculated, giving three different estimates of τ_adv, and their subsequent distribution over all experiments is shown in Fig 11(b). In other words, for a single integration experiment with 2,000 drifters, an estimate of τ_adv represents the minimum time spent in the study area of 70% of the drifters in the case of the blue histogram. By calculating means of these τ_adv distributions we arrive at three estimates of τ_adv = 28, 41, 52 days, that correspond to the averaged minimum residence time within the study area of 70, 50, and 30% of the drifters, respectively. Note that if any one of the years that drifter simulations were performed are excluded from the analysis, the mean values of τ_adv change by no more than ±3 days.

This analysis shows that the advective time scale provides the greatest control over the total amount of mixing, since it is generally much shorter than the time scale of the seasonal stratification (∼85 days), and also the majority of the mixing time scales presented in Table 3. However, there is a very large variability that is present in these distributions, and it can be seen that the difference in the means of the 70, 50, and 30% values of τ_adv are smaller than the variability of a single distribution. In the next section we provide estimates of the reduction in stratification that these time scales imply.

Discussion

Impacts on the large-scale stratification

As the main results of this study we have estimated both the mixing and advective time scales, τ_mix and τ_adv, and compared them to the observed rate of build up of stratification, dϕ/dt. It is interesting now to relate these two time scales by looking at the case of a wind farm of finite size in a steady residual current with all other stratification removing and building processes neglected, i.e., assuming a balance between advection and wind farm mixing. In formulating such a model, we make the assumption that a representative advection speed through the farm is given by L/τ_adv. Then, if the stratification entering the wind farm is ϕ₀ = ϕ(x = 0), the total change in stratification after passing through the farm, Δϕ ≡ ϕ₀ − ϕ(x = L), is given by 〈R_f〉P_str τ_adv. In this case, the fractional reduction of stratification can be represented by the ratio of our two fundamental time scales, i.e., (33) where here τ_mix ≡ ϕ₀/〈R_f〉P_str is the time scale to mix ϕ₀. This shows explicitly that the impact of the OWF structures on the stratification is dependent on the ratio of the advective and mixing time scales. Note that we could have also formulated this model in a Lagrangian perspective in terms of the time spent by a drifter within the elevated mixing region of the OWF; the result is the same. This idealised model can be used as a very rough guide to judge the spatial extent to which OWFs can be built before significant impacts to the stratification are expected.

Table 4 shows a summary of the reduction in stratification that will occur through OWFs from our two scenarios of (i) existing farms where L ≈ 8 km, and (ii) for the study region of Fig 11(a) with L ≈ 100 km. This is done also for the three different mixing models, for both low- and high-turbulence conditions, and for the three different drifter coverages of 70, 50, and 30%. In computing τ_mix for the large farm case in Table 4, we have used a P_str that is obtained from a spatial average over the study area, rather than the values typical of the Bard and Global Tech farms listed in Table 1. These values were found to be 2.8 and 13 mW m⁻² for the low- and high-turbulence cases, close to that found for the Bard and Global Tech farms. The difference in the reduction of stratification from the 70 and 30% drifter coverages is approximately a factor of two, and the inter-annual variability of τ_adv can be up to a factor of 4 or more. This shows that the impact of the OWF-induced mixing is sensitive to the large variability of the advection by residual currents through the German Bight, which is largely driven by the variability in wind forcing.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Estimates of the reduction in stratification Δϕ/ϕ₀.

https://doi.org/10.1371/journal.pone.0160830.t004

Note that the computed reduction in stratification assumes that no stratification building processes are acting. If we were to include these then we could modify Eq (33) by adding the term −τ_adv/τ_s to the right hand side (with a negative sign representing the addition of stratification), where τ_s is a typical time scale for the stratification to build up to the level ϕ₀, i.e., the seasonal stratification time scale. In addition, we have made the assumption that there are no OWFs upstream so that the water column entering the elevated mixing region of the OWFs is “pristine”. The results indicate that the high-turbulence scenario with extensive OWF development in the German Bight (L = 100 km) could significantly impact the large-scale stratification, whereas current construction levels have only a very small impact.

Conclusions

With the large-scale planning and construction of offshore wind farms (OWFs) in coastal seas that are deep enough to support the development of seasonal stratification, we have provided a first, order-of-magnitude estimate of the potential impact this development may have on the mixing of stratification. This mixing is induced by the turbulent wake of the OWF foundation structures as the tidal currents continuously move past. Using idealised modelling we have developed a series of estimates of the mixing time scale that characterises the time to complete mixing. This time scale was generally found to be larger, though comparable to, the summer stratification period, showing that the estimated mixing could be important for the development of stratification. However, these estimates were found to be sensitive to both the time evolution of the pycnocline thickness and the drag of the foundation structures—both of which are uncertain. In addition, estimates of the advective time scale, the duration a water parcel is likely to spend inside the enhanced mixing region of a wind farm of finite size, show that, for a significant impact on the stratification, extensive regions of the North Sea must be covered in OWFs. At the current construction levels (Fig 1) we should not expect any large-scale changes to the stratification of the North Sea. However, in future development scenarios where OWFs fill large portions of the German Bight (Figs 1 and 11a) we could expect significant reductions in the stratification. These impacts are expected to be highly variable due to the dependence of the advective time scale on wind forcing and its variability.

Given that a significant reduction in stratification could be possible from large-scale OWF development, it is still uncertain whether this impact could have positive or negative effects (see introduction and [10]). However, one possible mitigation strategy could involve the use of floating or semi-submerged OWF platforms [44] in order to minimise the effects.

Future work on the impact of OWFs on the mixing of stratification could focus on a number of different directions. Firstly, it is necessary to better understand the local turbulence production and induced mixing of the different OWF foundation structures so that more accurate mixing parameterisations can be developed. Second, these parameterisations could then be used in larger scale regional models to conduct ensemble averages over many seasons, thus providing a statistical measure of the spatial changes for given OWF development scenarios. Finally, it is also important to understand the effects of the enhanced mixing on the scale of an individual farm as they currently exist, i.e., approximately 10—20 km. Over this scale the enhanced mixing could show cascading effects on nutrient levels, ecosystems, and marine mammals. This is due to the strong interaction between turbulence levels and the growth of phytoplankton in general [45], as well as in the North Sea in particular [46]. The physical-biological interactions in OWFs are particularly unknown [19], and it is important to quantify these effects in order to have a complete understanding of the true impact of offshore wind farms.

Supporting Information

Acknowledgments

We would like to thank Dr. S.M. van Leeuwen (Cefas, UK) for her helpful review of the manuscript, and for providing the simulated stratification data, made avialable through the Open Government Licence v3.0, UK Crown Copyright. We also thank R. Kopetzky, M. Heineke, B. Peters, and the captains and crew of the Ludwig Prandtl for their help in the collection of the glider data. The Bundesamt für Seeschifffahrt und Hydrographie (BSH) and the Deutsches Ozeanographisches Datenzentrum have made the mooring data available to us. We are also grateful for the help from BSH in the deployment of the ADCP in 2014.