# Adjusting of Wind Input Source Term in WAVEWATCH III Model for the Middle-Sized Water Body on the Basis of the Field Experiment

Authors: Alexandra Kuznetsova; Georgy Baydakov; Vladislav Papko; Alexander Kandaurov; Maxim Vdovin; Daniil Sergeev; Yuliya Troitskaya
Journal: Advances in Meteorology
 (2016)

---

## Abstract

Adjusting of wind input source term in numerical model WAVEWATCH III for the middle-sized water body is reported. For this purpose, the field experiment on Gorky Reservoir is carried out. Surface waves are measured along with the parameters of the airflow. The measurement of wind speed in close proximity to the water surface is performed. On the basis of the experimental results, the parameterization of the drag coefficient depending on the 10 m wind speed is proposed. This parameterization is used in WAVEWATCH III for the adjusting of the wind input source term within WAM 3 and Tolman and Chalikov parameterizations. The simulation of the surface wind waves within tuned to the conditions of the middle-sized water body WAVEWATCH III is performed using three built-in parameterizations (WAM 3, Tolman and Chalikov, and WAM 4) and adjusted wind input source term parameterizations. Verification of the applicability of the model to the middle-sized reservoir is performed by comparing the simulated data with the results of the field experiment. It is shown that the use of the proposed parameterizationC
D
(
U
10
) improves the agreement in the significant wave height H
S from the field experiment and from the numerical simulation.


---

## Body

## 1. Introduction

Prediction of surface wind waves on the inland water bodies is recognized as an important problem involving many environmental applications, such as safety of the inland navigation and protection from the banks erosion. Lake waves also strongly affect the processes of exchange of momentum, heat, and moisture in the low atmosphere and thus determine microclimate of the adjacent areas, which should be taken into account in planning structure of recreation zones [1].The major physical problem of numerical wave modeling in inland water bodies is associated with short fetches, when parameters of wave excitation and development are strongly different from long-fetch condition typical for the open ocean [2]. Typically, in these conditions, the numerical description of waves in lakes and reservoirs is based on empirical models (see, e.g., [3, 4]). But the empirical relationships are based on the averaged characteristics that cannot predict the extremes important for many tasks of operational meteorology (storm conditions such as Great Lakes storm, discussed in [5]), and numerical wave models are required. Now there are a number of examples of application of third generation models for waves forecast in large lakes. So, WAVEWATCH III [2] is used successfully for the wave forecasts on the Great Lakes in the USA [6, 7]. The data for a current wave situation is presented on the open website and is updated every three hours [8]. Furthermore, WAVEWATCH III and SWAN [9] are applied to Caspian Sea and Ladoga Lake to analyze the wind and waves climate hindcasting [10]. Nevertheless, lakes and reservoirs of smaller sizes (less than 100 km linear size, the so-called middle-sized reservoirs) also have examples of hurricane-force wind and severe surface wave states conditions. The first attempt of application of a global wave model WAM [11] for the wave forecasting on a middle-sized reservoir was reported recently in [12].Among the peculiarities of the low-fetch waves at the middle-sized reservoirs is the stronger wind input, which is proportional to the ratio of wind friction velocity (or 10 m wind speed) to the wave phase velocity [2]. Another feature is the enhanced nonlinearity caused by higher steepness of the waves. Then the tuning of the ocean wave model to the inland water conditions should be twofold: adjusting of the wind source and “collision integral.” Dissipation due to wave breaking can be expected to be similar to wave conditions due to their universal nature.One more problem of tuning of numerical models to the conditions of middle-sized reservoirs and lakes is the small amount of experimental data that can be used for its verification. Rare examples of such experiments are studied [13, 14], which show the specificity of wind-wave interaction in the indicated circumstances. In this paper, we present a tuning of the wind input term in WAVEWATCH III model to the conditions of the middle-sized reservoir on an example of Gorky Reservoir belonging to the Volga Cascade. The tuning is based on the data of the field experiment held by our group. The methods of the experiment are different from those used in [13, 14], and they focus on the study of airflow in close proximity to water surface. The comparison of results of the numerical experiments with the results of the field experiments on Gorky Reservoir is presented.The paper has the following structure. In Section2, the basic parameterizations of wind and wave interactions in WAVEWATCH III v.3.14 are presented. In Section 3, the studied reservoir and the field experiment setup with its instrumentation and methods are described. In Section 4, the field data processing and the results of the field experiment are presented. In Section 5, numerical experiment within tuned to the conditions of the middle-sized reservoir WAVEWATCH III is set out. In Section 6, the adjustment of wind input source term in tuned WAVEWATCH III and the results of numerical experiments within tuned WAVEWATCH III and within a tuned WAVEWATCH III with the adjusted wind input are presented.


## 2. Wind Input Parameterizations in WAVEWATCH III

WAVEWATCH III [2, 15, 16] is based on the numerical solution of the equation for the spectral density of wave action N in the approximation of phase averaging:(1)
∂
N
∂
t
+
∇
x
x
˙
N
+
∂
∂
σ
σ
˙
N
+
∂
∂
θ
θ
˙
N
=
1
σ
S
in
+
S
dis
+
S
nl
.The left hand side of  (1) describes the kinematics of waves, σ is the radian frequency, and θ is the wave direction. In the right hand side, there are terms that describe the wind-wave growth S
in, dissipation mainly due to wave breaking S
dis, and 4-wave nonlinear interaction of waves S
nl.This paper focuses on the wind input parameterizations, whose parameters can be adjusted on the basis of the field measurements data. Generally, the term describing the wind input is determined as(2)
S
in
=
β
k
,
θ
N
k
,
θ
σ
,where β
k
,
θ is the dimensionless wind-wave growth rate parameter approximated in WAVEWATCH III v.3.14 by different parameterizations (WAM 3, Tolman and Chalikov, and WAM 4). Among them, the WAM 3 parameterization [17–19] is defined by two empirical formulas. The first one is for the wind-wave growth rate(3)
β
k
,
θ
=
C
in
ρ
a
ρ
w
max
⁡
0
,
28
u
∗
c
ph
cos
⁡
θ
-
θ
w
-
1
,where C
in = 0.25 is a constant, ρ
a
/
ρ
w is the ratio of the densities of air and water, u
∗ is the friction velocity, c
ph is the phase velocity, and θ
w is the main wind direction. The second one is a relation between the 10 m wind speed U
10 and a friction velocity u
∗
=
U
10
C
D provided by the empirical formula of the parameterization of the surface drag coefficient C
D, proposed in [19]:(4)
C
D
=
0.001
×
0.8
+
0.65
U
10
.In the Tolman and Chalikov parameterization [20], the dimensionless coefficient depends on the drag coefficient and on the dimensionless frequency of the spectral components. In the WAM 4 parameterization [21], the wind-waves interaction parameter β is presented by the adjusted Miles formula. To calculate the roughness parameter, the feedback of the wind-waves spectra is taken into account as well. So, in the considered parameterizations, the wind input is determined by the wind-wave interaction parameter β and by the dependence u
∗
U
10 which is defined by C
D.


## 3. Field Experiment

The tuning of WAVEWATCH III to the conditions of the middle-sized reservoir was performed for conditions of Gorky Reservoir on the basis of field experiment.Measurements were carried out in 2012–2014 from May to October in the waters of Gorky Reservoir. It has an elongated shape (Figure1), which allows studying wind waves of different fetch, depending on the wind direction. Gorky Reservoir is an artificial lake in the central part of the Volga River formed by the dam of Nizhny Novgorod Hydroelectric Station. Its lake part is 85 km long and up to 15 km wide. The mean depth of the reservoir varies in the range of 4–20 meters, and the depth in the area of measurements is 9–12 meters depending on the season and the point of measurement.Figure 1
(a) Gorky Reservoir (Google Earth data). (b) Zoom view of the measurements area.
(a)
(b)


### 3.1. Instrumentation

Instrumentation was placed on a buoy station with the original design based on the oceanographic Froude buoy. Froude buoy is a mast submerged in water and held in a vertical position by the float close to the surface and by the load on the depth (Figures2(a) and 2(b)). Buoy’s total length is 12 m, and the length of the part above the water is 5.3 m. The resonant frequency of the vertical oscillations is 0.25 Hz, which corresponds to a wavelength of 25 m. On the buoy mast, 4 ultrasonic speed sensors (WindSonic Gill Instruments Production) are located at heights of 0.85 m, 1.3 m, 2.27 m, and 5.26 m over the mean water surface. A fifth sensor is located on the float tracking waveform that allows measuring the wind speed in close proximity to the water surface. The distance from the float to the buoy mast is 1 m; the height of the wind speed measuring zone is 10 cm from the water surface. The buoy is also equipped with air temperature sensors (at heights of 0.1 m (float), 0.85 m, and 1.3 m), water temperature sensors, and three-channel wire wave gauge that allows us to retrieve the wave space-time spectra.Figure 2
Froude buoy: (a) real view of the operating state and (b) scheme.
(a)
(b)WindSonic is two-component ultrasonic sensor with 4% measurement accuracy and velocity resolution of 0.01 m/s. Operating range of wind speed measurements 0–60 m/s includes measurements in calm conditions. Resistive temperature sensors measure the environmental temperature with resolution of 0.01°C and 3% measurement accuracy. Wave gauge consists of three pairs of resistive wire sensors, located at the vertices of an equilateral triangle with a side of 62 mm, and the data sampling rate is 100 Hz.


### 3.2. Surface Wave Spectra

Three-dimensional frequency-wavenumber spectra were retrieved from wave gauge data by the algorithm similar to the wavelet directional method (WDM) suggested in [22]. The details of the applied method are described in [23]. Time series of water elevation from each pair of wire sensors were processed by the window FFT with the window width 2
N (N is an integer) with 50% overlapping. The complex amplitudes of harmonics at each frequency ω
A
ω
x
n
,
y
n
exp
⁡
i
ϕ
ω
x
n
,
y
n were calculated; here A
ω is the wave magnitude, ϕ
ω is the wave phase, n
=
1, 2, 3 is the number of the wave sensors, and x
n
,
y
n are the n-sensor’s Cartesian coordinates. Suppose that the wave field is a superposition of harmonic waves with the wavenumbers k
→
=
k
x
,
k
y and one harmonic wave dominates in each interrogation window, and then their phases are(5)
ϕ
ω
x
n
,
y
n
=
k
x
x
n
+
k
y
y
n
.Thus, the wavenumber components can be calculated by the phase difference at different wave sensors. To obtain the directional spectra, the Cartesian coordinatesk
x
,
k
y were transformed to the polar coordinates k
,
θ and then 3D spectrum E
σ
,
k
,
θ was obtained similar to [22] by binning the amplitudes squared into calculated bins in k and θ. Integration of E
σ
,
k
,
θ over wavenumber or frequency yields frequency E
σ
,
θ or wavenumber E
k
,
θ directional spectra, respectively. Integrating over θ gives the 1D frequency and wavenumber spectra correspondingly (examples are in Figure 3). The system allows estimating the parameters of the wave, whose length exceeds the double distance d between the sensors:(6)
k
max
⁡
=
π
d
=
0.5
cm
-
1
.Figure 3
Wave spectra forU
10 = 6-7 m/s averaged over 60 minutes: (a) 1D frequency spectrum and (b) 1D wavenumber spectrum.
(a)
(b)The developed algorithm is based on the supposition that the dominating wave field within the interrogation window at a given frequencyσ is a harmonic wave, which is correct for a rather short time interval due to grouping of the surface wave field. Several periods of the energy-wave should fit into the interrogation window. For typical values of the observed peak frequency f
p = 0.35–0.8 Hz (or the observed peak period T
p = 1.25–2.80 s), the size of interrogation window is selected to be 20.48 s (2048 points).Frequency and wavenumber spectra shown in Figure3 are obtained by processing of two consecutive recordings 60 minutes long. It should be noted that the asymptotic behavior of frequency spectrum and wavenumber spectrum (σ
-
5, k
-
3) correspond to Philips saturation spectrum (see, e.g., [24]), but they do not correspond to Toba spectrum (see, e.g., [25]) which is typical for the ocean conditions.For the comparison with the numerical modeling results, the values of significant wave heightH
S are estimated as four-standard deviation of the water surface elevation, which are calculated as the integral of the frequency spectrum: (7)
H
S
=
4
∫
E
σ
d
σ
1
/
2
.


### 3.3. Parameters of Low Atmosphere over the Inland Water Body

The location of wind velocity sensors corresponds to the structure of the airflow. It is well known (see, e.g., [26]) that the velocity profile in the constant flux layer (where the turbulent momentum flux τ
turb
=
ρ
air
u
x
′
u
z
′ does not depend on the vertical coordinate z; here ρ
air is the air density and u
x
′ and u
z
′ are the fluctuations of horizontal and vertical components of the wind velocity, resp.) has a logarithmic form:(8)
U
z
=
u
∗
κ
ln
⁡
z
z
0
,where u
∗
=
τ
turb
/
ρ
air
=
u
x
′
u
z
′ is the friction velocity determined by the turbulent momentum flux, z
0 is the surface roughness parameter, and κ is the von Karman constant. In [27], a distortion of the velocity field in the presence of a rough surface is considered, and it is shown that, in case of a monochromatic wave propagating along the wind, for the stream function φ averaged over the turbulent fluctuations, an equation can be written as(9)
U
-
c
d
2
φ
d
η
2
-
k
2
φ
-
d
2
U
d
η
2
φ
=
0
,where U
(
η
) is the dependence of wind speed on height above the surface (in curvilinear coordinates), c is phase velocity of the wave, and k is wavenumber. When the magnitude U
′
′
/
k
2
U
-
c is much greater or much less than 1, an approximate solution is the function φ
=
A
U
-
c
e
-
k
η, where A is wave amplitude. In the case of the logarithmic velocity profile, this condition takes the form (see [27])(10)
u
∗
/
κ
k
η
2
U
-
c
≪
1
o
r
u
∗
/
κ
k
η
2
U
-
c
≫
1
,and it is performed well enough at the height of the order of the wave amplitude and higher. Thus, the perturbation is the bending of the flow lines along the rough surface, and it decreases exponentially with altitude. Therefore, to fix speed sensor relatively to the mean streamlines, the speed at a distance from surface must be measured at a fixed horizon, and the measurements close to the surface should be carried out using a tracking waveform sensor located on the float.It is important that the lower sensor should not be located in the wave boundary layer. The magnitude of the wave boundary layerε can be estimated in accordance with [28](11)
k
~
u
~
∗
κ
ln
⁡
z
z
0
-
c
~
z
=
ε
=
κ
u
~
∗
z
ε
2
z
=
ε
,where k
~, c
~, and u
~
∗ denote typical values of the mentioned above variables. In the conditions of the Gorky Reservoir (k
~ = 2-3 m−1, u
~
∗ = 0.1–0.4 m/s), the evaluation of the height of the wave boundary layer gives a value ε
~
1 mm that is significantly less than the height of the lower sensor measuring wind speed location.It should be noted that the buoy is connected to the vessel by the cable and is located at a fixed distance of 30 meters; the cross section of the vessel is approximately equal to 3 × 3 m. Thus, the buoy is located at a distance of about 10 characteristic dimensions from the obstacle. According to the recommendations of Gill Instruments [29], it is enough to consider the airflow to be unperturbed, but the measurements show the presence of small (3–6%) deviations of the measured profile from the logarithmic form even in a steady wind and neutral stratification (Figure 4). This deviation consists of the speed increasing at a height of 5.26 m on 3–6% compared with the logarithmic approximation. Apparently, this is due to the effect of the vessel shielding on the 4 lower sensors.Figure 4
Examples of the deviation of the wind velocity profile (5 minutes averaging) from the logarithmic form under different conditions.To estimate the perturbations introduced by vessel, the formulas of the defect rate attenuation in the turbulent wake are used. It is known [30] that the perturbation introduced by the body in the unlimited turbulent flow decreases in proportion x
/
D
-
2
/
3 along the flow and in proportion exp
⁡
-
r
2
/
2
σ
~
2 along the cross section. Here, D is the characteristic size of the body, x and r are the cylindrical coordinates, σ
~ is a width of the Gaussian function, and(12)
σ
~
D
~
x
D
1
/
3
;the coefficients of the proportionality are determined empirically. In the case of airflow vessel shielding in a semibounded space, it is fair to expect that the nature of the turbulent wake is the same, but the appropriate factors must be chosen. Since the top speed sensor is located outside the vessel shielding zone, and the lower four are well approximated by a logarithmic dependence, the value σ
~ should be between the heights of the fourth and fifth sensors. The linear coefficient is selected so that the velocity profiles are well approximated by a logarithmic dependence in average. Finally, to determine the magnitude of the perturbation introduced by the vessel in the airflow, the dependence (13)
γ
=
U
′
U
0
=
0.3
x
D
-
2
/
3
exp
⁡
-
r
2
2
σ
~
2is used, where σ
~
/
D
=
0.4
x
/
D
1
/
3, U
0 is the estimated wind speed in the unperturbed flow, and U
′ is the magnitude of the perturbation. Note that the selected coefficients are similar to those obtained in [31] for the wake behind the sphere. To determine the wind velocity in the unperturbed flow, the wind speed measured by sensors is multiplied by the coefficient α
=
1
-
γ
-
1. The coefficients α for the various sensors are equal to 1.069, 1.065, 1.061, 1.047, and 1.009, respectively, from the bottom up to the top.


## 3.1. Instrumentation

Instrumentation was placed on a buoy station with the original design based on the oceanographic Froude buoy. Froude buoy is a mast submerged in water and held in a vertical position by the float close to the surface and by the load on the depth (Figures2(a) and 2(b)). Buoy’s total length is 12 m, and the length of the part above the water is 5.3 m. The resonant frequency of the vertical oscillations is 0.25 Hz, which corresponds to a wavelength of 25 m. On the buoy mast, 4 ultrasonic speed sensors (WindSonic Gill Instruments Production) are located at heights of 0.85 m, 1.3 m, 2.27 m, and 5.26 m over the mean water surface. A fifth sensor is located on the float tracking waveform that allows measuring the wind speed in close proximity to the water surface. The distance from the float to the buoy mast is 1 m; the height of the wind speed measuring zone is 10 cm from the water surface. The buoy is also equipped with air temperature sensors (at heights of 0.1 m (float), 0.85 m, and 1.3 m), water temperature sensors, and three-channel wire wave gauge that allows us to retrieve the wave space-time spectra.Figure 2
Froude buoy: (a) real view of the operating state and (b) scheme.
(a)
(b)WindSonic is two-component ultrasonic sensor with 4% measurement accuracy and velocity resolution of 0.01 m/s. Operating range of wind speed measurements 0–60 m/s includes measurements in calm conditions. Resistive temperature sensors measure the environmental temperature with resolution of 0.01°C and 3% measurement accuracy. Wave gauge consists of three pairs of resistive wire sensors, located at the vertices of an equilateral triangle with a side of 62 mm, and the data sampling rate is 100 Hz.


## 3.2. Surface Wave Spectra

Three-dimensional frequency-wavenumber spectra were retrieved from wave gauge data by the algorithm similar to the wavelet directional method (WDM) suggested in [22]. The details of the applied method are described in [23]. Time series of water elevation from each pair of wire sensors were processed by the window FFT with the window width 2
N (N is an integer) with 50% overlapping. The complex amplitudes of harmonics at each frequency ω
A
ω
x
n
,
y
n
exp
⁡
i
ϕ
ω
x
n
,
y
n were calculated; here A
ω is the wave magnitude, ϕ
ω is the wave phase, n
=
1, 2, 3 is the number of the wave sensors, and x
n
,
y
n are the n-sensor’s Cartesian coordinates. Suppose that the wave field is a superposition of harmonic waves with the wavenumbers k
→
=
k
x
,
k
y and one harmonic wave dominates in each interrogation window, and then their phases are(5)
ϕ
ω
x
n
,
y
n
=
k
x
x
n
+
k
y
y
n
.Thus, the wavenumber components can be calculated by the phase difference at different wave sensors. To obtain the directional spectra, the Cartesian coordinatesk
x
,
k
y were transformed to the polar coordinates k
,
θ and then 3D spectrum E
σ
,
k
,
θ was obtained similar to [22] by binning the amplitudes squared into calculated bins in k and θ. Integration of E
σ
,
k
,
θ over wavenumber or frequency yields frequency E
σ
,
θ or wavenumber E
k
,
θ directional spectra, respectively. Integrating over θ gives the 1D frequency and wavenumber spectra correspondingly (examples are in Figure 3). The system allows estimating the parameters of the wave, whose length exceeds the double distance d between the sensors:(6)
k
max
⁡
=
π
d
=
0.5
cm
-
1
.Figure 3
Wave spectra forU
10 = 6-7 m/s averaged over 60 minutes: (a) 1D frequency spectrum and (b) 1D wavenumber spectrum.
(a)
(b)The developed algorithm is based on the supposition that the dominating wave field within the interrogation window at a given frequencyσ is a harmonic wave, which is correct for a rather short time interval due to grouping of the surface wave field. Several periods of the energy-wave should fit into the interrogation window. For typical values of the observed peak frequency f
p = 0.35–0.8 Hz (or the observed peak period T
p = 1.25–2.80 s), the size of interrogation window is selected to be 20.48 s (2048 points).Frequency and wavenumber spectra shown in Figure3 are obtained by processing of two consecutive recordings 60 minutes long. It should be noted that the asymptotic behavior of frequency spectrum and wavenumber spectrum (σ
-
5, k
-
3) correspond to Philips saturation spectrum (see, e.g., [24]), but they do not correspond to Toba spectrum (see, e.g., [25]) which is typical for the ocean conditions.For the comparison with the numerical modeling results, the values of significant wave heightH
S are estimated as four-standard deviation of the water surface elevation, which are calculated as the integral of the frequency spectrum: (7)
H
S
=
4
∫
E
σ
d
σ
1
/
2
.


## 3.3. Parameters of Low Atmosphere over the Inland Water Body

The location of wind velocity sensors corresponds to the structure of the airflow. It is well known (see, e.g., [26]) that the velocity profile in the constant flux layer (where the turbulent momentum flux τ
turb
=
ρ
air
u
x
′
u
z
′ does not depend on the vertical coordinate z; here ρ
air is the air density and u
x
′ and u
z
′ are the fluctuations of horizontal and vertical components of the wind velocity, resp.) has a logarithmic form:(8)
U
z
=
u
∗
κ
ln
⁡
z
z
0
,where u
∗
=
τ
turb
/
ρ
air
=
u
x
′
u
z
′ is the friction velocity determined by the turbulent momentum flux, z
0 is the surface roughness parameter, and κ is the von Karman constant. In [27], a distortion of the velocity field in the presence of a rough surface is considered, and it is shown that, in case of a monochromatic wave propagating along the wind, for the stream function φ averaged over the turbulent fluctuations, an equation can be written as(9)
U
-
c
d
2
φ
d
η
2
-
k
2
φ
-
d
2
U
d
η
2
φ
=
0
,where U
(
η
) is the dependence of wind speed on height above the surface (in curvilinear coordinates), c is phase velocity of the wave, and k is wavenumber. When the magnitude U
′
′
/
k
2
U
-
c is much greater or much less than 1, an approximate solution is the function φ
=
A
U
-
c
e
-
k
η, where A is wave amplitude. In the case of the logarithmic velocity profile, this condition takes the form (see [27])(10)
u
∗
/
κ
k
η
2
U
-
c
≪
1
o
r
u
∗
/
κ
k
η
2
U
-
c
≫
1
,and it is performed well enough at the height of the order of the wave amplitude and higher. Thus, the perturbation is the bending of the flow lines along the rough surface, and it decreases exponentially with altitude. Therefore, to fix speed sensor relatively to the mean streamlines, the speed at a distance from surface must be measured at a fixed horizon, and the measurements close to the surface should be carried out using a tracking waveform sensor located on the float.It is important that the lower sensor should not be located in the wave boundary layer. The magnitude of the wave boundary layerε can be estimated in accordance with [28](11)
k
~
u
~
∗
κ
ln
⁡
z
z
0
-
c
~
z
=
ε
=
κ
u
~
∗
z
ε
2
z
=
ε
,where k
~, c
~, and u
~
∗ denote typical values of the mentioned above variables. In the conditions of the Gorky Reservoir (k
~ = 2-3 m−1, u
~
∗ = 0.1–0.4 m/s), the evaluation of the height of the wave boundary layer gives a value ε
~
1 mm that is significantly less than the height of the lower sensor measuring wind speed location.It should be noted that the buoy is connected to the vessel by the cable and is located at a fixed distance of 30 meters; the cross section of the vessel is approximately equal to 3 × 3 m. Thus, the buoy is located at a distance of about 10 characteristic dimensions from the obstacle. According to the recommendations of Gill Instruments [29], it is enough to consider the airflow to be unperturbed, but the measurements show the presence of small (3–6%) deviations of the measured profile from the logarithmic form even in a steady wind and neutral stratification (Figure 4). This deviation consists of the speed increasing at a height of 5.26 m on 3–6% compared with the logarithmic approximation. Apparently, this is due to the effect of the vessel shielding on the 4 lower sensors.Figure 4
Examples of the deviation of the wind velocity profile (5 minutes averaging) from the logarithmic form under different conditions.To estimate the perturbations introduced by vessel, the formulas of the defect rate attenuation in the turbulent wake are used. It is known [30] that the perturbation introduced by the body in the unlimited turbulent flow decreases in proportion x
/
D
-
2
/
3 along the flow and in proportion exp
⁡
-
r
2
/
2
σ
~
2 along the cross section. Here, D is the characteristic size of the body, x and r are the cylindrical coordinates, σ
~ is a width of the Gaussian function, and(12)
σ
~
D
~
x
D
1
/
3
;the coefficients of the proportionality are determined empirically. In the case of airflow vessel shielding in a semibounded space, it is fair to expect that the nature of the turbulent wake is the same, but the appropriate factors must be chosen. Since the top speed sensor is located outside the vessel shielding zone, and the lower four are well approximated by a logarithmic dependence, the value σ
~ should be between the heights of the fourth and fifth sensors. The linear coefficient is selected so that the velocity profiles are well approximated by a logarithmic dependence in average. Finally, to determine the magnitude of the perturbation introduced by the vessel in the airflow, the dependence (13)
γ
=
U
′
U
0
=
0.3
x
D
-
2
/
3
exp
⁡
-
r
2
2
σ
~
2is used, where σ
~
/
D
=
0.4
x
/
D
1
/
3, U
0 is the estimated wind speed in the unperturbed flow, and U
′ is the magnitude of the perturbation. Note that the selected coefficients are similar to those obtained in [31] for the wake behind the sphere. To determine the wind velocity in the unperturbed flow, the wind speed measured by sensors is multiplied by the coefficient α
=
1
-
γ
-
1. The coefficients α for the various sensors are equal to 1.069, 1.065, 1.061, 1.047, and 1.009, respectively, from the bottom up to the top.


## 4. Field Data Processing and Results

Investigation of the wind flow parameters is carried out by profiling. General recording of wind speed is up to 5-hour duration and is divided into 5-minute sections (300 measuring points) with a 50% of overlap. As a result of the averaging, 5 values of wind speed corresponding to five horizons of measurements are received for each section. The resulting averaged profile is approximated by function (8) with the parameters of the approximation u
∗ (friction velocity) and z
0 (surface roughness). The values of the 10 m wind speed U
10 and of the drag coefficient C
D are retrieved from the resulting approximation.The impact of data obtained from different horizons on the resulting approximation of the wind velocity profile is analyzed. Figure5 shows a comparison of the retrieved dependence C
D
U
10 for two combinations of speed sensors: with and without the lower sensor. The points on the plot are the result of the binning of the wind speed data within a cell with size of Δ
U
10 = 0.5 m/s. The error gates correspond to the standard deviation. Also the results of field experiments [13, 14] and empiric oceanic parameterization COARE 3.0 [32] are shown in Figure 5. It is evident that, without the lower sensor data, the values C
D
U
10 are characterized by a slightly greater spread and are located higher and more close to the results of [13, 14, 32], while taking the lower sensor into account shows lower values of the drag coefficient. In [13], the sensors are located at the heights of 0.89 m up to 10 m from the mean water level; in [14], the sensors are located at heights from 0.5 m up to 4 m. In both cases, all the sensors are firmly fixed on the masts, and the approximation is performed throughout the wind speed profile.Figure 5
Comparison of the retrieved dependenceC
D
U
10 with and without the lower sensor: gray diamonds denote the binning of data with the lower sensor (the standard deviation as the error gates); black solid circles are the binning of the data without the lower sensor (the standard deviation as the error gates); gray circles are the results of field experiment [13]; gray crosses are the results of field experiment [14]; dashed line is the empiric oceanic parameterization COARE 3.0.Figure6 shows a comparison of the retrieved dependence C
D
U
10 using two lower sensors only and using all five sensors. The use of two sensors only reveals significant differences in the wind parameters recovery in the weak winds area: the scatter and absolute value of C
D
U
10 decreases. In the field of moderate and strong winds, the retrieved values of drag coefficient differ slightly despite the small increase in the measurement error.Figure 6
Comparison of retrieved dependenceC
D
U
10 using two and five sensors: gray solid circles are the binning of the five-sensor data (the standard deviation as the error gates); black diamonds are the binning of the two-sensor data (the standard deviation as the error gates); gray circles are the results of the field experiment [13]; gray crosses are the results of the field experiment [14]; dashed line is the empiric oceanic parameterization COARE 3.0.These results can be explained by the deviation of wind velocity profile from the logarithmic form. It may be caused by the stratification of the surface layer of the atmosphere and impact of the coastline as well as gustiness of the wind, because the lower part of the profile adapts more quickly to the changing conditions of waves, and the airflow parameters determine the momentum transfer from the wind to the waves exactly at the water-air boundary. Consequently, further analysis of the dependenceC
D
U
10 was based on the measured data from sensors 1 and 2.Throughout the 2012–2014 years, the dataset consisting of approximately 100 hours of recordings in the range of wind velocities 1–12 m/s for different fetch values (L
fetch = 1–50 km) and stratification conditions (T
air
-
T
water
=
-
5
-
1
5
°
C) is obtained. The resulting dependence C
D
U
10 using two lower sensors is received and shown in Figure 7. The approximation of the obtained data is made by a function(14)
C
D
=
0.00124
U
10
-
1
+
0.00034
+
0.000049
U
10
.Figure 7
Detailed retrieved dependenceC
D
U
10: gray diamonds denote the dependence received using two lower sensors; solid line is an approximation of the obtained data by a function C
D
=
0.00124
·
U
10
-
1
+
0.00034
+
0.000049
·
U
10; black circles denote binning of the two-sensor data (the standard deviation as the error gates); dashed line is the empiric oceanic parameterization COARE 3.0.The result of binning of the wind speed data within a cell with sizeΔ
U
10 = 0.5 m/s and the standard deviation as the error gates are shown in Figure 7. The dependence (14) is used below to adjust the wind input in WAVEWATCH III.


## 5. Tuning of WAVEWATCH III: Numerical Experiment

The tuning of WAVEWATCH III to the conditions of the Gorky Reservoir consists of the following factors. In the open code, the minimum value of a significant wave heightH
S is adjusted in a number of blocks, where the lowest value of H
S is directly indicated. For a description of the reservoir, the topographic grid of the Gorky Reservoir with dimensions 72 × 108 and increments of 0.00833° (which corresponds to approximately 800 m by 900 m for the considered latitudes) is used. The grid is taken from the NOAA data “Global Land One-Kilometer Base Elevation (GLOBE).” Topographic grid of Gorky Reservoir is shown in Figure 8. There is no reliable information about the bathymetry of the considered area, although the navigational maps show that the depth is big enough to consider an approximation of deep water. Thus, the constant depth of 9 m is taken.Figure 8
Topographical grid of Gorky Reservoir. Computational cell with a size of 0.00833° is shown.The frequency range is changed to 0.2–4 Hz in accordance with the experimentally observed range, which is split in 31 frequencies in the simulation and is modeled by a logarithmic formula for the frequency growth(15)
σ
N
=
δ
N
-
1
σ
1
,where the growth rate is determined to be δ
=
1.1 in accordance with the recommendations of [2]; 30 angular directions of the wave field are considered. The initial seeding is triggered and evolved in the wind. In practice, to simulate wind waves on the surface of the seas and oceans, the reanalysis data is typically used as a wind forcing. In the middle-sized inland waters, this approach is not applicable because of its too low spatial resolution (2.5°). In addition, in this area there are only two weather stations (Volga GMO, Yuryevets), but they are on the coast, and it was found out that the wind speed on the coast is different from those over the waters of the reservoir. In this regard, the magnitude and direction of the wind data are taken from the field experiment and are considered to be homogeneous over the whole water area of the reservoir. In fact, the wind field is expected to be heterogeneous, as such factors like the elongated shape of the reservoir and the high banks can lead to a significant spatial variability of the wind field. It should be noted that this assumption of the homogeneity of the wind forcing over the pond can be a source of errors in the numerical experiment.As the wind above the reservoir is characterized by a strong mutability, the averaging of the wind speed in the experiment is performed in the interval of 15 minutes. Thus, the simulation is held with input data updated every 15 minutes, measured in field experiment: 10 m wind speed and direction, the water-air temperature difference. The comparison is made for the following output: 1D spectra elevations, significant wave height, and the average wave period. All data is obtained at the point corresponding to the point of observations and is averaged in the range of 15 minutes to match the similarity with the averaged data of the field experiment.


## 6. Adjusting of Wind Input Source Term in Tuned WAVEWATCH III

For the further tuning of WAVEWATCH III, a comparison of the used parameterizations of wind input is performed. For this, coefficients determining WAVEWATCH III parameterizations of wind input (C
D and β) are displayed at each step of the numerical simulation. Figure 9(a) shows the dependence β
(
σ
u
∗
/
g
) for three considered parameterizations, results of numerical simulations for Tolman and Chalikov parameterizations, and theoretically estimated curves for parameterizations WAM 3 and WAM 4. Figure 9(b) shows the dependence C
D
U
10 of WAVEWATCH III and the proposed parameterization C
D
U
10 obtained as a result of field measurements (14). It can be seen that different parameterization of β and C
D used in WAVEWATCH III v.3.14 gives similar values in the conditions of middle-sized reservoirs (σ
u
∗
/
g = 0.1–0.15, U
10 = 1–10 m/s). At the same time, the proposed parameterization of C
D gives significantly lower values (on ≈50%) for moderate and strong wind speeds U
10 > 4 m/s.Figure 9
(a) Dependence of the wind-wave interaction coefficient on the dimensionless friction velocity for different parameterizations. (b) Dependence of the surface drag coefficientC
D on the wind speed U
10. WAM 3 is indicated by a dashed line, Tolman and Chalikov by circles, and WAM 4 by crosses.
(a)
(b)Then, wind-waves regime is studied for built-in parameterizations WAM 3, Tolman and Chalikov, and WAM 4. In the tuned WAVEWATCH III, we also used the adjusted wind input parameterization consisting of the use of the proposed experimental parameterizationC
D within parameterizations WAM 3 and Tolman and Chalikov. As it was shown in Section 2, WAM 3 is based on the explicit formula, and Tolman and Chalikov is based on the implicit formula of C
D
U
10. These formulas are modified directly in the program code. Thus, instead of built-in parameterizations of C
D in the model, the new proposed parameterization (14) is used. Figures 10(a) and 10(b) show the results of the modeling and field measurements for the days 13.06.13 and 20.06.14, which are typical for the biggest part of the considered data. The lower plots show the measured values of the wind used in the simulation, and the plots on the top show a change for the retrieved values of H
S, obtained both from the field experiment and from the numerical experiments. In the model calculation, H
S is based on a formula(16)
H
S
=
4
E
.Figure 10
The upper graphs: dependence ofH
S on time. Results of field experiments are marked by crosses, the simulated values of H
S for parameterizations WAM 3 (dark gray dashed line), WAM 3 with new C
D (dark gray solid line), Tolman and Chalikov (black dashed line), Tolman and Chalikov with new C
D (black solid line), and WAM 4 (long dotted line). The lower graphs: evolution of wind given by the field experiment, specified as input for WAVEWATCH III (a) 13.06.13, (b) 20.06.14.
(a)
(b)The same value gives the calculation ofH
S in the experiment (7). As it can be seen from Figure 10, usually the values of the significant wave height in simulations with built-in parameterizations are overestimated. But it can be seen that the use of the proposed parameterization C
D
U
10 improves the agreement with the field experiment. The dependence H
S
(
t
) (Figure 10) shows that at the beginning of the time interval of the measurements (first 50 minutes) the wave regime is developing only, while the experimental values are already much higher. First of all, this is due to the fact that the wind prehistory for the simulations (before the start of the measurements) for both dates is taken from the weather stations, and, as it is mentioned in Section 5, the wind speed on the coast is different from those over the waters of the reservoir. We also associate the difference between the simulation output and experimental values with the inaccuracy of wind forcing: due to the fact that the wind is set to be homogeneous, the waves that come from other parts of the reservoir are not large enough. To better match the results, inhomogeneous wind field is required.Table1 describes the evaluation of the difference in the applying of different parameterizations of wind input source term for two test dates: 13.06.13 and 20.06.14. It can be seen that WAM 3 typically overestimates the values of H
S compared with the experimental data for both dates, whereas the use of the proposed new C
D improves the accordance very well.Table 1
Standard deviation ofH
S.
WAM 3 + newC
D
WAM 3
Tolman and Chalikov  + newC
D
Tolman and Chalikov
WAM 4
13.06.13
0.28
0.36
0.19
0.23
0.32
20.06.14
0.23
0.41
0.25
0.20
0.20At the same time, Tolman and Chalikov parameterization underestimates the part of the values for 13.06.13, but the use of newC
D improves the total standard deviation (STD) for all the day of 13.06.13. For 20.06.14, the use of new C
D makes the underestimation of Tolman and Chalikov in the beginning of the time interval bigger, because it decreases the value of the energy entering the system (the new C
D lies lower than C
D in Tolman and Chalikov as it is shown in Figure 9(b)). This worsens the STD for Tolman and Chalikov with new C
D for 20.06.14. It should be mentioned that built-in Tolman and Chalikov source term performs well enough for the conditions of the middle-sized water body. It may be due to the particular properties of the dependence of C
D on U
10 in Tolman and Chalikov parameterization (see Figure 9(b)). Tolman and Chalikov parameterization underestimates the growth rate at low wind speeds (U
10 < 3 m/s) and overestimates it at higher wind speeds (U
10 > 4 m/s) in comparison with the proposed new C
D in the adjusted WAM 3. Possibly, these effects compensate each other and give close results for the integral value of H
S. Nevertheless, new C
D is applied in the framework of this parameterization and shows different results for the dates of 13.06.13 and 20.06.14 (Table 1), but for the analysis of all the measured data with the simulation data it improves the coincidence (as shown in Figure 12(b)).The resulting values of WAM 4 application are usually located between the results of WAM 3 and Tolman and Chalikov, and here in both cases the same situation is realized (the analysis of all data is shown in Figure12(c)).Figure11 shows a comparison of the wave spectra at the point of measurements at a fixed time with the spectra obtained from a numerical experiment for built-in parameterizations WAM 3, Tolman and Chalikov, and WAM 4 and for the adjusted parameterizations WAM 3, Tolman and Chalikov with the new dependence C
D
U
10 as in formula (14). This improvement in the prediction of wave spectra is observed for the new parameterization C
D
U
10 for the biggest part of the considered data (Figure 11(a)). The spectra in the beginning of the time interval that correspond to the situation in Figure 10, where the experimental values exceed the estimated values of H
S, are shown in Figure 11(b). However, along with improving the prediction of wind-wave characteristics, there is still a situation in which the simulated values of H
S are overestimated compared with the experimental values; this situation is reflected in the spectra in Figure 11(c).Figure 11
1D wave spectra. The experimental spectrum is indicated by light-gray bold solid wide line, the simulated spectrum with the parameterizations WAM 3, gray dashed line, WAM 3 with the newC
D, gray solid line, Tolman and Chalikov, dark gray dashed line, Tolman and Chalikov with the new C
D, dark gray solid line, and WAM 4, a long dotted line. (a) The improved prediction of the wave spectra with the use of the new parameterization of C
D; (b) spectra of the beginning of the time interval, where the spectra from the field experiment are higher than the spectra from the numerical experiment; (c) the spectra from the numerical experiment are higher than the spectra from the field experiment.
(a)
(b)
(c)Figure 12
H
S (top graph) and T
m (bottom graph) in comparison with the data of field experiment for (a) parameterizations WAM 3 (diamonds) and WAM 3 with the new C
D (crosses), (b) parameterizations Tolman and Chalikov (diamonds) and Tolman and Chalikov with the new C
D (crosses), and (c) parameterization WAM 4 (diamonds).
(a)
(b)
(c)For all of the considered experiments, a comparison of the integral characteristics (H
S and the mean wave period T
m) is performed. The mean wave period T
m is simulated using the following formula: (17)
T
m
=
T
m
0
,
-
1
=
∫
f
min
⁡
f
r
E
f
d
f
-
1
∫
f
min
⁡
f
r
E
f
f
-
1
d
f
.On Figure12, the x-axis represents the values obtained in the field experiments, and y-axis represents the results of the numerical simulations. On the top plots in Figure 12, for all the considered parameterizations (both built-in and with the use of a new parameterization), values of H
S in the output of the numerical simulation are compared with those obtained from the experiment. The lower plots are for the values of T
m. For all built-in parameterizations, the overestimation of the significant wave height and the underestimation of the mean wave period are typical, and the STD of H
S for WAM 3 is 52%, for Tolman and Chalikov is 40%, and for WAM 4 is 46%. The use of the new parameterization reduces the STD of H
S for WAM 3 from 52% to 39% and for Tolman and Chalikov from 40% to 37%. This is an expected result, as in the numerical experiment with the use of new parameterization of C
D, the wave growth increment is defined more precisely that means that the amount of energy entering the system is simulated more accurately.However, the lower graphs in Figure12 show that the prediction of mean wave periods has significant discrepancy with the measured ones, and the use of the new parameterization of C
D
U
10 does not make sufficient changes. Perhaps this is due to the fact that the adaptation of WAVEWATCH III to marine environment is reflected not only in the function of the wind input, but also in taking into account the specific parameters of numerical nonlinear scheme DIA [33, 34], because nonlinear processes are responsible for the redistribution of the energy received from the wind in the spectrum. WAVEWATCH III considers the wave characteristic of marine and ocean conditions, which have a lower slope compared to the waves on the middle-sized inland waters. The coefficients of proportionality in the scheme DIA are adjusted to the sea conditions. Steeper waves of middle-sized reservoir may require a different adjustment of parameters corresponding to a situation with stronger nonlinearity, which should lead to more rapid frequencies downshift. Consequently, mean wave periods will decrease. At the same time, we can expect that such a tuning of the numerical nonlinear scheme should not affect the quality of the predictions of H
S, which indicates the amount of energy received by the system, but should lead to a better prediction of mean wave periods. This hypothesis will be tested in the subsequent numerical experiments.


## 7. Conclusions

The paper shows the tuning of WAVEWATCH III to the conditions of the middle-sized reservoir on the example of the Gorky Reservoir, which is specified in the model using real topographic grid NOAA “GLOBE.” In carrying out the calculations, the default values of model parameters are modified on the basis of field measurements on the reservoir. In particular, the minimum value of significant wave height is adjusted and frequency range changed to 0.2–4 Hz. The initial seeding is developing under the influence of unsteady uniform wind given by the experiment. The wave field is simulated using both the built-in parameterizations of wind input adapted to the conditions of the open ocean and the parameterization using the new form of the surface drag coefficient (14), which is obtained from a series of field experiments.Field experiments in the Gorky Reservoir show that the values ofC
D in moderate and strong winds are on ≈50% lower than those typical for the ocean conditions. In the course of the experiment, wave characteristics (frequency and wavenumber spectra, the mean wave period, and significant wave height) were obtained for different wind conditions. It is found out that the spectra have the asymptotic behavior similar to Phillips saturation spectrum. Field experiments in the Gorky Reservoir show that the values of C
D in moderate and strong winds are on ≈50% lower than those typical for the ocean conditions.The results of the numerical experiments are compared with the results obtained in the field experiments on the Gorky Reservoir. The use of the built-in parameterizations shows a significant overestimation of the simulatedH
S compared with the experimental results. We interpret it by the overestimation of the turbulent wind stress (friction velocity u
∗) and, accordingly, of the wind input. The use of the new parameterization C
D
(
U
10
), based on field measurements on the reservoir, reduces the values of u
∗ and hence wind increment of surface waves. That improves the agreement in H
S from the field experiment and from the numerical simulation. Comparison of the simulation with built-in ocean parameterizations of the wind input overestimates values of the mean wave period T
m as well. At the same time, the adjustment of the wind input does not affect significantly the agreement of T
m values in the results of numerical simulation and in the field experiment. We interpret this by the fact that the nonlinear scheme is also adjusted to the conditions of seas and oceans, and we plan to make the adjustment of the parameters of the numerical nonlinear scheme DIA to the conditions of the middle-sized reservoir.Another source of possible errors of numerical experiment should also be noted. Due to the lack of sufficient experimental data, the wind speed is assumed to be uniform over the entire water area of the reservoir with the temporal variability defined from the experiment. In fact, nonuniform distribution of the wind is expected, as factors such as the elongated shape of the reservoir and the high banks can lead to a significant spatial variability of the scale with order of 1 km or less. The use of the wind from the reanalysis data is also impossible because of too low spatial resolution (2.5°). Taking into account the high spatial variability is a challenging problem. To solve it, it is planned to use the atmospheric models of high and ultrahigh spatial resolution (e.g., atmospheric model Weather Research & Forecasting (WRF) with block LES (Large Eddy Simulation)).



---

*Source: 8539127-2016-01-04.xml*