Emergence of Large-Scale Hydrodynamic Structures Due to Atmospheric Offshore Wind Farm Wakes

The potential impact of offshore wind farms through decreasing sea surface wind speed on the shear forcing and its consequences for the ocean dynamics are investigated. Based on the unstructured-grid model SCHISM, we present a new cross-scale hydrodynamic model setup for the southern North Sea, which enables high-resolution analysis of offshore wind farms in the marine environment. We introduce an observational-based empirical approach to parameterize the atmospheric wakes in a hydrodynamic model and simulate the seasonal cycle of the summer stratification in consideration of the recent state of wind farm development in the southern North Sea. The simulations show the emergence of large-scale attenuation in the wind forcing and associated alterations in the local hydro- and thermodynamics. The wake effects lead to unanticipated spatial variability in the mean horizontal currents and to the formation of large-scale dipoles in the sea surface elevation. Induced changes in the vertical and lateral flow are sufficiently strong to influence the residual currents and entail alterations of the temperature and salinity distribution in areas of wind farm operation. Ultimately, the dipole-related processes affect the stratification development in the southern North Sea and indicate potential impact on marine ecosystem processes. In the German Bight, in particular, we observe large-scale structural change in stratification strength, which eventually enhances the stratification during the decline of the summer stratification toward autumn.


INTRODUCTION
With sustainable energy generation through offshore wind farms, kinetic energy is withdrawn from the atmosphere and consequently horizontal momentum reduces on the leeward side of the respective wind turbines. The consequences of this energy extraction are atmospheric wakes, which are characterized by the downstream reduction of the mean wind speed and the development of increased turbulence along the wind speed deficit (Lissaman, 1979;Fitch et al., 2012Fitch et al., , 2013Volker et al., 2015;Akhtar et al., 2021). In large wind turbine clusters, the individual wakes merge downstream into a single wind farm related wake structure, whereby wind direction and wind turbine layout play an important role (Frandsen, 1992;Li and Lehner, 2013). By using in-situ airborne measurements and satellite Synthetic Aperture Radar (SAR) data, recent studies observed atmospheric wakes behind offshore wind farms in the North Sea and derived information about spatial extension and intensity of wind farm wakes Hasager, 2005, 2006;Li and Lehner, 2013;Emeis et al., 2016;Djath et al., 2018;Platis et al., 2018Platis et al., , 2020Platis et al., , 2021Siedersleben et al., 2018;Djath and Schulz-Stellenfleth, 2019;Cañadillas et al., 2020). The measurements demonstrated a strong correlation between atmospheric stability and the dimension of wakes. Wind speed deficits propagate much further under stable conditions and were detected up to 70 km downstream of the respective wind farms (Djath et al., 2018;Cañadillas et al., 2020). At the same time, the superposition of wind farm wakes is a decisive factor for the wake dimension (Djath et al., 2018). The magnitudes of wind speed deficits depend for the most part on the mean wind speeds at wind farms and the wind farm drag by the wind turbines.
The impact of energy extraction by wind turbines cascades down to the sea surface and reduces the wind stress at the sea surface boundary. Although the strongest wind speed reductions emerge at hub height [∼20-30%, Cañadillas et al. (2020) and Platis et al. (2020)], the SAR measurements showed that the wind speed at 10 m height is still reduced by about 10% (Christiansen and Hasager, 2005;Djath et al., 2018). Since, not only residual currents but also mixing of the surface mixed layer are primarily shear-driven (Kantha and Clayson, 2015), anomalies in the wind field can have severe consequences for the upper ocean dynamics. This applies particularly to the shallow North Sea, where the general circulation is significantly controlled by the atmosphere (Sündermann and Pohlmann, 2011). In addition, thermodynamic processes in marginal seas, like the North Sea, are strongly sensitive to variations in atmospheric forcing (Schrum and Backhaus, 1999;Skogen et al., 2011;Schrum, 2017). Turbulent processes near the sea surface boundary determine vertical fluxes (e.g., heat, water, and momentum) between atmosphere and ocean (Kantha and Clayson, 2015). Thus, a reduction of shear-driven turbulent mixing can result in changes of the heat content in the upper ocean and associated surface heating or cooling.
Earlier studies started to investigate the oceanic response to wind stress anomalies in idealized model approaches (Broström, 2008;Paskyabi and Fer, 2012;Paskyabi, 2015). The studies demonstrated the occurrence of an elevation pattern and corresponding up-and downwelling at the sea surface, which results from the adjustment of wind-driven Ekman transport along the wake-impacted area (Ludewig, 2015). For the case of an idealized stratified water column, Ludewig (2015) showed that the convergence and divergence of surface water masses trigger dipole-related vertical transport and results in perturbations of the pycnocline. At this, Ludewig (2015) found associated vertical velocities of several meters per day in a German Bight scenario, which induced anomalies in the temperature and salinity distribution.
The latest results of recent studies on potential thermoand hydrodynamic effects by wind farm wakes raise concerns about substantial changes to the hydrodynamics of the North Sea. In particular, since the offshore development in the North Sea is growing continuously, questions about environmental consequences of offshore wind farming become crucial for prospective research. Broström (2008) already discussed the impact of the wake-induced dipoles on the marine temperature field and related regional nutrient availability and stated the need for more realistic investigations of the wake-induced forcing mechanisms.
In this manuscript, we address the potential mean impact of atmospheric wakes on the hydrodynamic system of the southern North Sea in consideration of the current state of fully commissioned offshore wind farms (Figure 1). For this purpose, we developed a new cross-scale hydrodynamic model setup, which allows high model resolution in areas of offshore wind farm production. In order to investigate the effect of wake-associated wind speed deficits, we introduce a top-down wake formulation build upon former wake models and SAR observational data. By implementing the simplified wake parameterization into the hydrodynamic model, we aim to expand the results of existing studies and complement missing knowledge about large-scale and long-lasting effects in a realistic North Sea scenario. In particular, we answer questions about the magnitude and the dimension of hydrodynamic changes and provide a first overview about the spatial perturbations due to offshore wind farm wakes.

ATMOSPHERIC WAKE PARAMETERIZATION
Over the years, several empirical and analytical models have been developed to parameterize the downstream wake effects FIGURE 1 | Bathymetry of the model domain (data downloaded from https://www.emodnet-bathymetry.eu/, June 2019) and the fully commissioned offshore wind farms in the southern North Sea, which were considered in this study (red polygons; status as of November 3, 2020, data provided by https://www.4coffshore.com/windfarms/). in the atmosphere (Jensen, 1983;Frandsen, 1992;Emeis and Frandsen, 1993;Frandsen et al., 2006;Emeis, 2010;Peña and Rathmann, 2014;Shapiro et al., 2019). While trying to describe wake deficit and wake length, most formulations focused on the relative difference between wind speed along the lee side of wind farms and the undisturbed wind flow. For instance, the analytical model proposed by Emeis (2010) considers wind turbines as additional roughness and is based on the equilibrium between momentum extraction by wind turbines and the replenishment of momentum through turbulent fluxes from above. At this, the change in roughness is assumed proportional to the wind turbine drag coefficient and the wind speed at hub height. However, existing wake models, with regard to wind speed reduction, have so far only been developed for the processes at hub height. Due to the vertical changes of the wake pattern (Frandsen et al., 2006;Emeis, 2010;Fitch et al., 2012;Akhtar et al., 2021), the former models can thus not directly be used for processes near the sea surface boundary. For this reason, we adjusted former model assumptions based on SAR measurements so that they become applicable for numerical ocean modeling.
In a strongly simplified first order approximation, which neglects dependencies arising from specific wind farm characteristics, we parameterized the wind speed deficits resulting from operating wind farms and reduced the mean wind speed in dependence of the respective wind direction. Other impacts, such as turbulence changes or effects on local weather conditions remain unconsidered.
The parameterization for the downstream wind speed reduction [Eq. (1)] is based on earlier studies (Frandsen, 1992;Frandsen et al., 2006) and alters the undisturbed wind field u 0 by the wind speed deficit u. The parameterization is applied in a reference coordinate system, where x is the downstream distance aligned to the prevailing wind direction and y defines the distance from the central wake axis. Accordingly, the deficit u(x,y) consist of two components describing the downstream wake recovery and the width of the wake structure.

Wake Recovery
The formulation of the downstream velocity deficit is based on the concept of the model by Emeis (2010) (E10 model). The model is a top-down approach model, meaning that the wind farms are considered as one unit of additional roughness, and describes the wake recovery by an exponential decay function. The E10 model was validated in recent studies, which showed that the exponential approach is able to reproduce airborne measurements of atmospheric wakes appropriately Platis et al., 2020Platis et al., , 2021. Similar to recent studies, we use the exponential approach of the E10 model to describe the wind speed deficits on the lee side of wind farms. Here, however, we apply the wake formulation in the spatial domain. At this, the wind speed magnitude is prescribed to decrease the strongest close to the offshore wind farms and recovers exponentially over the downstream distance. In summary, the velocity deficit u along the downstream distance x is given by: with α as the maximum relative deficit and σ as the exponential decay constant.
Given that the E10 model is established for turbine wakes at hub height, some adjustments are required for the determination of the wake deficit α and the decay constant σ. The individual values for α and σ, in particular near the sea surface, depend on complex relations between multiple aspects, such as the wind field, atmospheric stability, vertical momentum fluxes and wind farm density as well as the wind turbine drag. By the knowledge of the authors of this manuscript, there is so far no empirical or analytical formulation, which either considers these aspects or provides respective values for wind deficit and wake length. Thus, we selected typical mean values for α and σ based on measurements of recent studies ( Table 1) and computed SAR data statistics (Figure 2). Thereby, the decay constant σ is defined by the wake length.
For the SAR data statistics, we utilized latest data sets from Copernicus Sentinel 1 (downloaded April 2019) 1 , which consisted of satellite SAR Sentinel1-A and Sentinel1-B data from the period June 2016 to March 2019. We focused on the wind farm Global Tech in the German Bight, which is composed of 80 turbines, and considered wakes from all wind directions, except for cases of superposition of wakes from neighboring wind farms. The velocity deficits were computed using the 10-m wind field from the Synthetic Aperture Radar (SAR). Indeed, SAR technology demonstrated its ability of measuring the surface winds at fine resolution. Thereby, the changes in surface roughness caused by the perturbation of wind speed are captured by the SAR system and modify the Normalized Radar Cross Section (NRCS), which is measured by the satellite. With this, the 10-m surface wind speed is derived from the surface roughness using a geophysical model function CMOD5N (Hersbach et al., 2007) and velocity deficits can be calculated relative to the undisturbed wind field. A detailed description of how to obtain wind speed deficits from satellite SAR data can be found in Djath et al. (2018); Djath and Schulz-Stellenfleth (2019). 1 https://www.esa.int/Applications/Observing_the_Earth/Copernicus/ In total, the SAR data set resulted in 16 undisturbed wake cases exceeding 15 km length, which are shown in Figure 2. The wake measurements show the highest deficits in the first kilometers behind the wind farm, which decrease over the downstream distance. For the set of profiles, the mean velocity deficit was calculated and fitted to the exponential wake recovery function [Eq. (2)] using the least square method. The exponential fit led to values for α and σ of 7.5% and 32 km, respectively. However, the estimated values for the wake dimension apply primarily for the offshore wind farm Global Tech and might vary due to different wind farm arrays or atmospheric conditions. Thus, we ultimately selected reasonable values of α = 8% and σ = 30 km for the simplified wake parameterization of an undisturbed wake structure. This assumption is supported by earlier studies (Table 1), as the estimated wake deficit and length agree with the mean of previous sea surface wake observations.

Wake Cross-Sectional Shape
The cross-sectional shape of an atmospheric wake can be described by a symmetric exponential function, which is scaled by the characteristic wind farm width L. This assumption was validated in recent airborne measurements, where several wind speed profiles were measured perpendicular to the wake at hub height and fitted to a Gaussian function . However, Djath and Schulz-Stellenfleth (2019) indicated that the cross-sectional shape of a wake at 10-m height is likely more distinct toward the wake edges. Hence, we chose an exponential decay constant of γ = L/3, in order to narrow the wake cross section at sea surface height in comparison to the hub height assumption. As for the former studies, L describes the wind farm width with respect to the respective wind direction and is calculated for each wind farm individually.
The adjusted Gaussian function [Eq.
(3)] produces similar cross-sectional shapes to the Kaiser window function, proposed by Djath and Schulz-Stellenfleth (2019), and strikes a balance between the recent cross shape proposals.

Model Setup
Due to their spatial extension, atmospheric wakes affect both the small-scale processes near the wind farms and the largerscale processes in the vicinity of the wind farms. In order to model the cross-scale wake effects, we developed a model setup, which enables the simultaneous resolution of large-scale ocean dynamics and smaller-scale processes near offshore wind farms. For this purpose, we used the Semi-implicit Cross-scale Hydroscience Integrated System Model [SCHISM, Zhang et al. (2016b)]. The SCHISM model is a hydrostatic model grounded on unstructured horizontal grids, which uses Reynolds-averaged Navier-Stokes equations based on the Boussinesq approximation. The model solves transport equations with a second-order Total Variation Diminishing (TVD) advection scheme and applies a higher-order Eulerian-Lagrangian method for momentum advection. At this, it uses semi-implicit time stepping and a hybrid finite-element/finite-volume formulation. A more detailed description of the SCHISM model and its capacities can be found in Zhang et al. (2016a,b).
Our model region covers the North Sea extending laterally from the British Channel in the South to the Norwegian trench in the North and thereby includes the major areas of current offshore wind energy production in the North Sea (Figure 1). The model setup includes the highly resolved coastlines as well as major islands and river estuaries of the surrounding mainland. We used a depth-depending horizontal grid cell resolution for the triangular unstructured grid, ranging from 500 m in shallow coastal areas to 5,000 m in the deep open sea. In total, the grid has approximately 161 K nodes and 311 K triangles. In the vertical direction, we applied a flexible LSC 2 [Localized Sigma Coordinates with Shaved Cells, Zhang et al. (2015)] grid, which is composed of a depth-depending number of localized sigma coordinates and considers in total 40 depth levels. The layer thickness increases gradually from about 2 m in the surface layers to 10 m in levels below 100 m depth.
At the open boundaries, the model was forced by the North-West European shelf ocean physics reanalysis data from the Copernicus Marine Service (downloaded July 2019) 2 , which we interpolated for time series of elevation, velocity, temperature, and salinity. Additionally, we prescribed tidal amplitudes and phases of eight tidal constituents (M 2 , S 2 , K 2 , N 2 , K 1 , O 1 , Q 1 , and P 1 ) from the data-assimilative HAMTIDE model (Taguchi et al., 2014). For the atmospheric forcing, the coastDat-3 COSMO-CLM ERA interim atmospheric reconstruction (HZG, 2017) was applied, while daily river discharge values were provided by the mesoscale Hydrologic Model [mHM, Samaniego et al. (2010), Kumar et al. (2013)] using the E-OBS18 temperature and precipitation data set (Cornes et al., 2018). Werner et al. (in prep.) will provide a more detailed description of the discharge data. Bottom roughness was set to 0.1 mm and we used a vertical background diffusivity of 10 −6 m 2 /s. Initial conditions were interpolated from the same data sets as the boundary forcing data.

Model Simulation
After an initial spin-up simulation of two and a half years, we simulated offshore wind farm operation for the summer period of 2013 (May to September) in a subsequent model run. As the characteristic "flushing time" of the North Sea is about one year (Otto et al., 1990), the short spin-up period can be deemed sufficient for the model area. While the period of the simulations was generally chosen based on data availability, the simulation period for the wind farm simulations was chosen to ensure mostly stable atmospheric conditions, which enable the formation of large-scale wake structures. Thus, the simulations allow the investigation of the impact of wind farm wakes on the summer stratification. For all simulations, an implicit time step of 120 s was used and the daily mean was calculated for the output data.
The wind farm simulation took into account the recent state of offshore wind farm production within the model domain (see Figure 1). However, the wind farms were only considered via the implemented wake parameterization and were not represented physically in the horizontal grid. For the simulation, we applied the wake parameterization to the wind field interpolated onto the horizontal model grid and reduced the horizontal components of the wind speed respectively. Specifically, the wind speed on downstream grid points was iteratively reduced for each wind farm. While the parameterization modified the downstream wind field at each time step of the simulation, we considered a wind speed limitation (cut-in and cut-out speed) for wind farm operation in the wake model. For this, we used a typical wind speed range of 3-25 m/s at hub height, which inhibits wake generation for too weak and too strong wind speeds.
We enhanced the horizontal grid cell resolution in the wind farm simulation to resolve hydrodynamic processes near the offshore wind farms sufficiently and to provide highly resolved wake effects. The grid optimization involved a boundary layer of 30 km around each wind farm polygon with a grid cell size of 1,000 m (see Figure 3). In addition, a high-resolution zone with a grid cell size of 500 m was defined in the first five kilometers around each wind farm, where the wind deficits are most pronounced. In total, the adjusted horizontal grid of the wind farm simulation resulted in approximately 278 K nodes and 544 K triangles.
In order to compare differences between the wake-exposed North Sea and the reference case we ran the summer period twice, with and without wind farm parameterization. For both cases, the enhanced grid cell resolution was applied.

Model Validation
In order to validate the model simulations, temperature and salinity values were downloaded from the ICES database (downloaded March 2020) 3 for the period from January 2012 to December 2013. The observational data set consists of about 3,250 stations, at which temperature and salinity samples were collected at different water levels. For the validation, we calculated the depth-averaged data to avoid uncertainties in depth levels and compared it to the respective daily-mean values of the spin-up simulation, which forms the basis of the wind farm scenarios.
We utilized a Taylor diagram (Taylor, 2001) to display the statistics of our model simulations. The diagram is based on the law of cosine and quantifies differences between observation and model by the centered root mean square difference E' and the correlation coefficient R. The root mean square difference and the correlation coefficient are calculated as: with σ X = 1 N N n = 1 X − X 2 as the standard deviation, M as model data and O as observation data. Similar to Daewel and Schrum (2013), we subdivided our model domain into different validation areas based on their physical characteristics ( Figure 4A). Here, we focused on five main characteristic regions: The Atlantic inflow area (D), the British coast (E), the shallow tidal-impacted coastal areas connected to the Wadden Sea (F&G), the Baltic in-and outflow area (H) and the seasonally stratified central North Sea (K&L). For the statistics, we treated each region individually.
The calculated Taylor diagram ( Figure 4B) shows high correlation (R > 0.9) between the spin-up simulation and the observations for sea surface temperatures. Seven out of eight subareas exhibit correlation coefficients above 0.95, from which the areas connected to the Wadden Sea (F&G) show the highest values (R > 0.99). Only region K, the deep central North Sea area, shows a slightly lower correlation, which is probably due to the influence of the large northern open boundary of the model domain. Surface salinity has a lower, but still decent, overall agreement between the model and the observational data. Most coefficients range between 0.8 and 0.9, whereas only region E, the British coast region, and again region K show lower correlations. Generally, the lower correlation in salinity likely results from the comparison of daily-mean model data and instantaneous observational data, as we observe strong fluctuations in the fresh water inflow due to the tidal signal. Averaging the model data results in attenuated tidal fluctuations and thus lower correlation to the instantaneous values of the observational data. The British coastal area in particular, has strong tidal amplitudes, which could explain the particularly low correlation in region E.
In general, the model proves a sufficient correlation to the insitu measurements and is able to reproduce seasonal stratification in the southern North Sea ( Figure 4A). The pattern of the mean potential energy anomaly, which is a gravitational-based measure of stratification strength of the water column (Simpson and Bowers, 1981), shows expected characteristics, such as strong stratification in the deep central North Sea as well as mixing fronts around the shallow Dogger Bank and along the tidalimpacted Wadden Sea area. Thereby, the stratification patterns (exemplarily shown for August, Figure 4A) agree well with results of former North Sea studies (Schrum et al., 2003;Holt and Proctor, 2008) considering that magnitude and spatial dimension of summer stratification vary between years (Schrum et al., 2003).

RESULTS AND DISCUSSION
The analysis of wake effects in this manuscript concentrates on general processes related to the transfer of atmospheric momentum into the water column and associated mean changes in horizontal currents and stratification. In order to emphasize the induced changes during the summer period of 2013, wake effects are shown as the differences between the case of wind farm operation and the case without offshore wind farms in the southern North Sea. In the following, the data is primarily depicted as monthly means, while only a selection of figures is depicted in the analysis. For further insights on the monthly means, see Supplementary Figures A1-A6.

Spatial and Temporal Variability in Hydrodynamics
The implemented wake parameterization iteratively generates leeward wind speed reductions at each time step of the simulation. Despite the constant values for the wake deficit and wake length prescribed by the wake model, the resulting wake structures vary individually in size and intensity due to the different sizes of the respective wind farm arrays and due to superposition of neighboring wakes. A snapshot of a strong wind event during the simulation and an example of the associated downstream velocity deficits are depicted in Figure 5A. The emerging wakes extend along the present wind field and show comparable patterns to recent satellite SAR wake observations (Djath et al., 2018;Siedersleben et al., 2018). At this, the simulated wakes exhibit relative deficits around 10% in short downstream distances, which agrees with the observations of recent satellite SAR data ( Table 1). The strongest reductions in wind speed are observed in densely built areas of the domain, in particular in the German Bight. Although the wake model prescribed a constant deficit of 8%, superposition of neighboring wakes in densely built areas results in even higher velocity deficits of up to 15%, as multiple wind farms affect the wind speed.
As a result of constantly changing wind directions, pronounced wake patterns disappear when averaging over time and shape continuous zones of reduced wind speed around the respective wind farms. Thereby, the German Bight is the most affected region of the model domain, as a large number of current offshore wind farms is located off the Danish, German and Dutch coast. The mean deficits in wind speed of the 5 months of simulation are shown in Figure 5B. The deficits were calculated for the absolute values of the wind vectors. Figure 5B shows that the mean wind speed deficits primarily adapt to the predominant wind field during the simulation period, which here points into eastward direction. Furthermore, strong wind events, which might differ from the mean wind directions, can influence the shape of the wind speed anomalies particularly in areas of weaker winds (e.g., the British coast). Between May and September, the total absolute mean changes in wind speed range on average around 1-2% relative to the mean undisturbed wind field and can reach even beyond 5% in areas of clustered wind farms, like the Southern Bight or German Bight. Compared to the mean undisturbed wind speeds, which are between 5.0 and 7.0 m/s, the deficits account for an order of magnitude of about 0.1 m/s and range between 0.1 and 0.4 m/s. The order of magnitude is in agreement with recent atmospheric modeling results of Akhtar et al. (2021), where deficits near the sea surface range around 0.5 m/s during the summer months, considering that Akhtar et al. (2021) applied more influential future wind farm scenarios. Nevertheless, here, with regard to the interpretation of the magnitude of mean wind speed perturbations, the supposed over-and underestimation of the wake dimension due to the prescribed constant wake parameters should be taken into account.
The wake-induced anomalies in wind speed cover large areas around the actual wind farm locations. Consequently, wide areas are affected by reduced shear forcing and reduced vertical momentum fluxes at the sea surface boundary. This becomes particularly apparent looking at a cross section through the strongly influenced southeastern part of the German Bight ( Figure 5C). Along the cross section the percentage change in wind speed and wind stress is shown relative to the undisturbed values from the scenario without wind farm operation. Both, wind speed and wind stress, exhibit continuous deficits over several tens of kilometers, which do not fully recover in between adjacent wind farms and are particularly strong in the case of superimposed wakes. Thereby, the impact on the wind stress, which is the decisive factor for vertical air-sea exchange, is about twice as high in percentage as for the mean wind speed. Within the wind farm areas deficits slightly recover, since the wake parameterization reduces wind speed only downstream of wind farms. A similar pattern of the continuous influence arising from superimposed wakes was also shown by recent high-resolution atmospheric modeling of offshore wind farm effects (Akhtar et al., 2021), where moreover the effect on offshore wind energy production itself was addressed.
In this study, the focus is on the large regions of attenuated wind forcing and related impact on the marine environment. In the shallow southern North Sea, where for the most part wind determines the general circulation (Sündermann and Pohlmann, 2011), the reduction of wind speed and associated wind stress results in substantial changes of the residual circulation. As less shear force is acting on the sea surface in areas around the wind farms, vertical momentum fluxes become attenuated and less momentum is transferred from the atmosphere into the ocean. This entails potential changes in wave formation and horizontal surface velocity but also potentially reduces turbulence within the surface mixed layer and impacts stratification. Figure 6A depicts the mean changes in horizontal surface velocity after the first month of wind farm operation. The pattern shows expected lateral gradients in the horizontal currents as a result of the reduced wind stress behind offshore wind farms. At this, negative anomalies arise at wind farms in direct correlation to the downstream wind speed deficits of the predominant mean wind field (Supplementary Figure A1). In consideration of the prescribed wake intensity (8% deficit, 30 km length) and the prevailing wind speeds, the deficits in surface current speed range around −0.0025 m/s and exhibit peak values beyond −0.005 m/s. While the order of magnitude agrees with former wind farm studies (Ludewig, 2015), the velocity changes account for up to 5% compared to the mean surface residual velocity in May (about 0.1 m/s). The induced changes are up to 10-25% of the interannual and decadal surface velocity variability, as the mean surface velocity in the southern North Sea ranges between 0.11 and 0.13 m/s (Daewel and Schrum, 2017). Consequently, wake-induced anomalies in surface velocity can be quite substantial, considering that instantaneous changes are generally even stronger, and might affect the horizontal transport and eventually the residual currents in areas of wind farm operation significantly.
In order to evaluate the impact of the wind stress reduction on the vertical structure of the horizontal momentum, mean vertical profiles at the wind farm locations are depicted in Figures 6B-D. For the vertical profiles, we used grid points only located within a radius of about 5 km around each wind farm and calculated the average monthly mean changes within the water column. Thereby, shallow (z < 25 m) and deep (z ≥ 25 m) waters were treated separately, as they show different dynamics related to the wind shear. Figure 6B shows the average horizontal velocity profile at the offshore wind farm locations and the temporal variability during the month of May. As for the entire southern North Sea, the mean surface velocity at the wind farm locations is about 0.12 m/s, which agrees with the results of Daewel and Schrum (2017). However, the variability of surface velocity itself reaches up to 0.15 m/s within May and is much stronger than the interannual variability. In terms of wake-induced alterations in the horizontal velocity, the wind farm profiles show the strongest changes occurring at the sea surface boundary, where the impact by the wind speed reduction is most powerful (Figure 6C). Here, the maximum values range around -0.002 m/s at the sea surface, which is lower than the maximum changes in Figure 6A, since the wind farm locations are less affected by wakes than the areas in downstream direction. The order of magnitude of the changes accounts for about 2% of the mean horizontal velocity profiles at the wind farms as well as of the velocity variability over the month of May (Figure 6B). The wake-induced changes in horizontal velocity are more pronounced in deeper waters than in shallow coastal areas. Thereby, deep waters exhibit a strong vertical gradient between the surface and deeper layers, which emphasizes that particularly the wind-driven surface mixed layer is affected by the wind stress deficits. Nevertheless, the impact of the momentum deficit also cascades downward into deeper layers and diminishes toward the sea floor.
In addition to the perturbations in velocity, the attenuation of vertical momentum fluxes results in the associated reduction of turbulent mixing inside the surface mixed layer ( Figure 6D). Here, as a measure for the mixing rate, we depict the vertical eddy diffusivity K v . The profiles exhibit a reduced mixing rate over the entire water column. As for the horizontal velocity, the deficits in mixing are more pronounced in deep waters than in well-mixed shallow waters, which is likely favored by influence of the bottom mixed layer in shallow depths. In both cases, the strongest deficits occur near the pycnocline depth. While, in general, the apparent decrease in mean turbulent mixing implies a less turbulent water column, the maximum deficits near the mixed layer depth indicate a shallower surface mixed layer at offshore wind farms. This implies that the wind wakes counteract to the recently identified effects of Kármán vortices and turbulent wakes created by the pile structures of wind turbines, which are responsible for additional mixing in downstream direction (Carpenter et al., 2016;Grashorn and Stanev, 2016;Schultze et al., 2020). Thus, offshore wind farms induce counteracting processes, namely the reduction of horizontal momentum and the generation of turbulence, due to the wind turbine rotors and wind turbine foundations, respectively. However, the counteracting processes emerge on different spatial scales, since the pile effects remain primarily within the wind farm areas (Schultze et al., 2020).
Earlier studies showed the impact of the wind stress reduction on the sea surface elevation at offshore wind farms (Broström, 2008;Paskyabi and Fer, 2012;Ludewig, 2015). The studies demonstrated the formation of up-and downwelling dipoles under constant wind directions, which extend several kilometers around the wind farms. In a realistic model setup, such as we used, however, the wind field changes continuously over time and individual dipole patterns are expected to superimpose or mitigate, depending on the respective wind field. Figure 7A depicts the resulting mean surface pattern in August. Indeed, the changes in elevation superimpose into large coherent patterns of positive and negative anomalies, which smoothly merge into each other. Thereby, large-scale dipoles with spatial dimensions of up to hundreds of kilometers emerge in the wind farm regions and in particular in the German Bight area, where superposition and thus amplification of wake effects occurs most frequently. Compared to the mean wind field in August, which points into eastward direction (see Supplementary Figure A1D), the surface elevation dipoles exhibit their positive amplitude roughly on the windward side of wind farm areas and their negative amplitudes to the lee side. This might indicate a damming effect of winddriven horizontal flow by the offshore wind farms.
The changes in sea surface elevation show a clear correlation to the anomalies in mean depth-averaged velocity, which resemble the dipole pattern in the German Bight in August (Figure 7B). At this, the mean depth-averaged velocity changes range between ± 0.001 and 0.005 m/s and reach on maximum nearly 0.01 m/s. These changes are up to twice as strong as the wind stress related velocity changes at the sea surface in May ( Figure 6A) and account for up to 10% of the mean horizontal flow velocity in the German Bight. In percentage terms, such perturbations of the flow pattern are quite significant, as they are comparable to the potential changes due to climate change (Mathis and Pohlmann, 2014;Schrum et al., 2016) and may affect the residual currents as well as the gradients of the mean sea surface elevation in the German Bight area.
The relation between horizontal flow perturbations and sea level anomalies becomes more apparent by zooming in on the dipole structure in the German Bight specifically (Figures 8A,B). The figures illustrate the vector field of the mean depth-averaged horizontal velocity changes, in addition to the mean sea level and velocity changes in August. Figure 8A shows that the positive and negative mean sea level anomalies are clearly related to attenuated or enhanced horizontal flow velocities. At this, the sea level sinks in regions of increased flow velocity, while the sea level rises due to damming in regions of reduced flow velocities. This applies specifically to the large-scale dipole pattern. In the center of the German Bight, where the perturbations are most pronounced, the mean changes in sea surface elevation reach orders of magnitude of around ± 2 mm and result in spatial gradients of up to 4 mm over a distance of several tens of kilometers. These magnitudes are weak, compared to the local tidal variability and account for about 1-10% of the mean surface elevation variability in the German Bight, which varies by around 0.05-0.10 m over the respective spatial scales. However, the occurrence of such strong sea level gradients can in turn affect the geostrophic balance in the German Bight. Considering the apparent sea level differences of 2-4 mm, the resulting horizontal pressure gradients would affect the predominant circulation pattern by an order of magnitude of about 10 −3 m/s, which is similar to the wakerelated changes in horizontal velocity. In fact, these secondary induced perturbations of the geostrophic flow can be seen in the changes of the mean horizontal flow pattern (Figure 8B), where eddy-like structures with associated magnitudes of ± 0.006 m/s emerge in the center of the German Bight. However, these eddylike structures might also be the result of mean depth-averaged circulation changes, which in turn favor the emergence of up-and downwelling cells.

Impact on Stratification
The alterations of the sea level and the horizontal transport suggest changes in the lateral and vertical density fields. In a previous modelling study, Ludewig (2015) observed the occurrence of vertical upward and downward currents at offshore wind farms, which were inversely correlated to the displacements of the sea surface dipole. With sufficient temperature and salinity stratification, the induced vertical flow is expected to result on the one hand in advection of colder and saltier toward the surface and on the other hand in the downward advection of warmer and fresher water within the dipole area. Consequently, the changes in sea surface elevation are accompanied by changes in the vertical temperature and salinity distribution. Figure 7C depicts the mean sea surface salinity changes in August, where the magnitudes on average range between ± 0.05 g/kg. In the Southern Bight, the presumed correlation between the surface elevation anomalies and the surface salinity changes becomes most apparent. Positive and negative changes in surface salinity agree with the locations of the sea level perturbations and extend several tens of kilometers around the wind farms. At this, saltier surface water occurs in regions of lowered sea level, while surface salinity decreases in regions of elevated sea surface height. In the German Bight however, the large-scale pattern is more complex and does not clearly resemble the specific pattern of the mean sea surface elevation. Here, mean surface salinity decreases within the center of the large sea level dipole and along the shallow Jutland coast. Between the two negative salinity anomalies, salinity increases and the positive anomaly spreads along the fresh water plume from the Elbe estuary.
The mean surface temperature changes show much less agreement to the anomalies in sea surface height ( Figure 7D).
Instead, it appears that the surface temperature primarily increases in the vicinity of offshore wind farms, which confirms earlier simulations by Ludewig (2015). The apparent surface heating presumably results from the reduction in mixing and is hence directly related to the large-scale wind speed reduction in the atmospheric forcing. Over time, advection and lateral transport likely cause the temperature anomalies to spread over a large scale in the affected areas. Therefore, coherent patterns of increasing mean sea surface temperature are present in areas of wind farm development. At this, mean changes in sea surface temperature are on average around ± 0.02-0.05 • C, whereas changes in the German Bight reach even beyond 0.1 • C. The large-scale surface heating of up to 0.1 • C imitates the effects of climate change, in which an increase in sea surface temperature is also to be expected as a result of the warming of the earth's atmosphere (Schrum et al., 2016). However, the wake-related changes are about one order of magnitude smaller than the average perturbations due to climatic changes (Schrum et al., 2016). Furthermore, the changes account for a maximum of 10% of the annual and interannual variability of surface temperature in the southern North Sea, which are at least 1.0-1.5 • C (Daewel and Schrum, 2017).
Despite the deviations in the surface anomalies, the mean depth-averaged changes in temperature and salinity confirm the assumed correlation between the sea level changes and the changes in the density field in the German Bight (Figures 8C,D). The figures show again the vector plot of the mean depthaveraged horizontal velocity changes, in addition to the mean depth-averaged changes in temperature and salinity in August. For the depth-averaged data, salinity and particularly temperature exhibit stronger changes than for the sea surface with up to ± 0.1 g/kg and ± 0.3 • C, respectively. At this, the anomalies show a clear agreement with the vector plot of the depth-averaged velocity changes and with the mean sea level anomalies (see Figure 8A). In particular, the mean temperature changes strongly resemble the sea level alterations ( Figure 8D). Nevertheless, for both, temperature and salinity, a notable amplification of the changes occurs due to the eddylike structures in the horizontal flow anomalies and hence due to the amplified changes in sea level. These changes agree with the presumed correlation between the sea level anomalies and changes in the vertical density distribution and thus with the results of Ludewig (2015). However, associated changes in the mean vertical velocity field are not detectable.
The identified changes in temperature and salinity suggest an impact on the seasonal summer stratification in the southern North Sea. In order to investigate the impact, we calculated the potential energy anomaly as a measure for the stratification strength. Figure 9 depicts monthly mean changes in potential energy anomaly and the position of the tidal mixing fronts for the months of June to August 2013. For the tidal mixing fronts, the transition from mixed to stratified waters was defined by a surface-to-bottom temperature difference of 0.5 • C, similar to Skogen et al. (2011). Alterations in the potential energy anomaly can be seen in all summer months, which occur most strongly in the German Bight (Figure 9). The coherent changes show again large dipole patterns, which are related to the anomalies in sea surface elevation and consequently to temperature and salinity anomalies. What is notable here is that the changes in potential energy anomaly primarily occur in stratified waters. The apparent sensitivity of the wake effects to stratification was already seen in the influence on the mixing (Figure 6C), where the impact was much stronger in deeper and hence likely stratified waters.
The comparison of the monthly means shows a seasonal trend in the intensity of the impact on the potential energy anomaly (Figure 9). Initially, as the summer stratification evolves (Figures 9A,B), the patterns have a strong correlation to the spatial changes in the sea surface temperature and salinity (see Supplementary Figures A4, A5). At this, the changes in potential energy anomaly range on average around ± 2 Jm −3 . However, during the decline of stratification toward autumn ( Figure 9C), strong magnitudes of beyond ± 6 Jm −3 arise in the German Bight, which account for about 10-20% of the actual mean stratification strength ( Figure 4A). Thereby, the strongest changes in the potential energy anomaly arise in the center of the surface elevation dipole and in regions of declining stratification toward the Frisian coast. While the alterations within the dipole are expected to result from the dipole-related processes, changes near the tidal mixing fronts are presumably associated to the reduction of wind stress and hence of the reduction of mixing in regions of declining stratification. As the mixing reduces, the wake effects sustain stratification in areas of wind speed reduction and thus lead to the amplification of stratification. Figure 10 depicts the vertical profiles at wind farms, similar to Figures 6B-D, for the mean mixing rate K v for the months of June to September. The vertical profiles support the assumption about the impact of the mixing alterations on the development of the summer stratification. As seen for the month of May (Figure 6B), the mixing at wind farms generally decreases in all months of the simulation as a result of the wind speed reduction. Thereby, the most pronounced changes occur around  the mixed layer depth and in deep waters. In comparison of the different months, however, it is notable that the strongest reductions in mixing occur during the formation (May and June) and the breakdown (August and September) of the summer stratification. Especially toward autumn, the changes are up to five to ten times stronger than in the other months. At this, the maximum deficits generally occur near the mixed layer depth, which emphasizes the enhancement of stratification and the shallowing of the surface mixed layer at offshore wind farms due to atmospheric wake effects.

CONCLUSION
As wind turbines extract kinetic energy from the wind field, atmospheric wakes result on the leeward side of offshore wind farms and affect the far field wind speed near the sea surface boundary. However, wind farm wakes and their impact on the hydrodynamic system represent a multivariate system, which is highly sensitive to external influences in atmosphere and ocean. In this study, we developed an empirical approach to describe the wake-related wind speed deficits in a hydrodynamic model. With this, we provided first insights about the interference of single wake effects due to the recent offshore wind energy production and resulting larger-scale disturbances in hydroand thermodynamics in the southern and central North Sea. The implemented wake model generated wakes with constant dimensions, which, in terms of the different wind farm properties (e.g., the number of turbines and turbine density) or the changing atmospheric conditions, lead to over-and underestimation of the actual wake dimensions. Although the actual wake dimensions might vary in reality, the basic impact of wind farms, namely the downstream reduction in wind speed, is sufficiently considered to investigate the general consequences of wind speed deficits due to offshore wind farms. However, a more specific and realistic formulation will be important for more detailed future studies, since the wake dimensions determine the intensity of the associated impacts.
Over time, the extraction of energy by offshore wind farms results in extensive areas of reduced wind speed and subsequently the decrease of the shear-driven forcing at the sea surface boundary. As this reduces the momentum transfer from the atmosphere into the ocean, horizontal velocities and turbulent mixing initially decrease several tens of kilometers around offshore wind farms. Thereby the induced perturbations imply significant changes for the residual currents in the respective areas. Furthermore, convergence and divergence of water masses lead to the formation of sea surface elevation dipoles, which over time merge into large coherent structures. As shown here, these large-scale anomalies in the sea surface elevation are one of the main drivers of wake-related processes in the ocean. In addition to the general reduction of turbulent mixing, the largescale sea level alterations trigger lateral and vertical changes in the temperature and salinity distribution and affect the hydrodynamics in areas covered by offshore wind farms. However, the magnitude of these changes is rather small compared to the long-term variability of temperature and salinity and can hardly be distinguished from the interannual variability. A severe overall impact by the wake effects on the ocean's thermodynamic properties is thus not expected but rather large-scale structural change in the stratification strength and unanticipated mesoscale spatial variability in the mean current field. Nevertheless, further investigations are necessary to assess possible feedback on the air-sea exchange and thus potential impact on the regional atmospheric conditions, since surface heating along with the reduction in turbulent mixing influences the upward heat and momentum fluxes from the ocean into the atmosphere.
In this study, the structural changes in stratification become noticeable in a couple of ways. Firstly, we observed large dipole-related changes in the potential energy anomaly, as the geostrophic and baroclinic changes alter the temperature and salinity distribution. Secondly, the reduction of mixing at offshore wind farms results in the enhancement of the stratification strength, in particular, during the decline of the summer stratification. While the structural changes in stratification are minor in shallow mixed waters, the pronounced alterations in stratified waters can translate to the mixed layer depth, which likely increases or decreases depending on the respective stratification changes. This, in turn, might be crucial for marine ecosystem processes (Sverdrup, 1953). During the stratified summer months, the mixed layer depth is acting as barrier for nutrients and phytoplankton and plays a major role for the ecosystem dynamics. Therefore, induced fluctuations of the mixed layer depth can entail the intrusion of nutrients from the pycnocline into the surface mixed layer or the spreading of the nutrient-poor surface layer, respectively. The alterations in the nutrient availability, in turn, might affect local primary production and the nutrient balance. Thus, further studies are required to elucidate the impact on marine ecosystems and organisms in the North Sea, with regard to current and future wind farm scenarios.
Although the simulation was limited by the summer season of the year 2013, this study provides already important knowledge about the magnitude of implications of extensive offshore wind farming in the marine environment. Ultimately, the results showed that, in addition to the initial large-scale reduction in dynamics, wake effects cause internal processes, which affect the lateral and vertical transport of heat and salt even far beyond the associated wind farms. At this, changing wind directions inhibit severe local impact, as the varying internal processes mitigate each other. For the current state of wind energy production in the southern North Sea, the impact by the offshore wind farms entails structural changes in stratification strength and perturbations of the residual currents. Thereby, the cascading hydrodynamic processes particularly affect areas of clustered wind farms, like the German Bight. However, future wind energy development includes large installations in the Southern Bight and especially near the Dogger Bank, which will increase number of wind farms and hence possibly the impact in those areas (see Akhtar et al., 2021). Therefore, further investigations and improvements of the wake model are crucial, in order to assess the actual impact of future wind farm development in the North Sea environment. Especially since stratification anomalies entail perturbations of the pycnocline, assessment of potential biogeochemical impacts is of major interest.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
NC generated the data set, performed the data analysis, and wrote the manuscript. BD performed statistical analysis and contributed to the parameterization section of the manuscript. UD and CS provided essential background knowledge and adviced to the study. CS initiated the idea of the study. All authors contributed to conception and design of the study, manuscript revision, read, and approved the submitted version.