Spatial variability of long-term trends in significant wave height over the Gulf of Gdańsk using System Identification techniques

Since the proposed approach to simulation of significant wave height fields depends largely on the nature and shortcomings of the available data, these datasets are described in the following section and listed in Table 1.

As wind velocity fields of quality comparable with outputs of regional numerical weather forecast models spanning the entire last century are nonexistent, global wind velocity fields available from the 20^th Century Reanalysis Project 2 (20CRv2) at NOAA/OAR/ESRL PSD were used (Compo et al. 2011). These data were generated using a global, coupled atmosphere-land model with surface pressure assimilation from the International Surface Pressure Databank, monthly sea surface temperature and ice concentration from the Hadley Centre (Compo et al. 2011). The available 3 h mean wind velocity fields are gridded with approx. 2° (approx. 120 NM) spatial resolution and span the period of 1871–2013 as presented in Table 1. However, due to the nature of the data, grid size and time resolution, the 20CRv2 global wind data are not suitable as a direct input into regional or local wave spectral models, which was demonstrated by smoothed out storm maxima and changed wind directions, when compared to the UMPL-ICM wind data. To correct the most persistent errors, m-ARX data-based models were used to identify the link between the 20CRv2 wind velocities and the high quality UMPL-ICM winds as denoted on the diagram in Figure 1.

The UMPL-ICM wind fields for the 1998–2001 period were generated at the ICM (Interdisciplinary Centre for Mathematical Computing, Warsaw University) using the UMPL (Unified Model for Poland) model, a derivative of the UM (Unified Model) developed by the ECMWF (Herman-Iżycki et al. 2002). The UMPL-ICM wind velocity fields have a relatively high temporal and spatial output resolution of 3 h and 9 NM.

Even after applying this simplified “scaling down” procedure, the quality of simulated wind velocities was not high enough to apply the WAM model directly. The final use of data-driven models, as represented by the lower row of the diagram in Figure 1, was supposed to artificially compensate for the most relevant deficiencies, while linking still deficient wind data with high quality wave height data.

The significant wave height fields generated mostly by Cieślikiewicz & Herman (2002) with the third generation spectral wave model WAM (WAMDI 1988; Komen et al. 1994) were used as the output data to train and verify parametric models. The WAM model was run on a nested grid with 1 NM resolution and 3 h model output, shallow options enabled (frequency table, depth refraction, dissipation including bottom friction, specific nonlinear interactions), and high-resolution (9 NM, 3 h) UMPL-ICM winds (Herman-Iżycki et al. 2002) as the driving input. The WAM dataset consists of model output for 29 storms which happened in 1998–2001.

Grids for WAM, UMPL-ICM and 20CRv2 datasets are shown in Figure 2. The striking contrast between the 20CRv2 global grid (approx. 120 NM) and the local 1 NM WAM grid is easily observed.

Figure 2

Models described in this section are used either to simulate the wind velocity components (m-ARX) or, using the simulated wind velocity, to model the significant wave height (TSK, optimized BH), according to the diagram in Figure 1. They are presented here using the generic notation for the multidimensional input u(t) and the single output y(t) to show the similarities between the TSK and m-ARX as well as due to the fact that m-ARX is used to simulate both velocity components, which would call for a generic notation anyway.

One of the most common ways to link the I-dimensional input u(t) and single output y(t) data under N different regimes is to associate a linear model with each of the regimes. These n = 1,2,…,N regimes can be described in numerous ways: providing their characteristic values ^cu_n, subsets or ranges of values in the input vector space. Then, the input data vector u(t) can be divided along with the corresponding output values into N regimes to provide enough data to estimate N linear models. The resulting output takes the form of an interpolation between N linear models:

Y(t)=∑<!-- ? -->n=1NWn(t)∑<!-- ? -->i=1I∑<!-- ? -->k=1Kbiknui(t−<!-- - -->k) $$\begin{array}{} \displaystyle \boldsymbol{Y(t)}=\sum_{n=1}^{N}\boldsymbol{W}_{n}\boldsymbol{(t)}\sum^{I}_{i=1}\sum^{K}_{k=1}\boldsymbol{b}^{n}_{ik}\boldsymbol{u}_{i}\boldsymbol{(t-k)} \end{array}$$(1)

with interpolation coefficients w_n(t) based on the inverse distance between u(t) and the regime centers ^cu_n; w_n(t)=w_n(∥u(t)-^cu_n∥). For the n-th regime, bikn $\begin{array}{} \displaystyle \boldsymbol{b}^{n}_{ik} \end{array}$ are typical coefficients of the ARX model (Ljung 1987) estimated using relevant data with i indexing components of the input data vector u(t) and k denoting the backward time shift. The input vector u(t) was simply equivalent to the locally interpolated 20CRv2 wind velocity vector, whereas the UMPL-ICM wind velocity components were used as the output y(t) providing one model (1) for each output wind velocity component.

While it is perfectly possible to use one of the known identification (or optimization) algorithms to identify the values of w_n or to use TSK to scale the wind velocities, better results were achieved using fixed, predetermined regimes with centers ^cu_n and with distance-based coefficients w_n(t). Local knowledge leads to an obvious choice for regime definitions as prevalently northerly and southerly winds with the mean wind velocity for both sectors used as the regime centers.

Note that within the m-ARX framework given by (1), the regime characteristic values ^cu_n are not optimized. Optimization is, however, possible within the frames of the Takagi-Sugeno-Kang (TSK) neuro-fuzzy system where the regime of the n-th linear system is defined as a fuzzy rule (Takagi & Sugeno 1985), associating the n-th linear model with fuzzy set centers ^cu_ni, i = 1,2,…,IK, with each u_i(t – k), k = 1, …, K treated as a separate input, following the code by Jang (1993).

For each rule, ^cu_ni centers can constitute a (pseudo) vector ^cu_n = [^cu_n1, ^cu_n2, …, ^cu_nlK] whose generalized distance from the input vector u(t) serves as a basis for the interpolation between the regimes. The generalized distance between u_i(t) and the center of the fuzzy set A_ni takes the form of a membership function μ_ni(t)=μ_ni(∥u_i(t)–^cu_ni∥) (Jang 1993), usually of a Gaussian type. Within a TSK system, interpolation between different linear systems is provided by rule weights w_n(t):

Y(t)=∑<!-- ? -->n=1NWn(t)∑<!-- ? -->i=1IKbinui(t) $$\begin{array}{} \displaystyle \boldsymbol{Y(t)}=\sum_{n=1}^{N}\boldsymbol{W}_{n}\boldsymbol{(t)}\sum^{IK}_{i=1}\boldsymbol{b}^{n}_{i}\boldsymbol{u}_{i}\boldsymbol{(t)} \end{array}$$(2)

w_n(t) = w_n(μ_n1, μ_n2, …, μ_nIK), which gives an opportunity to increase the influence of the i-th input. Note that the real difference between m-ARX and TSK systems is the automatic identification of both the location of fuzzy set centers ^cu_ni and parameters of membership functions μ_ni. While the coefficients of linear models bin $\begin{array}{} \displaystyle b_{i}^{n} \end{array}$ are estimated using the same methods as the ARX model, the ^cu_ni values are estimated using the backpropagation algorithm (Jang 1993), which supports the claim to the “neuro-” part of the neuro-fuzzy system name. Furthermore, the Radial Basis Function Networks, which are considered a simplified class of neural networks, can be shown to be equivalent to a TSK system if a radial function is used as a membership function (Jang & Sun 1993).

TSK systems are mainly used to model processes occurring under separate regimes (Bakhtyar 2011; Chaudhuri 2011) and have been applied to significant wave height hindcasting (Mahjoobi et al. 2008) and prediction (Zanageh et al. 2009) and were shown to perform similarly to artificial neural networks (Malekmohamadi et al. 2011; Mahjoobi et al. 2008). In our case, the TSK systems were chosen to simulate significant wave height for their inherent ease of identifying the underlying different regimes of the system. Because regimes are not as easily discernible, as in the case of wind scaling, an automated solution using TSK seemed advantageous.

Within our study, the TSK systems used for H_s simulation employ a simulated wind velocity modulus and simulated normalized wind velocity components as an input u(t), while the WAM modeled H_s was used as an output y(t) to train the TSK. The precise composition of the input vector varied between the points of the Gulf of Gdańsk and contained one or two weighted mean wind velocity moduli from past 0–6 hours and normalized wind velocity components. For each point, up to four fuzzy sets with parameters identified automatically marked the regimes and the Gaussian membership function was used as an membership function to interpolate between them.

Significant wave height H_s parameterizations of Breugem-Holthuijsen (Breugem & Holthuijsen 2006) and the Coast Engineering Manual (USACE 2008) used the simulated wind velocity as an input (Fig. 1) and were optimized to fit the previously described high quality WAM H_s field. Both data-optimized engineering formulas (referred to as BH) proved to be equally useful as TSK, especially in conditions characterized by unique and rare wind velocities or calm weather. The output of the optimized BH parameterizations and TSK simulation were joined, depending on meteorological conditions (typical breezy to gale conditions were covered by TSK; calmer weather, extreme storms and winds from infrequent directions were covered by optimized BH parameterizations), to provide maximum accuracy and stability without the need to restrict the TSK simulation output by log-shaped restriction functions. Obviously, this procedure works on each grid point separately.

Most engineering parameterizations combine H_s time series generated under two idealized regimes of wave generation: limited by fetch X and depth ^XH_S or time ^tH_S (Massel 1996). The first regime corresponds to a scenario when waves evolving under the wind speed U₁₀ influence reach the maximum wave height permitted by depth d or fetch X, whereas the second one describes wave height evolution in time using the group velocity to convert wave evolution time into travelled distance. Because Breugem & Holthuijsen (2006) present ^XH_S parameterization as:

XHS=0.24U102gtanhk3gdU101.14tanh0.0941gXU1020.79tanhk3gdU1021.140.572 $$\begin{array}{} \displaystyle {}^{X}H_{S}=\frac{0.24U^{2}_{10}}{g}tanh\left(k_{3}\left(\frac{gd}{U_{10}}\right)^{1.14}\right) \left(tanh\left(\frac{0.0941\left(\frac{gX}{U^{2}_{10}}\right)^{0.79}} {tanh\left(k_{3}\left(\frac{gd}{U^{2}_{10}}\right)^{1.14}\right)}\right)\right)^{0.572} \end{array}$$(3)

with k₃ denoting a depth-limiting coefficient, it must be augmented by time-limited ^tH_S parametrization, which is adopted from the Coast Engineering Manual (USACE 2008). It should be noted that although these parameterizations do not directly include all swell effects, it is possible to account for a slow decay of the wave height when wind direction shifts or wind speed decreases, using a longer H_s memory. Furthermore, these parameterizations are optimized to reflect the behavior of WAM simulated fields, which is hoped to include typical swell conditions as modeled by the WAM, at least to some extent.

The above optimization procedure is carried out in stages, separately for each grid point. First, the fetch and depth limited parameterization (3) and regime relevant data are used to optimize geometric fetch values to enable realistic fetch X(t) values, which allows to include waves incoming from directions significantly different than the wind direction during the event and partly accounts for directional wave spreading. Additionally, the depth limiting coefficient k₃ needs to be optimized at that stage. During the second stage, the wind velocity is adjusted U_l0(t)=a^s ∥u(t)∥^b_s, s = 1,2 separately for the two wind velocity sectors with sector division and both a_s and b_s values determined automatically for each grid point. Subsequently, coefficients of time-limited parameterization are estimated using wind and wave data corresponding to both regimes. The time-limited parameterization is augmented with an optimized parameter for the amount of change in wind direction necessary to trigger the transition to the regime of time-limited generation, which results in the inclusion of the incoming waves from directional sectors associated with former wind directions.

Finally, both ^XH_S and ^tH_S time series are joined and coefficients for the H_s maxima enhancement are calculated separately for winds from the northern and southern sector. The difference between the definition of the sectors used for the adjustment of U₁₀ and H_s is their automatic determination in the case of U₁₀ while in the case of H_s the sectors are fixed.

Unfortunately, this procedure usually fails to provide a good estimate of H_s in swell-dominated conditions, either due to topography-induced differences between the wind directions and wave propagation or fast changing meteorological conditions. This is especially pronounced in partially sheltered areas, where a small difference in wave direction is associated with a huge fetch difference. However, the optimization of selected coefficients of the BH parameterizations, using WAM-generated significant wave heights, can reduce most consistent errors. For instance, the optimization of fetch values to fit the WAM-generated H_s may account for typical wave conditions (including the amount of swell) associated with a given wind direction, provided the data are consistent.

When a wave generation regime changes from fetch- to time-limited (e.g. due to wind direction change), a joining procedure must be applied to allow for a simultaneous slow decay of H_s values associated with the old regime and H_s built-up under a new regime, to some extent including swell.

As the optimization is carried out in stages, modifying one group of parameters at a time, there is no theoretical guarantee that a global minimum is found. However, as the search path and constraints for particular optimizations were chosen with care, using physical insight and previous runs of the genetic algorithm for a few selected points to determine typical parameter bounds, this approach may seem a reasonable compromise to decrease the computation time. A modified root mean square distance between simulated H_s and the validation set of H_S(t) is used as the optimization criterion.

To optimize the continuous parameters, the standard Matlab constrained solver (fmincon) with the interior-point algorithm is employed, using a relatively large number of iterations and a small required change of parameter values to stop the optimization process.

Having obtained all necessary coefficients, a simulation is carried out for the validation data set, and then for the entire 1871–2008 period.

The outline of the simulation follows the diagram in Figure 1. First, the 20CRv2 wind velocity components for each wet grid node were scaled down using the identified m-ARX models with better results achieved when scaling down the meridional velocity component. Next, simulated wind fields were used to simulate significant wave heights. As the TSK simulation proved to be more efficient in simulating typical stormy conditions, both the TSK and the optimized BH simulation results were joined, depending on the local wind velocity to generate the significant wave height field over the Gulf of Gdańsk for the entire 1871–2008 period.

All simulation models were identified and validated using disjoint data sets with data corresponding to WAM cold start omitted for both the identification and the validation. A short outcome of the validation process is presented in the next section “Results”.

All models were validated against a disjoint, verification-only data set of 14 storms out of 28 (1998–2001), using the root mean square error (RMSE) and the scatter index (SI) defined as the RMSE value, normalized by the mean value of validation data. The use of SI and SI-based quality measures seems to be a common technique for wave model validation as demonstrated by its application in the validation of the third generation spectral wave models: WAM (WAMDI Group, 1988), SWAN (Ris et al. 1999) and Wavewatch III (Tolman 2002). In the case of significant wave height H_s, due to the small values of mean H_s in coastal and nearshore waves, SI tends to underestimate the quality of fit (Ris et al. 1999).

The m-ARX simulated winds were validated against UMPL-ICM data. While the results of zonal wind component downscaling showed only a slight improvement, the scatter index (SI) values for the meridional component were significantly reduced from the maximum of 51% to a somewhat more acceptable level of 31%.

For brevity, only the validation results for the significant wave height simulation using TSK will be presented here, as values generated by this method are most likely associated with typical stormy conditions.

The scatter indices (SI) and the root mean square errors (RMSE) for each point are presented in Figure 3 and the comparison of significant wave height modeled with WAM and simulated with TSK during selected storm maxima is presented in Figure 4. Note that within the most western, sheltered part of the Gulf (Puck Bay), the scatter index values range from 5 to 27%, whereas in the open part of the Gulf, the SI values are almost constant within the range of 21–31% as observed in Figure 3. The highest SI values are found in the area occasionally protected by the Hel Peninsula and reach 50% in the nearshore zone. On the other hand, the comparison of significant wave height during storm maxima (Fig. 4) shows that the simulation retains spatial patterns, though mostly underestimates the H_s values.

Figure 3

Figure 4

As the quality of both simulation methods is similar, the simulation error using both methods simultaneously does not exceed 30% of the simulated value in most areas of the Gulf of Gdańsk. Taking into account the discrepancy between the downscaled 20CRv2 wind fields and the UMPL-ICM wind used to run the WAM model generating the validation H_s, such an error is acceptable, especially if one is interested in the long-term statistics of the dataset and not in the actual values of H_s. It is worth noting that to capture a change in storminess, it is enough to obtain a good estimate of N_H – the monthly number of hours with H_s exceeding a given threshold that is specific for each point.

In this paper, the number of hours N_H(t) with H_s exceeding the point-specific threshold for each month is used as a measure of storminess. For each gridpoint, the threshold value is chosen as the H_s value exceeded with a probability of 0.05 in the WAM-generated data (1998–2001). The N_H values are calculated for each month and form a single time series mostly to filter out the daily variability and focus on the frequency of more significant stormy events. In this approach, the seasonal (or monthly) variability of N_H over a typical year is not considered, neither is the long-time variability for a given month. The cumulative trend is calculated as a difference between the final and the initial values of the trend line fitted to N_h.

To compare the cumulative changes in N_H attributed to the linear trend at different points more easily, these changes are normalized to the mean N_H value at a given point to form a map in Figure 5. Unfortunately, the N_h, as defined above, is not an outstandingly consistent measure of the frequency of heavy sea conditions across all points of the Gulf of Gdańsk, especially nearshore where waves are depth limited. However, it gives similar estimates of N_H for points with a similar shape of the cumulative distribution function for H_s. In these cases, the estimates are more consistent than when N_H was defined using a fixed, common threshold for all points of e.g. 1 m. In this way, changes for each point in the Gulf of Gdańsk, when normalized with a point-mean value of N_H can be compared with each other at the expense of information on absolute values. The authors believe that the presentation method applied in this study is advantageous as the commonly used change in mean monthly H_s values, though interesting and readily understood, cannot be used to present spatial patterns in long-term trends. For example, a change of 0.2 m has a completely different meaning for points with a mean H_s value of 1 m and, for example, 0.6 m.

Figure 5

As shown in Figure 5, the 1871–2008 cumulative change in N_H is negative in the western part of the Gulf of Gdańsk (min. 25% of the mean value) and positive in the eastern part (max 80% of the mean value).

It is expected that further research will elaborate on long-term oscillation spatial patterns and provide a more detailed and coherent spectral description of N_h.