the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Seismic radiation from wind turbines: observations and analytical modeling of frequencydependent amplitude decays
Fabian Limberger
Michael Lindenfeld
Hagen Deckert
Georg Rümpker
In this study, we determine spectral characteristics and amplitude decays of wind turbine induced seismic signals in the far field of a wind farm (WF) close to Uettingen, Germany. Average power spectral densities (PSDs) are calculated from 10 min time segments extracted from (up to) 6 months of continuous recordings at 19 seismic stations, positioned along an 8 km profile starting from the WF. We identify seven distinct PSD peaks in the frequency range between 1 and 8 Hz that can be observed to at least 4 km distance; lowerfrequency peaks are detectable up to the end of the profile. At distances between 300 m and 4 km the PSD amplitude decay can be described by a power law with exponent b. The measured b values exhibit a linear frequency dependence and range from b=0.39 at 1.14 Hz to b=3.93 at 7.6 Hz. In a second step, the seismic radiation and amplitude decays are modeled using an analytical approach that approximates the surface wave field. Since we observe temporally varying phase differences between seismograms recorded directly at the base of the individual wind turbines (WTs), source signal phase information is included in the modeling approach. We show that phase differences between source signals have significant effects on the seismic radiation pattern and amplitude decays. Therefore, we develop a phase shift elimination method to handle the challenge of choosing representative source characteristics as an input for the modeling. To optimize the fitting of modeled and observed amplitude decay curves, we perform a grid search to constrain the two model parameters, i.e., the seismic shear wave velocity and quality factor. The comparison of modeled and observed amplitude decays for the seven prominent frequencies shows very good agreement and allows the constraint of shear velocities and quality factors for a twolayer model of the subsurface. The approach is generalized to predict amplitude decays and radiation patterns for WFs of arbitrary geometry.
 Article
(12300 KB)  Fulltext XML

Supplement
(600 KB)  BibTeX
 EndNote
In recent years, debates on the emission of seismic waves produced by wind turbines (WTs) and its potential effects on the quality of seismological recordings have led to increased research efforts on this topic. The main objectives are the characterization of WTinduced seismic signals, the definition of protection radii around seismological stations, and the modelingbased prediction of WT effects on seismological recordings in advance of the installation of WTs. Styles et al. (2005) reported about discrete frequency peaks in seismic noise spectra that increase with wind speed and the rotation rate of a nearby WT and assigned the observed peaks to vibration modes of the WT tower and rotor rotation. Zieger and Ritter (2018) and Stammler and Ceranna (2016) confirmed discrete frequency peaks between 1 and 10 Hz and analyzed signal amplitude decays with distance to the WTs described by a power law. Saccorotti et al. (2011) observed seismic signals with a frequency of about 1.7 Hz that were associated with WTs at distances of up to 11 km. Friedrich et al. (2018) used a migration analysis to identify seismic signals from nearby wind farms (WFs) and were able to distinguish between various WFs based on differences in frequency content. Polarization analyses was used by Westwood and Styles (2017) to show that Rayleigh waves dominate the wave field emitted from WTs. This observation was confirmed by numerical simulations (Gortsas et al., 2017). The increase of the noise amplitude with the square root of the number of WTs ($\sqrt{N}$) was observed by Neuffer et al. (2019) based on WT shutdown tests. Lerbs et al. (2020) proposed an approach to define protection radii, e.g., 3.7 km around the Collm Observatory (CLL) in Germany, using a power law to describe the spatial wave attenuation. Furthermore, the ground motion polarization near a single WT was analyzed and provided insights into the interaction of WT nacelle movement and emitted seismic signals.
Approaches to model the seismic radiation from WTs are rare and focus mostly on modeling the ground vibration of a single WT (Gortsas et al., 2017) or its operational components only (e.g., Zieger et al., 2020) but not on wave field propagation considering superimposed wave fields and amplitude decay with distance to multiple WTs simultaneously. However, Saccorotti et al. (2011) used an analytical approach to model the observed amplitude decays by summing up the calculated noise amplitudes produced by several WTs, including an intrinsic attenuation law, but they did not study possible effects of multiple WTs on the interference of the emitted wave fields.
In this paper, we present an analytical approach to model frequencydependent seismic radiation and amplitude decays with distance in comparison to robust longterm observed decay curves, measured at a WF in Uettingen (Bavaria, Germany). In a first step, we derive distancedependent noise spectra from recordings of up to 6 months in duration and characterize the relation between signal frequency and amplitude decay. We face the challenge of handling phase differences between multiple source signals that have strong effects on the seismic radiation due to significant changes in interference pattern of the superimposed wave fields. We apply the phase shift elimination method (PSE method) to generate representative source signals as an input for the analytical modeling of the observed amplitude decays. The comparison between modeled and observed amplitude decays also allows the constraining of the parameters of a simple twolayer model of the subsurface. We further show how it is possible to generalize the approach to predict radiation patterns for arbitrary WF geometries.
Our surveys were conducted in the neighborhood of a WF in Uettingen, about 9 km west of Würzburg in Bavaria. The WF consists of three WTs positioned in a NW–SE line with a spacing of 350 m and 450 m, respectively. The Nordex N117 type WTs have 2400 kW rated power and a tower height of 141 m. Their maximum rotation rate is about 12 rpm (rotations per minute). To measure the amplitude decay of the seismic WT signals we deployed 19 seismic stations along a profile of 8.3 km length, starting at the easternmost of the three WTs and running in NE direction approximately perpendicular to the geometrical layout of the WF (Fig. 1). Additionally, we placed three stations in the WT basements in order to record the seismic source signal of each WT. The instruments were installed between July and November 2019, and data recording will extend until August 2021. All stations are equipped with Trillium Compact posthole sensors (20 s) and Centaur data loggers (Nanometrics) recording continuously at a sampling frequency of 200 Hz. To improve the signal and noise conditions the sensors of the profile stations were placed in shallow boreholes of 1–2 m depth.
The local nearsurface geology is defined by Triassic sedimentary rocks. Beneath a thin soil layer, limestones of the Muschelkalk are situated over clastic sediments of the Buntsandstein, mainly terrestrial quartzite, sandstone, and claystone layers. Geologic cross sections suggest that the lower Muschelkalk under the topographic surface reaches a thickness of up to several tens of meters (Bayerisches Geologisches Landesamt, 1978). However, at some seismic stations the Muschelkalk–Buntsandstein boundary is only a few meters below the surface. In topographic depressions the Muschelkalk can be completely missing, i.e., thin quaternary soft sediments directly cover Upper Buntsandstein rocks.
2.1 Calculation of average power spectral densities (PSDs)
We analyzed a continuous dataset between September 2019 and March 2020, covering a range of 159 to 207 d depending on the exact station installation date. We associate the measured amplitudes in the seismic waveform data with the corresponding WT parameter (in this case “rotor speed”) at a resolution of 10 min. For this reason, the recordings of each profile station were split into 10 min segments that were transformed to power spectral density spectra using the method of Welch (1967). Each of these spectra were then sorted according to the respective rotor speed into bins of 1 rpm width. With this procedure we generated close to 10 000 single PSDs within the bin of maximal rotor speed (11–12 rpm) called “full power” status for each station, and about 2000 single PSDs for the “zero power” status of the WT (0–1 rpm). In order to reliably remove outliers and reduce the impact of local transient noise (e.g., traffic on nearby roads), we excluded 75 % of the largest PSD amplitudes and used only 25 % of the single PSDs to calculate the final average PSD spectra. This seems to be a relatively strong limitation of the dataset. However, due to the long observation period there are still enough data left to calculate robust average spectra. We think that this approach provides a reliable and conservative estimate of the spectral WT amplitudes with a minimized influence of interfering transient signals. Figure S1 in the Supplement illustrates the influence of different percentiles on the calculated average PSD at station F01.
Figures 2 and 3 show the resulting average PSDs (25 % percentile) for the full power and the zero power WT status, respectively. Besides the strong microseismic peak at about 0.2 Hz, we identified nine peaks of significant energy centered at 1.14, 1.7, 2.3, 3.5, 4.8, 6.0, 7.6, 10.5, and 17.2 Hz. All of them show a systematic amplitude decrease with increasing station distance, indicating that their origin is located at the WT. For peaks 1 to 7 we fitted the observed amplitude decay with a power law model (see next section). Because of the rapid amplitude decay at frequencies >10 Hz, we were not able to reliably fit peaks 8 and 9. For comparison we show the respective average PSDs recorded during zero power status (0–1 rpm) in Fig. 3. In this case the observed spectral peaks 1 to 9 have completely disappeared. The remaining (sharp) peaks show no systematic dependence with distance to the WF or the rotation rate of the WTs, which is an indication that their origin is not related to the WTs. Highfrequency signals at >4 Hz have a large amplitude in the near field only. Similar peak distributions have been observed by Neuffer et al. (2019) and Lerbs et al. (2020).
2.1.1 Powerlaw fitting of the observed amplitude decay
To quantify the PSD amplitude decay, the respective peak maxima of the full power PSDs (Fig. 2) were picked at each station. Figure 4 shows the resulting attenuation curves for peak 1 (1.14 Hz) to peak 7 (7.6 Hz) using a doublelogarithmic representation, i.e., the logarithm of peak amplitude is shown versus the logarithm of the station distance. If the PSD amplitude decay corresponds to a power law, which is the basic assumption, there should be a linear correlation between log(amplitude) and log(distance). The attenuation factor, b, can then be calculated as the slope of a linear fit of the attenuation curves. As Fig. 4 shows, the measured PSD amplitude decays can be described in good approximation with a power law between station F02 and F12, which corresponds to a distance range of 300 to 4000 m. Beyond F12, the measured PSD amplitudes increase with larger distances, except for peak 1 (1.14 Hz), and it was not possible to identify clear peak maxima, since the background noise dominates the spectra. Towards the end of the seismic profile the stations get closer to a highspeed railway track and populated areas with raised ambient noise conditions, which might explain the observed excessive PSD amplitudes in this region. However, the first two stations of the profile – F01 (194 m) and F01a (239 m) – also show deviations from a power law attenuation. Due to the proximity of these stations to the WT, the amplitudes may also be affected by nearfield effects.
For these reasons we decided to restrict the analysis of the amplitude decay to the distance range between 300 and 4000 m and to estimate the attenuation factor, b, using a linear leastsquares fit between station F02 and F12 for all seven peak frequencies. The results show a systematic increase with frequency and yield values from b=0.39 at 1.14 Hz up to b=3.93 at 7.6 Hz. In Fig. 5 we show the frequency dependence of the attenuation factor b. It exhibits a nearly perfect linear relationship between b and frequency, at least within the analyzed frequency range from 1.14 to 7.6 Hz. The comparison to results from other authors is discussed below (Sect. 5).
The b values quantify the decay of PSD peaks in the vicinity of the Uettingen wind farm. Since PSD amplitudes are proportional to the square of amplitudes in the time domain it is possible to estimate the corresponding b values for time signals by applying a factor of 0.5. Table 1 lists both types of attenuation factors: b_{PSD} for the PSD decay, and b_{AMP} for the corresponding time amplitude decay. This should enable a better classification of our results.
2.2 Observation of phase shifts between multiple WT vibrations
Each WT can be considered a seismic source. By analyzing the seismograms measured simultaneously in the three WTs (I01, I02, and I03) of the WF Uettingen, we observe phase shifts between the individual wave forms (Fig. 6a). As an example, three time series (vertical component), each recorded in one of the WT cellars during a rotation rate of about 11.5 rpm, are filtered to a narrow bandwidth around frequency peak 1 with 1.14 Hz (1.10–1.18 Hz) and are compared within a time window of 22 s. In the first 2 s, the signal phase of seismic station I03 is shifted by π compared to signal I01 and I02, which are in phase. Between 10 and 13 s, all signals are almost in phase, which consequently means that the WTs are vertically vibrating in phase. After 15 s, all three signals are shifted to each other and are not in phase anymore. For a longer time period of 1 h, the phase shifts between signals measured at I01 and I02 are determined using a crosscorrelation analysis with a moving time window of 5 s (=720 time segments) along the 1 h time segment. The temporal shift is converted to the corresponding phase shift between −π and π for each window. The distribution of all 720 resulting phase shifts is almost uniform (Fig. 6b) and shows no systematic behavior with time (Fig. 6c), which leads to the conclusion that phase differences between source signals appear rather randomly, especially over longer time periods.
In the following section, we model the observed amplitude decays and set up a mathematical formulation that includes a source function, attenuation factors, geological properties, and the superposition of multiple wave fields (produced by multiple WTs). In view of the observation that the source signals of neighboring WTs are not in phase, we study the influence of possible signal phase differences on the amplitude decay and propose a solution as to how to account for or “eliminate” this effect in the calculation.
3.1 Surface wave field approximation
Previous research suggests that mainly vertically polarized Rayleigh waves are emitted from WTs and that they dominate the WTinduced seismic noise (Westwood and Styles, 2017; Neuffer and Kremers, 2017; Gortsas et al., 2017). However, recent studies indicate that both Rayleigh and Love waves are emitted from WTs (see Lerbs et al., 2020; Neuffer et al., 2021). In our models, we assume that surface wave amplitudes decay proportionally to ${r}^{\mathrm{1}/\mathrm{2}}$ (with distance r to the source) due to geometrical spreading of the surface wave front on a cylindrical area in the 2D surface plane
where r_{0} is a reference (minimal) distance (Bugeja, 2011).
Geometrical spreading is independent of wave frequency. In addition, attenuation due to intrinsic absorption reduces the wave amplitude with distance to its source
The damping factor, D, depends on frequency w=2πf, seismic wave velocity c, and again the travel distance r of the wave (Bugeja, 2011). Furthermore, D is a function of the seismic quality factor Q, which describes the loss of energy per seismic wave cycle due to anelastic processes or friction inside the rock during the wave propagation. The damping of the wave is decreasing with increasing Q. The source signal S(t) itself is approximated by a continuous periodic cosinefunction to simulate the periodic motion at the base of the WT in vertical direction
where S(t) is a function of time t, signal frequency w=2πf, amplitude calibration factor A, wave number $k=w/c$, and signal phase Φ.
Assuming a homogeneous halfspace, the wave amplitude can be calculated for any distance r to the source (Fig. 7). Considering N source points (WTs), the amplitude at each point and hence the total wave field is derived by summation over all N wave fields
where index i corresponds to source point i and the relative radial distance to the source points is given by ${r}_{\mathrm{0},i}/{r}_{i}$. Shear quality factor Q_{S} and seismic shear wave velocity c_{S} are model parameters and define the properties of the material the wave is traveling through. Z(t) is the superposition of the individual wave fields and can be calculated at any time t.
By modeling the interference wave field as a function of time, this approach allows to derive rootmeansquare amplitudes (rms amplitudes) at any point at the surface. For the calculation, the amplitude calibration factor A_{i} will be set to 1 for every source signal, since all three WTs of the Uettingen WF are of the same type. It should be noted that body waves (P, S) are not considered in this modeling approach, since the simulated wave field approximates a wave field dominated by Rayleigh waves. The velocity of Rayleigh waves c_{R} is generally slightly lower than the shear wave velocity c_{S}, whereas the ${c}_{\mathrm{R}}/{c}_{\mathrm{S}}$ ratio depends on the Poisson ratio ν (e.g., Rahman and Michelitsch, 2006). Assuming theoretical values of ν from 0.0 to 0.5, the ratio ${c}_{\mathrm{R}}/{c}_{\mathrm{S}}$ can reach values between 0.87 and 0.95 (Leiber, 2003; Hayashi, 2008), which means that the Rayleigh wave velocity is maximally about 13 % lower than the shear wave velocity. However, it is possible to approximate surface wave fields using the shear wave velocity (Kumagai et al., 2020). The penetration depth of Rayleigh waves, which is influenced by the physical properties of nearsurface geological layers and the wavelength, plays an important role in approaching the modeling of WTinduced seismic wave propagation. The correct quantification of the penetration depth of Rayleigh waves is widely studied; however, so far there is no general consensus on their penetration depth in relation to the seismic wavelength λ. Based on results of Hayashi (2008), Kumagai et al. (2020) claim that surface wave velocity reflects the average Swave velocity of the geological layers down to a depth between $\frac{\mathrm{1}}{\mathrm{4}}\mathit{\lambda}$ and $\frac{\mathrm{1}}{\mathrm{2}}\mathit{\lambda}$, whereas $\frac{\mathrm{1}}{\mathrm{3}}\mathit{\lambda}$ is often chosen to be the most suitable assumption (e.g., Larose, 2005). Moreover, it is common to derive depth information from observed wave attenuation applying modeling or tomography methods to seismological data (e.g., Siena et al., 2014). Due to Rayleigh wave dispersion, it is known that lowfrequency surface waves reach deeper into the subsurface and thus will travel through materials with likely higher Q and seismic velocities c. Consequently, the damping is reduced compared to highfrequency surface waves (Karatzetzou et al., 2014; Farrugia et al., 2015). Taking this into account, we use the following relation for wavelength depth conversion
In this study, we take advantage of the link between frequencydependent amplitude decays (depending on Q_{S} and c_{S}, Eq. 4) and surface wave penetration depth to derive information about shear wave velocities and quality factors in the subsurface (Eq. 5).
3.2 Effect of sourcesignal phase on seismic radiation and amplitude decay
Since we observe significant changes of the phase shifts between signals measured at the three WTs (Sect. 2.3), we aim to study its effect on the wave field that is emitted by the three WTs in Uettingen. Hence, three wave fields are calculated using three different source phase compositions assuming a 1.14 Hz source signal frequency, 1500 m s^{−1} wave velocity and equality factor of 30 as an exemplary model. Source points are located at ${x}_{\mathrm{1}}={x}_{\mathrm{2}}={x}_{\mathrm{3}}=\mathrm{0}$ m and y_{1}=350 m, y_{2}=800 m, and y_{3}=1150 m (Fig. 7). Amplitudes are calculated along a profile extending from source signal S1 (see Fig. 7) and perpendicular to the WF line, which approximates the real geometry of the WF and seismic profile in Uettingen. The results show a clear dependence of the amplitude decays on the source signal phase composition (Fig. 8). In addition, amplitudes at the end of the 5000 m long profile differ significantly from each other. In the third scenario (Fig. 8c), the expected amplitude is a quarter of the amplitude that is reached if the WTs are vibrating in phase (Fig. 8a). Furthermore, strong effects appear in the first 2 km of the profile. Scenario (a) shows increased amplitudes due to constructive wave interference in the near field, whereas scenario (c) indicates a rapid decay of amplitudes within the first 1000 m of the profile. Scenario (b) shows a smoother and steadier decay of the amplitudes and reaches an amplitude at the end of the profile that is reduced by a fifth compared to scenario (a). These exemplary scenarios demonstrate only three out of infinite possibilities of different source signal phase compositions. Taking this into account, the seismic radiation of a WF is affected by phase differences of the source signals which can lead to strong changes in the wave field interference.
3.3 Phase shift elimination and data fitting
In this section we propose a method of how to handle the observed timevarying source signal phases and their effect on the seismic radiation using the assumption of a random appearance of signal phase constellation of multiple WTs, especially regarding long time periods. To define representative source signals, we developed a phase shift elimination method (PSE method). Within this PSE method, 500 radiation patterns (i.e., wave fields) are calculated using random signal phases Φ (between 0 and π) for each individual source signal. All 500 calculated wave fields are then averaged, and the average amplitude decay is extracted along the profile, which in turn is independent of the individual source signal phases. We experienced that the wave field averaging process and hence final amplitude decay calculation is sufficiently stable after 500 wave field simulations, whereas a number of <100 seems too low to generate a reproducible result. We apply this method to the Uettingen WF setup and compare the modeling results with the observed amplitude decays in Uettingen during rotation rates between 11 and 12 rpm (all WTs under full power). Since the observed PSD amplitudes are proportional to squared ground motion amplitudes in the time domain (rms amplitudes), we compare our modeling results (rms amplitudes) with the square root of the observed PSD amplitudes. The analysis is performed for signals with center frequencies of 1.14, 1.69, 2.26, 3.5, 4.85, 5.98, and 7.6 Hz, representing the 7 PSD peaks in Uettingen (see Sect. 2.1). For comparison, all decay curves are normalized to the amplitude measured in 300 m distance (seismic station F02), to be consistent with the attenuation analysis presented in Sect. 2.2. The calculated radiation pattern covers an area of 6000 m in length (x) and 1500 m in width (y) with a grid space of 10 m.
Calculated and observed data are fitted by a Q_{S}–c_{S} grid search to find the best model parameters. The data is grouped into signals with low frequencies <4 Hz (1.14, 1.69, 2.26, and 3.5 Hz) and high frequencies >4 Hz (4.85, 5.98, and 7.6 Hz) to distinguish between shallow and deep geological effects on the amplitude decay due to the frequencydependent penetrating depth of surface waves. All amplitude decays per group are fitted with one Q_{S}–c_{S} model. To set up the grid search, model parameter c_{S} are varied from 400 to 3000 m s^{−1} using steps of 20 m s^{−1}, and the parameter Q_{S} is varied between 6 to 250 using a step size of 2. An averaged (500 decay curves) amplitude decay with distance is calculated for each combination of Q_{S} and c_{S} and is compared to the observed data by calculating the rootmeansquare error (RMSE):
where M represents the 14 seismic stations along the profile that are included in the fitting process. This process is performed for 16 114 different models per frequency (Fig. 9).
Moreover, the normalized rootmeansquare error (NRMSE) is obtained to quantify the fitting quality. To determine the NRMSE, each RMSE is divided by the range (maximum value − minimum value) of the observation amplitudes for each frequency in order to scale the comparison between the datasets. The total NRMSE is then given by the mean of all normalized RMSE:
To fit modeled and observed amplitudes we performed a separate gridsearch for both the group with highfrequency signals and with lowfrequency signals. During each individual fitting process, the PSE method was applied to ensure results that are independent of source signals phases. Regarding the group of lowfrequency signals (<4 Hz), we obtain Q_{S}=40 and c_{S}=960 m s^{−1} (Table 2) as the best model parameters. The values of the 20 best models range between 36 and 60 for Q_{S} and 920 m s^{−1} and 1040 m s^{−1} for c_{S}. Regarding the group of highfrequency signals (>4 Hz), we obtain Q_{S}=16 and c_{S}=540 m s^{−1} as the best parameters. Results of the 20 best models range between 12 and 32 for Q_{S} and between 540 and 660 m s^{−1} for c_{S} (Fig. 10). By fitting two frequency groups, we can derive a twolayer model (one layer with a halfspace below), after converting the frequencydependent wavelength to the corresponding penetration depth (Eq. 5). Thus, we expect a shear wave velocity of 540 m s^{−1} down to 37 m depth and 960 m s^{−1} until 280 m depth (Fig. 11). However, transition between the two layers (37 m to 91 m) is not clearly defined due to missing information for frequencies between 3.5 and 4.85 Hz. Furthermore, we can only gain information about the attenuation of signals down to a frequency of 1.14 Hz. Hence, we are limited concerning the conversion of wavelength to depth information, and we can only derive values (c_{S} and Q_{S}) to a estimated depth of about 280 m. We cannot give information about the properties of deeper layers. The velocity error is approximated by the range of the 20 best models (Fig. 10). During the fitting process, we noticed that it was not possible to fit all seven amplitude decays with only one c_{S}–Q_{S} model successfully, especially regarding signals with frequencies >4 Hz. A homogeneous model is consequently not reasonable in this case. However, the corresponding results are given in the Supplement (Fig. S2). Modeled and observed data are generally in very good agreement for each of the seven analyzed frequencies (Figs. 12 and 13). The very slow decrease of observed amplitudes, especially at 1.14 Hz (Fig. 12a), and the relatively strong decrease of signals with 7.6 Hz (Fig. 13c) are simulated correctly and confirm a higher attenuation with higher frequencies, as expected. For a frequency of 1.14 Hz, between x=2000 and 4000 m, modeled amplitudes are underestimated in comparison to the observations. Minor deviations between modeled and observed data for frequencies >4 Hz might be explained by local effects that are not represented in our laterally homogeneous models. Interestingly, a local increase of amplitude with distance is observable in the real data, especially for the 3.5 and 4.85 Hz signals, as well as in the simulated data (Figs. 12d and 13a). This undulation is likely caused by superimposed wave fields of multiple WTs, as indicated by the modeled radiation pattern. Moreover, the sensitivity concerning the source signal phase compositions decreases clearly with increasing frequency, which indicates that effects of phase differences between source signals are more significant for lower frequencies. The radiation patterns off the profile are quite symmetrical since the WTs are positioned in a clear geometry (in a line with similar distances between the WTs). This might not be given if the WTs are arranged in more complex layouts that could lead to locally increased or decreased amplitudes due to wave field interferences.
The aim of this study is to present reliable amplitude decays of seismological signals produced by multiple WTs and to model these amplitude decays with an analytical approach. The propagation of WTinduced seismic signals has been the subject of numerous studies. Many authors found that the amplitude decay with increasing distance (r) between WT and observation point can be described by a power law of the form $\mathrm{1}/{r}^{b}$. In general, the absorption factor, b, increases with increasing frequency. Results found in our study show a nearperfect linear increase of the b value with frequency and range from b=0.39 at 1.14 Hz up to b=3.93 at 7.6 Hz. By converting the b values from the spectral domain (PSD) into the time domain (Table 1), we find b values <0.5 at frequencies ≤2.3 Hz, which is lower than the theoretical value for geometrical spreading of a single source surface wave. An explanation for this is the interference of the three wave fields from the three WTs. Although the decay of signals from a single WT would likely not be lower than the geometrical spreading attenuation, we can show that the superposition of wave fields could significantly increase the amplitudes along the profile. The b factors derived by the various authors cover a broad range of values, even for similar frequency ranges. Flores Estrella et al. (2017) published b values from 0.73 to 1.87 for frequencies between 2.7 and 4.5 Hz. Zieger and Ritter (2018) derived values from 0.78 to 0.85 at 1–4 Hz and b=1.59 at 5.5 Hz. The results from Lerbs et al. (2020) range between 0.7 and 1.3 at 1–4 Hz and b=2.3 at 5 Hz. Neuffer et al. (2019) derive b values of 2.4 at 3 Hz and values of b>5 at frequencies of 6–7 Hz. In Fig. 5, we compare our results with the b values of Neuffer et al. (2019) and Lerbs et al. (2020). These studies yield a similar frequency dependence. However, the results of Neuffer et al. (2019) show systematically higher b values. This observation could be due to different geological conditions with stronger attenuation effects during wave propagation. Furthermore, Neuffer et al. (2019) used socalled “differential” PSD spectra to measure the peak amplitude decay. These amplitudes are calculated from the difference between the PSD peaks at full power and the PSD peaks at zero power, which could lead to an overestimation of the amplitude decay. Lerbs et al. (2020) get similar b values; however, compared to our results the scatter is significantly larger. Most authors explain the observed bvalue scattering using different local geological conditions that influence the attenuation of the emitted seismic WT signals. It should be noted, however, that some of the abovementioned studies use relatively short time windows to estimate the spectral amplitudes at increasing distances. Flores Estrella et al. (2017) analyze time series of 2 h lengths, and Lerbs et al. (2020) use 6 h. In this case the measured amplitudes could be affected by transient signals such as earthquakes or local anthropogenic noise sources, which may result in uncertain bvalue estimates. In contrast, Neuffer et al. (2019) extend the analysis to 6.5 weeks. Since the knowledge of the amplitude decay plays a fundamental role in the modeling of the WT signals, we decided to use significantly longer time windows (6 months) in order to derive robust average PSD spectra at the installed profile stations.
In terms of modeling approaches, most of the recent publications focus on modeling the seismic signals that are emitted by one single WT (e.g., Gortsas et al., 2017) or the whole WF is considered as one emitting source. Since we observe timevarying phase differences between the signals that are measured directly at the three individual WTs of the WF Uettingen, we propose that this effect must be included in the modeling of WFs. Our observations confirm the significance of phase differences between the seismic signals from the WTs of a wind farm and that the signal phase of a single WT is not stable over time. Hence, we expect that phase differences between source signals vary randomly, which was already presumed by Saccorotti et al. (2011). Superimposed wave fields lead to constructive and destructive interferences (which depend on, e.g., signal phases) and affect the spatial amplitude decay, as we can show in this study (Fig. 8). Similar to our approach, Saccorotti et al. (2011) modeled amplitude decays on the basis of superimposed wave fields and attenuation laws but did not include phase shift variations between signals of the WTs. However, they noticed that the increase of noise depends on WT number, which was later shown by Neuffer et al. (2019). Saccorotti et al. (2011) suggest that more accurate results can be derived by considering WTs that are not vibrating in phase. Here, we can prove the randomness of these phase differences between WTs and propose a solution by applying the PSE method to the modeling. Only with this consideration we can reproduce the observed amplitude decay. The PSEmethod (averaging 500 wave fields calculated with random signal phases) is generally difficult to apply if full wave form propagation simulation is needed (e.g., FEM, finite element methods), since the required computation time would increase rapidly. Within our modeling approach, the source amplitude is chosen to be uniform for the three WTs. Previous studies (e.g., Lerbs et al., 2020) showed that WTs emit signals with timevarying amplitude and azimuths. In terms of modeling radiation patterns for very short time periods, this should be considered when choosing representative source characteristics. To model radiation patterns that represent long time periods (quasistatic processes), a uniform source amplitude should be sufficient, provided that the WTs are of the same type.
For the Uettingen WF the discrepancy between observed and simulated amplitude decays for 1.14 Hz in distances larger than 2000 m to the WTs are likely due to the other nearby WFs A and B (Fig. 1). We assume that the lowfrequency signals of these WFs travel farther compared to higherfrequency signals and are measured in addition to the signals from the targeted three WTs in Uettingen (Fig. 1). This could lead to an overestimation of the signal amplitudes, especially in the far field of the WF Uettingen. However, since we observe peaks at identical frequencies in the near and far field of the WF, it is reasonable to assign these signals mainly to the wave field produced by the WTs in Uettingen. Signals from various WFs can generally be distinguished using, e.g., a migration approach (Friedrich et al., 2018). However, detailed analysis of the effect of additional WFs around Uettingen is beyond the scope of this study, but their impact should be considered in future analysis. Interestingly, the sensitivity to source signal phases (gray lines in Figs. 12 and 13) is significantly higher for 1.14 Hz signals than 7.6 Hz and is generally decreasing with increasing frequency. This indicates that the signal phases are not as important for higher frequencies than for lower frequencies (e.g., 1.14 Hz). It should be noted that some of the individual input source signal phase compositions led to decay curves that could not fit the observation data at all. This is solved using the PSE method.
Lerbs et al. (2020) proposed a solution that describes the wave attenuation with distance using an attenuation model solely based on a power law assumption (b values). This approach does not allow a more universal application to other WFs or regions since b values are not directly assigned to geological properties. The approach used in our study includes the intrinsic attenuation factor, which depends on two geological parameters, the seismic wave velocity and quality factor. We find frequencydependent Q values of 16 and 40 and c_{S} of 540 and 960 m s^{−1}. The local geology is dominated by sedimentary rocks of the Buntsandstein, which could explain the relatively low Q values (high damping). The attenuation is very likely dominated by intrinsic attenuation and not scattering, since the topography around Uettingen is relatively smooth and large damaging zones or faults are missing. A homogeneous halfspace is therefore the basic assumption within our model. However, we show that the effect of layered media in the underground should be considered assuming frequencydependent velocity and quality factors, due to significant dispersion effects of surface waves. It is generally an advantage to include geological properties in the model: (1) to consider actual physical properties of the medium the waves are traveling through and (2) to enable the possibilities of studying the effect of various geological conditions on the seismic radiation and amplitude decays.
To demonstrate the capability and possible application of the modeling approach used in this study, we modeled the radiation pattern of the original WF in Uettingen for 1.14 and 7.6 Hz signals and compared the results with the case where three imaginary WTs are arbitrarily added to the existing WF (Fig. 14). Model parameters are c_{S}=960 m s^{−1} and Q_{S}=40 for 1.14 Hz and c_{S}=540 m s^{−1} and Q_{S}=16 for 7.6 Hz. The pattern of the radiation for 1.14 Hz signals is clearly affected by adding three WTs to the WF in Uettingen, whereby amplitudes are significantly increased, even in remote areas, for example in the NNW of the WF (Fig. 14a and c). The effect on the characteristic radiation of 7.6 Hz signals is negligible since the signal amplitude is damped rapidly in both cases, i.e., modeling three WTs or six WTs (Fig. 14b and d). As the demonstration shows, the modeling approach allows the estimation of the characteristic seismic radiation pattern of an arbitrary WF in order to identify locations of low or high noise amplitudes or to evaluate WF geometry effects. Furthermore, the source locations, source signal frequencies and amplitudes, and the expected local underground are free for researchers to choose (with limitations regarding their complexity), which enables an approximation of the surface wave field emitted by WFs with various layouts.
We recorded the seismic signals emitted from a threeturbine WF in Uettingen, Bavaria, over a period of 6 months and analyzed the spectral characteristics and spatial amplitude decays. During the full power operation mode of the WTs we identify seven prominent spectral peaks in the frequency range from 1.14 to 7.6 Hz. The attenuation of the peak amplitudes with respect to the WT distances can be described by a power law with exponent b. We find that the calculated b values increase linearly with increasing peak frequency and range between 0.39 and 3.93. Due to the relatively long observation period, the calculated values provide a stable basis for the analytical simulation of the emitted wave field.
An analytical approach was developed to model the seismic radiation of the WF. From measurements we observe that WTs are not vibrating in phase and that the phase differences vary randomly over time. Furthermore, the results of the simulation show a strong influence of phase differences between single WT source signals on the radiation pattern and hence on the spatial amplitude decays. We applied a phase shift elimination method (PSE method) to eliminate this effect with the aim of deriving a representative seismic wave field. Modeling results were compared to the observed frequencydependent amplitude decays to derive model parameters (Q_{S} and c_{S}) for a twolayer model that provides information about the local geology. Concerning the modeling of WTinduced seismic signals, we can show that the signal phases of multiple source signals (multiple WTs) have significant influence on the seismic radiation of the WFs. This effect should be carefully considered when selecting suitable source signals to avoid misleading simulation results.
The code and data used in this research are currently restricted.
The supplement related to this article is available online at: https://doi.org/10.5194/se1218512021supplement.
FL and ML performed the field work and data analysis. FL set up the modeling approach and performed model calculations. GR participated in data interpretation and model development and supervised the article outline. HD and GR initiated the project and provided the computational framework. FL, ML, GR, and HD edited the article.
The authors declare that they have no conflict of interest.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We would like to thank ESWE Versorgungs AG for providing access to the wind farm facilities in Uettingen and their operational data, and we would especially like to thank Ulrich Schneider for his support in initiating the project. We appreciate the support of the mayors of Uettingen, Greußenheim, Remlingen, and Leinach as well as their respective municipal administrations. Special thanks are given to Wolfgang Reinhart, Joachim Palm, and Ana Costa for their help during fieldwork and station maintenance and to Gabriele Schmidt (ESWE), who was our contact person regarding technical matters of the WF Uettingen.
This research is part of the project KWISS and has been supported by the German Federal Ministry for Economic Affairs and Energy (FKZ no. 0324360) and ESWE Innovations und Klimaschutzfonds.
This openaccess publication was funded by the Goethe University Frankfurt.
This paper was edited by Charlotte Krawczyk and reviewed by Joachim Ritter and Hortencia Flores Estrella.
Bugeja, R.: Crustal Attenuation in the region of the Maltese Islands using CodaWave Decay, PhD thesis, Department of Physics, University of Malta, 2011.
Bayerisches Geologisches Landesamt: Geologische Karte von Bayern 1:25.000, Blatt 6124, Remlingen, München, 1978.
Farrugia, D., Paolucci, E., D'Amico, S., Galea, P., Pace, S., Panzera, F., and Lombardo, G.: Evaluation of seismic site response in the Maltese archipelago, in: Establishment of an Integrated Italy–Malta Cross–Border System of Civil Protection, Geophysical Aspects, 79–98, Chapter: V, 2015.
Flores Estrella, H., Korn, M., and Alberts, K.: Analysis of the Influence of Wind Turbine Noise on Seismic Recordings at Two Wind Parks in Germany, Journal of Geoscience and Environment Protection, 5, 76–91, https://doi.org/10.4236/gep.2017.55006, 2017.
Friedrich, T., Zieger, T., Forbriger, T., and Ritter, J. R. R.: Locating wind farms by seismic interferometry and migration, J. Seismol., 22, 1469–1483, https://doi.org/10.1007/s1095001897790, 2018.
Gortsas, T. V., Triantafyllidis, T., Chrisopoulos, S., and Polyzos, D.: Numerical modelling of microseismic and infrasound noise radiated by a wind turbine, Soil Dyn. Earthq. Eng., 99, 108–123, https://doi.org/10.1016/j.soildyn.2017.05.001, 2017.
Hayashi, K.: Development of the Surfacewave Methods and Its Application to Site Investigations, PhD thesis, Kyoto University, Japan, 2008.
Karatzetzou, A., Negulescu, C., Manakou, M., François, B., Seyedi, D. M., Pitilakis, D., and Pitilakis, K.: Ambient vibration measurements on monuments in theMedieval City of Rhodes, Greece, B. Earthq. Eng., 13, 331–345, https://doi.org/10.1007/s1051801496492, 2014.
Kumagai, T., Yanagibashi, T., Tsutsumi, A., Konishi, C., and Ueno, K.: Efficient surface wave method for investigation of the seabed, Soils Found., 60, 648–667, https://doi.org/10.1016/j.sandf.2020.04.005, 2020.
Larose, E.: Lunar subsurface investigated from correlation of seismic noise, Geophys. Res. Lett., 32, L16201, https://doi.org/10.1029/2005gl023518, 2005.
Leiber, C.O.: Assessment of safety and risk with a microscopic model of detonation, Elsevier, Amsterdam Boston, ISBN 0444513329, 2003.
Lerbs, N., Zieger, T., Ritter, J., and Korn, M.: Wind turbine induced seismic signals: the largescale SMARTIE1 experiment and a concept to define protection radii for recording stations, Near Surf. Geophys., 18, 467–482, https://doi.org/10.1002/nsg.12109, 2020.
Neuffer, T. and Kremers, S.: How wind turbines affect the performance of seismic monitoring stations and networks, Geophys. J. Int., 211, 1319–1327, https://doi.org/10.1093/gji/ggx370, 2017.
Neuffer, T., Kremers, S., and Fritschen, R.: Characterization of seismic signals induced by the operation of wind turbines in North RhineWestphalia (NRW), Germany, J. Seismol., 23, 1161–1177, https://doi.org/10.1007/s10950019098667, 2019.
Neuffer, T., Kremers, S., Meckbach, P., and Mistler, M.: Characterization of the seismic wave field radiated by a wind turbine, J. Seismol., 25, 825–844, https://doi.org/10.1007/s10950021100036, 2021.
Rahman, M. and Michelitsch, T.: A note on the formula for the Rayleigh wave speed, Wave Motion, 43, 272–276, https://doi.org/10.1016/j.wavemoti.2005.10.002, 2006.
Saccorotti, G., Piccinini, D., Cauchie, L., and Fiori, I.: Seismic Noise by Wind Farms: A Case Study from the Virgo Gravitational Wave Observatory, Italy, B. Seismol. Soc. Am., 101, 568–578, https://doi.org/10.1785/0120100203, 2011.
Siena, L. D., Thomas, C., Waite, G. P., Moran, S. C., and Klemme, S.: Attenuation and scattering tomography of the deep plumbing system of Mount St. Helens, J. Geophys. Res.Sol. Ea., 119, 8223–8238, https://doi.org/10.1002/2014jb011372, 2014.
Stammler, K. and Ceranna, L.: Influence of Wind Turbines on Seismic Records of the Gräfenberg Array, Seismol. Res. Lett., 87, 1075–1081, https://doi.org/10.1785/0220160049, 2016.
Styles, P., Stimpson, I., Toon, S., England, R., and Wright, M.: Microseismic and Infrasound Monitoring of Low Frequency Noise and Vibrations from Windfarms, PhD thesis, Keele University, UK, 2005.
Welch, P.: The use of fast Fourier transform for 335 the estimation of power spectra: A method based on time averaging over short, modified periodograms, IEEE Transactions on Audio and Electroacoustics, 15, 70–73, https://doi.org/10.1109/tau.1967.1161901, 1967.
Westwood, R. F. and Styles, P.: Assessing the seismic wavefield of a wind turbine using polarization analysis, Wind Energy, 20, 1841–1850, https://doi.org/10.1002/we.2124, 2017.
Zieger, T. and Ritter, J. R. R.: Influence of wind turbines on seismic stations in the upper rhine graben, SW Germany, J. Seismol., 22, 105–122, https://doi.org/10.1007/s1095001796949, 2018.
Zieger, T., Nagel, S., Lutzmann, P., Kaufmann, I., Ritter, J., Ummenhofer, T., Knödel, P., and Fischer, P.: Simultaneous identification of wind turbine vibrations by using seismic data, elastic modeling and laser Doppler vibrometry, Wind Energy, 23, 1145–1153, https://doi.org/10.1002/we.2479, 2020.