




\section{Introduction}

The initial motivation for developing the methods presented in this chapter was the search for a practical solution to the air gap problem for offshore structures. Over the years several reports had been filed detailing damage of the deck structure of offshore platforms due to wave impacts. In fact, such damage seemed to occur more frequently than could be expected from predictions based on standard theory. These wave impacts are highly undesirable events as they may, in the worst case, compromise the structural integrity of the platform deck. Hence, considerable attention was being paid to this problem over some time, and now reasonably good predictive tools have become available, as will be demonstrated in this chapter. 

It was eventually realized that the standard procedure for predicting the extreme wave crest height to be expected at a given platform location is largely based on extreme value statistics for a single point \citep{forr06}. It is clear that if this is used as a basis for predicting the probability of wave impact on the deck structure, then the fact that there is also a significant area effect on the extreme crest height distribution over the deck area is neglected. Hence, for a proper solution of the air gap problem, it is necessary to model the ocean surface as a random field so that also the spatial aspect can be correctly dealt with \citep{forr06}. For the case of a homogeneous Gaussian random wave field, this aspect is also discussed by \citet{juglard05,bary06}. 

Recently, a new, simplified method for predicting the space-time extreme value statistics of homogeneous Gaussian random fields has been proposed \citep{naba10pem}. In the present chapter it is shown that this simplified method can be extended to deal with homogeneous non-Gaussian random fields, and, in particular, with the special case of a second-order homogeneous ocean wave field \citep{naba12pem}. Additionally, a new semi-parametric representation of the space-time extreme value distribution for quite general homogeneous random fields over rectangular domains is proposed. It is demonstrated that both the simplified method as well as the new semi-parametric method makes it possible to make good predictions of the extreme crest heights over the deck area for jacket platforms, which are assumed to be transparent for the big waves. Even if the second-order ocean wave field is intrinsically non-Gaussian, its deviation from a Gaussian field is not very significant. Therefore, to implement a more stringent test of the two methods proposed, an example of a strongly non-Gaussian field has also been included. Both methods will be shown to provide accurate extreme value predictions for this case as well.  

The proposed simplified approach, as opposed to the new semi-parametric method, can be easily modified to also provide good solutions of the air gap problem for a TLP or a semi-submersible platform, provided measured or simulated time series of the ocean surface elevation or air gap at a sufficient number of points below the deck structure are available, which properly account for the motions of the platform and the effect of the structure on the wave field. Due to the interaction effects between the wave field and the structure, the wave field is not homogeneous in space under the platform deck, even if the incoming wave field is of this type. Hence, a prediction method that can deal with a random field that is nonhomogeneous in space while stationary in time will be required. The simplified method proposed in this chapter can provide such a prediction tool.  

Over the last 20 years or so, several approximations for the extreme value distributions of a real-valued random field
$X(\mathbf{x},t)$, where $\mathbf{x} = (x_1, \ldots, x_d) \in
\mathcal{D} \subseteq \mathbb{R}^d$ and $t \in \mathcal{T}
\subseteq \mathbb{R}$ have been proposed, of which a few are mentioned, \citep{adle81,vanm83,sun93,mabr97,ditl04}. Of particular interest are
two books which provide substantial theoretical results, mainly for Gaussian random fields \citep{pite96,adta07}.




\section{Spatial-Temporal Extremes for Gaussian Random Fields}\index{spatial-temporal extremes!Gaussian random fields}


Let $X = X(\mathbf{x},t)$ denote a zero mean random field defined
on $\mathbb{R}^d \times \mathbb{R}$, and let,
\begin{equation}\label{eq:st1}
\hat{X} = \hat{X}(\mathcal{D} \times \mathcal{T}) =
\max\{ X(\mathbf{x},t); \mathbf{x} \in \mathcal{D}\,,
\: t \in \mathcal{T}\} \,,
\end{equation}
where $\mathcal{T}$ denotes a time interval, specifically, $\mathcal{T}
= (0,T)$ for a given $T$.

In Chapter~14 of \cite{adta07}, the following result is proved for a smooth
Gaussian random field,
\begin{equation}\label{eq:st2}
\big|\, \Prob (\hat{X} \ge \xi ) - \EE[\phi(A_{\xi}(X,
\mathcal{D} \times \mathcal{T}))] \,\big| < O \big( e^{- \alpha
	\xi^2/(2 \sigma^2)} \big) \,,
\end{equation}
where $\phi(A_{\xi}(X,
\mathcal{D} \times \mathcal{T}))$ denotes the Euler characteristic of $A_{\xi}(X,
\mathcal{D} \times \mathcal{T})$, which is the excursion set of $X$ over
$\mathcal{D} \times \mathcal{T}$. That is, $A_{\xi}(X, \mathcal{D}
\times \mathcal{T}) = \{ (\mathbf{x}, t) \in \mathcal{D} \times
\mathcal{T} : X(\mathbf{x},t) \ge \xi  \}$, $\sigma^2$ is the
variance of $X$ (assumed constant), and $\alpha > 1$ is a constant
that can be identified as discussed by \cite{adta07}. Certain restrictions apply to the geometry of the domain $\mathcal{D}$. Instead of going into details about
this formula, a largely equivalent result will be given, which is due
to \citet{pite96}. This result applies to a homogeneous Gaussian field, which
has a simple geometry for $\mathcal{D}$, typically a nice
convex domain e.g. like a rectangle or a sphere. Any numerical
representation of the homogeneous random field $X$ by necessity
leads to a truncation of the spectral density (assumed to exist).
Hence, there is no practical limitation in assuming that the
spectral density of $X$ has bounded support, which leads to
smoothness of any order. Piterbarg's (asymptotic) formula, cf. Chapter 1 of \citep{pite96}, can then
be written as ($\xi = u \sigma$, $u \to \infty$ ),
%
\begin{equation}\label{eq:st3}
\Prob(\hat X \leq \sigma u)  \simeq  P_A(u) =
\exp\left\{-\varphi(u) {(2\pi)^{(d+1)/2}H_{d}(u)} V_d\right\}
, %\nonumber
\end{equation}
%
where $\simeq$ denotes asymptotically equal, $H_k(u)=(-1)^k \varphi(u)^{-1} ({\rm d}/{\rm d}u)^k \varphi(u)$ is
a Hermite polynomial of order $k$, $\varphi$ is the standard
normal probability density, that is, $\varphi(u) = (\sqrt{2\pi})^{-1} e^{-u^2/2}$, 
$V_d= m({\mathcal D}\times{\mathcal T})/m({\mathcal V})$, $m({\mathcal
V})=(2\pi\sigma)^{d+1}/\sqrt{\det\Lambda}$ ($\Lambda$ is defined by Eq.~\eqref{eq:lambda} below), where $m({\mathcal A})$ denotes the Lebesgue measure (volume) of the domain $\mathcal A$.

Before the calculation of $V_d$ is discussed any further, a slight detour is made and
some basic results are presented as detailed by \cite{krog04}. Denoting
the spectral density of $X$ by $\Psi(\mathbf{k}, \omega)$, its
covariance function $\gamma(\mathbf{x}, t)$ is given as,
\begin{equation}\label{eq:st4}
\gamma(\mathbf{x}, t) = \int_{\mathbb{R}^{d+1}} e^{\ii (\mathbf{k}
	\mathbf{x}' + \omega t)} \, \Psi(\mathbf{k}, \omega)  \, d
\mathbf{k} d \omega,
\end{equation}
where $\mathbf{k}$ denotes the wave number vector.

Denote by $\Lambda=(\Lambda_{ij})$ the essentially nonsingular
covariance matrix of the gradient $\nabla X({\bf
	x},t)=\left(\partial X/\partial{x_1}, \ldots, \partial X/\partial{x_d}, \partial
X/\partial{t}\right)$ evaluated in the point ${\bf x=0}$, $t=0$,
\begin{align}\label{eq:lambda}
\Lambda_{ij} &= {\rm Cov}\left(\nabla_{\!i}\, X({\bf 0},0),
\nabla_{\!\!j}\, X({\bf 0},0)\right)\nonumber \\
&= - \nabla_{\!i}\, \nabla_{\!\!j}\, \gamma({\bf 0},0)
= \int_{\mathbb{R}^{d+1}} r_i r_j \Psi ({\bf r}) \, {\rm d}{\bf r},
\end{align}
where ${\bf r}=(r_1, \ldots, r_d, r_{d+1})=(k_1, \ldots, k_d, \omega)$.


A linear space transformation $A:{\mathcal
	D}\times{\mathcal T}\to \tilde{\mathcal D}\times\tilde{\mathcal T}$ is now performed with the
help of a nonsingular $(d+1)\times (d+1)$ matrix $A$,
\begin{equation}
{\bf \tilde{y}}= A{\bf \tilde{x}}\,,
\end{equation}
where ${\bf \tilde{x}} = ({\bf x}, t) = (x_1, \ldots, x_d, t)$, then a new Gaussian random field $\tilde X({\bf \tilde{y}})$ defined as,
\begin{equation}
\tilde X({\bf \tilde{y}}) = X(A^{-1}{\bf \tilde{y}})/\sigma
\label{Y}
\end{equation}
has unit variance, $\tilde \sigma=1$, and the covariance matrix of the gradient of $\tilde X$ is
given by,
\begin{equation}
\tilde \Lambda = \frac{1}{\sigma^2} (A^{-1})^T \Lambda A^{-1} \,.
\label{lambda_Y}
\end{equation}

Since $\Lambda$ is a symmetrical matrix, the transformation
$A = \sqrt{d+1} \sigma^{-1} \, \Lambda^{1/2}$ is well defined and
results in $\tilde\Lambda = (d+1)^{-1}I$. For this special case
of the Gaussian random field $\tilde X({\bf \tilde{y}})$ corresponding to the
mentioned covariance matrix $\tilde \Lambda$, a complementary
extreme value distribution can be obtained for the particular case of a rectangular domain $\Delta$ with edges
oriented along the coordinate axes, that is, $\Delta=\otimes_{i=1}^{d+1} [0,
L_i]$, where $\otimes$ denotes a Cartesian product. This probability can be written as,
\begin{eqnarray}\label{eq:st5}
\Prob(\hat{ \tilde X} > u )  &\simeq&
\varphi(u) \sum_{k=1}^{d+1} \frac{H_{k-1}(u)}{ [2\pi (d+1)]^{k/2}} \Sigma_k({\bf L}) + \int_u^\infty{\varphi(x) {\rm d} x} \,,
\end{eqnarray}
where ${\bf L}=(L_1,\dots, L_{d+1})$ and $\Sigma_k({\bf L})$ is
the $k$-th elementary symmetric polynomial, defined as,
\begin{equation}\label{Sigma}
\Sigma_k({\bf L})=\sum\limits_{1\leq j_1<j_2<\dots<j_k\leq d+1}{
	L_{j_1}\cdot L_{j_2}\cdots L_{j_k}}.
\end{equation}
This result is obtained as a special case of Theorem~5.1 in the book by \cite{pite96}.

Now, let us rewrite Eq.~(\ref{eq:st5}) for the case of the Gaussian
field $X({\bf x}, t)$ having variance $\sigma$ and diagonal
covariance matrix of general form, $\Lambda_X={\rm
	diag}(\lambda_1,\dots,\lambda_{d+1})$. This case is important since an
arbitrary Gaussian field $Z({\bf z}, t)$ with covariance matrix
$\Lambda_Z$ can be brought into this form by performing a rotation $A$
($A^T=A^{-1}$) of the coordinate axes to the principal directions. As a consequence, only rectangular domains with edges oriented along principal directions are considered. 

It is now easy to check that for $X({\bf x}, t)$,
\begin{eqnarray}\label{eq:st6}
\Prob(\hat{  X} > u\sigma)  &\simeq& \varphi(u) \sum_{k=1}^{d+1}
{(2\pi)^{k/2}H_{k-1}(u)} \Sigma_k({\bf q}) +
\int_u^\infty{\varphi(t) {\rm d} t},
\end{eqnarray}
where ${\bf q}=\left( q_1, q_2, \cdots, q_{d+1} \right)$,
$q_i = L_i/ \ell_i^0$, $\ell_i^0={2\pi\sigma/\sqrt{\lambda_i}}$,
$i=1,\dots,d+1$.

For high levels $u$, and, correspondingly, small exceedance
probabilities, one may write,
\begin{eqnarray}\label{eq:st7}
\Prob(\hat{ X} \leq u\sigma)  &=& 1 - \Prob(\hat{X} >
u\sigma) \simeq \exp\left\{- \Prob(\hat{X} >  u\sigma)\right\}
\\ \nonumber
&\simeq& P_B(u) = \exp\left\{ - \varphi(u) \sum_{k=1}^{d+1}
{(2\pi)^{k/2}H_{k-1}(u)} \Sigma_k({\bf q}) \right\} \, .
\end{eqnarray}

This formula (referred to in the following as PB) is in
agreement with the ``conventional'' Piterbarg's asymptotic formula\index{Piterbarg's asymptotic formula}
Eq.~\eqref{eq:st3} (referred to as PA) if only the leading term of the sum in Eq.~(\ref{eq:st7}) are taken into
account, as $\Sigma_{d+1}({\bf
	q})=V_d$. The leading term in the exponent of PB (and the only one
in PA) provides the probability of the exceedance from inside the
volume of the domain ${\mathcal D}\times{\mathcal T}$, while other terms
of PB represent probability of exceedance for the
lower-dimensional manifolds associated with the boundaries of the
domain, like faces, edges, and vertices in the case of 3
dimensional $\Delta$. Therefore, if one wants to quantify the
exceedance probability for a rectangle $\Delta$, one of whose spatial sizes, say
$L_1$, tends to zero, PA must be rewritten for the lower-dimensional
case. This demands in particular using $V_{d-1}$ instead of
$V_{d}$, as the last tends to zero, while the same formula PB will
work for all ranges of $\Delta$, that is, there is no need to change it upon
a ``squeeze'' of one of the dimensions.

To check the validity of both PA and PB, a
comparison of them will be made against numerical results for (1+1)- and
(2+1)-dimensional Gaussian fields. For the (2+1)-dimensional Gaussian
field, the explicit expression for PB is,
\begin{eqnarray}\label{eq:st8}
%P(\hat{ X} \leq u\sigma)  &\sim&
P_B(u)  &=&
\exp\big\{ - \exp(-u^2/2)\big[q_x + q_y + q_t \nonumber \\
&&+\sqrt{2\pi}\,u\,(q_x \,q_y + q_x \, q_t + q_y \, q_t) + 2\pi \, q_x\, q_y \,q_t \,(u^2-1)
\big] \big\} \,. 
\end{eqnarray}
Here $q_x = q_1 = L_1/ \ell_1^0$, $q_y = q_2 = L_2 / \ell_2^0$, $q_t = q_3 = L_3 / \ell_3^0$, and 
$\ell_1^0=2\pi\sigma/\sqrt{\Lambda_{xx}}$,
$\ell_2^0=2\pi\sigma/\sqrt{\Lambda_{yy}}$,
$\ell_3^0=2\pi\sigma/\sqrt{\Lambda_{tt}} = T_z = 2\pi/\omega_z$,
where,
\begin{equation}\label{eq:st9}
\left(\begin{array}{c}\Lambda_{xx} \\ \Lambda_{yy}
\end{array}\right) = \int{ \left(\begin{array}{c} k_x^2  \\ k_y^2
	\end{array}\right)
	\Psi({\bf k},\omega)\, {\rm d}\omega\, {\rm d}{\bf k} } \,,
\end{equation}
\begin{equation}\label{eq:st10}
\Lambda_{tt} = \omega_z^2 \sigma^2 = \int{ \omega^2\,
	\Psi({\bf k},\omega)\, {\rm d}\omega\, {\rm d}{\bf k} } \,.
\end{equation}
The last equality in Eq.~\eqref{eq:st10} defines the zero-upcrossing (circular) frequency $\omega_z$.
A more general result corresponding to the special case of Eq.~\eqref{eq:st8} with $q_t = 0$ was derived already by \citet{ditl71}.


It can be shown that the rhs of Eq.~\eqref{eq:st3} has the following Gumbel
distribution as its max-stable asymptotic limit \citep{krog04},
\begin{equation}\label{eq:st130}
G(u \sigma) = \exp \big\{ - \exp \big[ - h_V (u - h_V) \big] \big\} \,,
\end{equation}
where the parameter $h_V$ is obtained by solving the equation,
\begin{equation}\label{eq:st50}
V_d \, h_V^2 \, e^{- h_V^2/2} = 1 \,,
\end{equation}
which has the approximate solution,
\begin{equation}\label{eq:st60}
h_V = \sqrt{2 \log V_d + 2 \log(2\log V_d)} \,.
\end{equation}



\section{A Simplified Approach}


The first step is rewriting Eq.~(\ref{eq:st1}) as follows,
\begin{equation}\label{eq:st11}
\hat{X} = \underset{0 \le t \le T}{\max} \, \{ \underset{\mathbf{x} \in \mathcal{D}}{\max}\, X(\mathbf{x},t) \} \,.
\end{equation}
A stochastic process consisting of spatial extremes is then introduced,
\begin{equation}\label{eq:st12}
\hat{X}_{\mathcal{D}}(t) = \underset{\mathbf{x} \in \mathcal{D}}{\max}\, X(\mathbf{x},t) \,.
\end{equation}
It is now assumed that a finite set of points $\{ \mathbf{x}_i \}_{i = 1}^M$ can be chosen so that
\begin{equation}\label{eq:st13}
\hat{X}_{\mathcal{D}}(t) \approx \underset{ 1 \le i \le M }{\max} X(\mathbf{x}_i, t) \,,
\end{equation}
for each $t \in (0, T)$, and that,
\begin{equation}\label{eq:st14}
\hat{X} \approx \underset{1 \le n \le N}{\max} \hat{X}_{\mathcal{D}}(t_n) \approx \underset{1 \le n \le N}{\max} \, \{ \underset{ 1 \le i \le M }{\max} X(\mathbf{x}_i, t_n) \} \,,
\end{equation}
for a suitable discretization of the time interval $(0, T)$, $0 \le t_1 < \ldots < t_N \le T$. For each $n$, $n = 1, \ldots, N$, denote by $\mathbf{x}_{i_n}$ a point so that $\underset{ 1 \le i \le M }{\max} X(\mathbf{x}_i, t_n) = X(\mathbf{x}_{i_n}, t_n)$, and let $\tilde{X}_n = X(\mathbf{x}_{i_n}, t_n)$. Then
\begin{equation}\label{eq:st15}
\hat{X} \approx \underset{1 \le n \le N}{\max} \, \tilde{X}_n \,,
\end{equation}

If the sampling times $t_n$ are sufficiently dense in $(0,T)$, then the
time series $\tilde{X}_n$ can be considered as quasi-continuous in time,
and an accurate discrete representation of a smooth stochastic process
$\tilde{X}(t)$ so that $\hat{X} \approx \underset{0 \le t \le T}{\max} \, \tilde{X}(t)$.

The derivations in this section applies also to the case of a nonhomogeneous field, and the
extreme value distribution $F_{\hat{X}}(\xi)$ is then with good approximation given by the equation,
\begin{equation}\label{eq:st16}
F_{\hat{X}}(\xi) \approx \exp \left\{ - \int_0^T \nu_{\tilde{X}}^+(\xi;t) \, dt \right\} \,,
\end{equation}
where $\nu_{\tilde{X}}^+(\xi;t)$ denotes the average upcrossing rate of the level $\xi$ by $\tilde{X}(t)$ at time $t$.
This approximation is contingent on the assumption that the probability of initial exceedance can be neglected. 
For estimation purposes, Eq.~\eqref{eq:st16} is rewritten as,
\begin{equation}\label{eq:st17}
F_{\hat{X}}(\xi) \approx \exp \{ - \overline{\nu}_{\tilde{X}}^+(\xi) \, T \} \,,
\end{equation}
where the time average $\overline{\nu}_{\tilde{X}}^+(\xi)$ is given as,
\begin{equation}\label{eq:st18}
\overline{\nu}_{\tilde{X}}^+(\xi) =  \, \frac{1}{T} \, \int_0^T \nu_{\tilde{X}}^+(\xi;t) \, dt \,,
\end{equation}
which is in a form suitable for empirical estimation from time series. However, it should be clear that in the case of a nonhomogeneous field, the extreme value predictions can only refer to the initial space-time domain $\mathcal{D} \times \mathcal{T}$, if no additional modelling or structure is introduced.  

The comparison between Eqs.~(\ref{eq:st3}), (\ref{eq:st7}), and
(\ref{eq:st17}) reveals that $-\log(P_B(\xi/\sigma))$ is not a
multiplicative function with respect to $T$ (or, equivalently,
$q_t$) in contrast to $-\log(P_A(\xi/\sigma))$ and
$-\log(F_{\hat{X}}(\xi))\approx \overline{\nu}_{\tilde{X}}^+(\xi)
\, T$. More precisely,
%
%
%
\begin{align}\label{eq:st19} %%
P_A(\xi/\sigma) &= \exp \{ - \nu_{A}^+(\xi)\, T \} \,,
\end{align}
\begin{align}\label{eq:st20} %%
P_B(\xi/\sigma) &= \exp \{ - \nu_{0}^+(\xi) - \nu_{B}^+(\xi)\, T
\} \,,
\end{align}
%
where, for the 3-dimensional case,
%
\begin{align}\label{eq:st21}%%
\nu_{A}^+(u\sigma) &= 2\pi \, q_x \, q_y \,(u^2-1) \exp(-u^2/2) / T_z
,\nonumber\\
%\nu_{B}^+(u\sigma)&=&\left[1+\sqrt{2\pi}\, u(n_x+n_y) + 2\pi
%n_xn_y(u^2-1)\right] \exp(-u^2/2) / T_z
%\\
\nu_{B}^+(u\sigma) &= \nu_{A}^+(u\sigma)+ \left[1+\sqrt{2\pi}\,
u(q_x + q_y) \right] \exp(-u^2/2) / T_z ,
\\ \nonumber%%
\nu_{0}^+(u\sigma) &= (q_x + q_y +\sqrt{2\pi}\, u \,q_x\, q_y) \exp(-u^2/2)\, ,
\end{align}
%
as it immediately follows from Eq.~(\ref{eq:st8}) since $q_t = L_3/\ell_3^0 = T/T_z$. The term
$\nu_{0}^+(\xi)$, the only one which survives in PB when $T\to 0$,
gives the probability of an initial exceedance event on the 'boundary'
$\mathcal{D}$ of the $\mathcal{D\times T}$ domain. Since the most
interesting cases in practice would typically have domains where the
time dimension is much greater then the spatial ones, which leads to $\max(q_x, q_y) \ll q_t$, 
$\nu_{0}^+(\xi)$ will be safely neglected for the numerical examples in this chapter. A similar approximation was also implemented in Eq.~\eqref{eq:st16}.


From  Eqs.~(\ref{eq:st3}) and (\ref{eq:st130}), it is seen that the mean
upcrossing rate tail, say for $\xi \ge \xi_0$, behaves in a manner
largely determined by a function of the form $\exp\{-a ( \xi -
b)^c\}$ ($\xi \ge \xi_0$), where $a$, $b$ and $c$ are suitable
constants. It is therefore assumed that the mean upcrossing rate function of $\tilde{X}(t)$ can be represented as,
\begin{align}\label{eq:st22}
\nu^+(\xi) \approx q(\xi) \, \exp\{- a (\xi - b)^c \} \, , \: \xi
\ge \xi_0\,,
\end{align}
for a suitable choice of $\xi_0$, where the function $q(\xi)$ is slowly varying compared with the
exponential function $\exp\{- a (\xi - b)^c \}$ for tail values of
$\xi$, cf. \cite{naga08jem}. Now, typically, the function $q(\xi)$ can be largely considered as a constant
for tail values of $\xi$. This suggests an extrapolation
strategy obtained by replacing $q(\xi)$ by a suitable constant
value, $q$ say.

The adopted procedure for identifying appropriate values for the parameters
$a, b ,c, q$, assuming a constant $q$, follows closely the optimization method detailed in Chapter~\ref{chap:acer}. It is based on minimizing the following
mean square error function with respect to the four arguments,
\begin{equation}\label{eq:st23}
F(q, a, b, c) = \sum_{j = 1}^N w_j\, \big| \log \hat{\nu}^+(\xi_j)
- \log q + a (\xi_j - b)^{c} \big|^2\,,
\end{equation}
where $\hat{\nu}^+(\xi_j)$ is the empirical estimate of the upcrossing rate at
the level $\xi_j$ ($\xi_0 \le \xi_1 \le \ldots \le \xi_N$), $w_j$ denotes a weight factor that puts more emphasis on the
more reliable estimates of $\hat{\nu}^+(\xi_j)$. The choice of weight factor is to some
extent arbitrary. Here, a weight factor will be used based on the
confidence interval associated with the empirical estimate of the upcrossing
rate $\hat{\nu}^+(\xi_j)$.




\section{Spatial-Temporal Extremes for Non-Gaussian Random Fields}\index{spatial-temporal extremes!non-Gaussian random fields}\label{sec:stengrf}

In this section, the discussion is limited to random wave fields with two spatial dimensions. However, it should be quite obvious how to extend the methods discussed to higher dimensional spaces. Accordingly, let $X = X(\mathbf{x},t)$ denote a zero mean random field defined
on $\mathbb{R}^2 \times \mathbb{R}$, and let
\begin{align}\label{eq:st26}
\hat{X} = \hat{X}(\mathcal{D} \times \mathcal{T}) =
\max\{ X(\mathbf{x},t); \mathbf{x} \in \mathcal{D}\,,
\: t \in \mathcal{T}\} \,,
\end{align}
where $\mathcal{T}$ denotes a time interval, specifically, $\mathcal{T}
= (0,T)$ for a given time $T$, and $\mathcal{D}$ is the area of interest. 

Except for the case of homogeneous Gaussian fields, and memoryless transformations of them, the asymptotic distribution function of $\hat{X}$ is unknown \citep{pite96,adta07}. In this section, the focus is on the case of second-order homogeneous wave fields. Since they are non-Gaussian, and cannot be obtained by a simple transformation of a Gaussian field, the standard theory of Gaussian fields does not apply. However, even if no simple formulas can be derived for extremes of second-order wave fields, the simplified empirically based approach described in the previous section, may still be applied, cf. \citep{naba10pem}. 


The points $\{ \mathbf{x}_i \}_{i = 1}^M$ of Eq.~\eqref{eq:st13} for the case of a wave field are typically chosen as the nodal points of a rectangular grid with a mesh size determined by the smallest significant wave length of the field. A mesh size equal to one tenth of this smallest wave length would be sufficient. Then to ensure the approximation in Eq.~\eqref{eq:st14}, the (equidistant) discretization of the time interval $(0, T)$, $0 = t_0 < t_1 < \ldots < t_N \le T$. $\Delta t = t_{j} - t_{j-1}$ can in most practical cases be chosen as one tenth of the spectral peak period. The robustness of the results with respect to discretization can be verified by refining the mesh size in space and time. 


Although the main focus of this section is on second order ocean waves, the applicability of the methods discussed are far more wide ranging. However, for more general cases the validity of the approximation in Eq.~\eqref{eq:st14} may require a more careful consideration. For a homogeneous random field which has a spectral density of bounded support, which can often be assumed in applications, then the field is infinitely smooth and there will be uniform bounds on the rate of change of the field. Hence, a finite grid will always exist that secures the validity of the approximation in Eq.~\eqref{eq:st14} within a given accuracy. For more general cases, assuming that the field has continuously differentiable realizations with uniformly bounded derivatives on bounded domains, then again the validity of Eq.~\eqref{eq:st14} within a specified level of accuracy is guaranteed for a suitable choice of the grid. Usually in practice, physical rather than mathematical arguments can be used to justify this approximation. 



An alternative procedure is now introduced based on a detailed study of the extreme value statistics of Gaussian random fields described by \cite{pite96,adta07}. Consider a rectangular domain $\mathcal{D} = \mathcal{L}_x \times \mathcal{L}_y$, where $\mathcal{L}_x$ denotes an interval of length $L_x$, and similarly, $\mathcal{L}_y$ denotes an interval of length $L_y$. Neglecting the probability of initial exceedance, which is acceptable for long time intervals, it is then postulated that the distribution of $\hat{X}(\mathcal{D} \times \mathcal{T})$ for a homogeneous random field $X$ can be written as,
\begin{align}\label{eq:st35}
F_{\hat{X}}(\xi) =  \exp \big\{ - \big( \nu_0^+(\xi) + \nu_x^+(\xi) L_x + \nu_y^+(\xi) L_y + \nu_{xy}^+(\xi) L_x L_y \big) \, T \big\} \,.
\end{align}
Here $\nu_0^+(\xi)$ denotes the one-point upcrossing rate of the level $\xi$. $\nu_x^+(\xi)$ ($\nu_y^+(\xi)$) represents corrections to the one-point upcrossing rate due to extensional effects in the $x$ ($y$) direction, while $\nu_{xy}^+(\xi)$ corrects for area effects. In particular, $\nu_0^+(\xi)$ is the mean upcrossing rate of the stochastic process $X(\mathbf{x},t)$ for an arbitrary, fixed point $\mathbf{x}$. $\nu_x^+(\xi)$ is determined by introducing the stochastic process $\hat{X}_{\mathcal{L}_x}(t) = \underset{\mathbf{x} \in \mathcal{L}_x \times \{y_0\}}{\max} X(\mathbf{x},t)$, for a suitably fixed value $y_0$, and writing the mean upcrossing rate of this process as $\nu_0^+(\xi) + \nu_x^+(\xi) L_x$. This defines the mean normalized upcrossing rate $\nu_x^+(\xi)$. 
Similarly, $\nu_y^+(\xi)$ is determined by introducing the stochastic process $\hat{X}_{\mathcal{L}_y}(t) = \underset{\mathbf{x} \in \{x_0\} \times \mathcal{L}_y}{\max} X(\mathbf{x},t)$, for a suitably fixed value $x_0$, and writing the mean upcrossing rate of this process as $\nu_0^+(\xi) + \nu_y^+(\xi) L_y$. This, then, defines the mean normalized upcrossing rate $\nu_y^+(\xi)$. Finally, $\nu_{xy}^+(\xi)$ is determined by writing the mean upcrossing rate of the process  $\hat{X}_{\mathcal{D}}(t) = \hat{X}_{\mathcal{L}_x \times \mathcal{L}_y}(t) = \underset{\mathbf{x} \in \mathcal{L}_x \times \mathcal{L}_y}{\max} X(\mathbf{x},t)$ as $\nu_0^+(\xi) + \nu_x^+(\xi) L_x + \nu_y^+(\xi) L_y + \nu_{xy}^+(\xi) L_x L_y $. Combining this with the previous relations provides a way to obtain $\nu_{xy}^+(\xi)$. These functions can be estimated in a manner similar to the procedure to be described next. The advantage of the representation provided by Eq.~\eqref{eq:st35} is that the calibration of the right hand side needs to be done for only one rectangular domain. After that, the formula can then be used for any other rectangular domain (with the same orientation). 

While the procedure for estimating the upcrossing rate functions described above, is one way to achieve this, a somewhat different approach is used here. Instead of describing the procedure in full generality, it will be detailed as applied to the specific example in the section on numerical examples. The edges of the chosen rectangular domain for the calibration, which was 300m$\times$300m, were divided into subintervals of length 25m. An aggregate of rectangular subdomains $\mathcal{D}_{mn}$, $m, n = 0, 1, \ldots, 12$, where $\mathcal{D}_{mn} = (0,\, 25 \cdot m) \times (0, \, 25 \cdot n)$, is then created. Note that this includes the degenerate domains $(0,0)$, $(0, \, 25 \cdot m ) \times (0, 0)$ and  $(0, 0) \times (0, \, 25 \cdot n)$.  A linear regression approach based on the upcrossing rate of the area extremes process for each of the resulting rectangular domains $\mathcal{D}_{mn}$ obtained for $m,n = 0, 1, \ldots, 12$ was used to estimate $\nu_0^+(\xi_j), \nu_x^+(\xi_j), \nu_y^+(\xi_j), \nu_{xy}^+(\xi_j)$ for each value of $\xi_j$, $j = 1, \ldots, J$, which denotes a preassigned range and number of $\xi$-levels leading to meaningful estimates of the various upcrossing rates from the data. 


In the following, $\nu^+(\xi)$ will be used as a generic notation denoting either the mean or the time averaged upcrossing rate, as the case may be. A key issue for prediction of extreme values is the estimation of $\nu^+(\xi)$. Also for the case considered in this section, it is assumed that the mean upcrossing rate tail, say for $\xi \ge \xi_0$, behaves in a manner largely determined by a function of the form $\exp\{-a ( \xi - b)^c\}$ ($\xi \ge \xi_0$) where $a$, $b$ and $c$ are suitable constants. Consequently, it is assumed that the mean upcrossing rate function of $\tilde{X}(t)$ can be represented as,
\begin{align}\label{eq:st36}
\nu^+(\xi) \approx q(\xi) \, \exp\{- a (\xi - b)^c \} \, , \: \xi
\ge \xi_0\,,
\end{align}
for a suitable choice of $\xi_0$, where the function $q(\xi)$ is slowly varying compared with the
exponential function $\exp\{- a (\xi - b)^c \}$ for tail values of
$\xi$, cf. \citep{naga08jem}. The function $q(\xi)$ can be largely considered as a constant
for tail values of $\xi$. This points to an extrapolation
strategy based on replacing $q(\xi)$ by a suitable constant value, $q$ say.

The adopted procedure for identifying appropriate values for the parameters
$a, b ,c, q$, assuming a constant $q$, is largely identical to that of the previous section.




\section{Empirical Estimation of the Mean Upcrossing Rate}\label{transf}\index{mean upcrossing rate!empirical estimation}


In the previous section it was shown that the key to providing
estimates of the extreme values of the response process $X(t)$ on
the basis of simulated response time histories, is the estimation
of the mean upcrossing rate. By assuming the requisite ergodic
properties of the response process for a short-term condition, the
mean upcrossing rate is conveniently estimated from the ergodic
mean value. That is, it may be assumed that,
\begin{align}\label{eq:emcr2}
\nu^+(\xi) = \lim_{t \to \infty} \frac{1}{t} \, n^+(\xi;0,t) \,,
\end{align}
where $n^+(\xi;0,t)$ denotes a realization of $N^+(\xi;0,t)$, that
is, $n^+(\xi;0,t)$ denotes the counted number of upcrossings
during time $t$ from a particular simulated time history for which
the starting point $t=0$ is suitably chosen. In practice, $k$ time
histories of a specified length, $T_0$ say, are simulated. The
appropriate ergodic mean value estimate of $\nu^+(\xi)$ is then
\begin{align}\label{eq:emcr3}
\hat{\nu}^+(\xi) = \frac{1}{k \,T_0} \sum_{j=1}^k  \,
n_j^+(\xi;0,T_0) \,,
\end{align}
where $n_j^+(\xi;0,T_0)$ denotes the counted number of upcrossings
of the level $\xi$ by time history no. $j$. This will be the
approach to the estimation of the mean upcrossing rate adopted in
this chapter.

For a suitable number $k$, e.g. $k \ge 20$, and provided that
$T_0$ is sufficiently large, a fair approximation of the 95 \%
confidence interval for the value $\nu^+(\xi)$ can be obtained as
$\mbox{CI}_{0.95}(\xi)= \Big( C^-(\xi)\,, \: C^+(\xi) \Big)$,
where
\begin{equation}\label{eq:cb}
C^\pm(\xi) = \hat{\nu}^+(\xi) \pm 1.96 \,
\frac{\hat{s}(\xi)}{\sqrt{k}} \,,
\end{equation}
and the empirical standard deviation $\hat{s}(\xi)$ is given as
\begin{align}\label{eq:emcr4}
\hat{s}(\xi)^2 = \frac{1}{k-1} \sum_{j=1}^k  \, \Big(
\frac{n_j^+(\xi;0,T_0)}{T_0} - \hat{\nu}^+(\xi) \,
\Big)^2\,.
\end{align}

Note that $k$ and $T_0$ may not necessarily be the number and
length of the actually simulated response time series. Rather,
they may be chosen to optimize the estimate of
Eq.~(\ref{eq:emcr4}). If initially, $\tilde{k}$ time series of
length $\tilde{T}$ are simulated, then $k=\tilde{k}k_0$ and
$\tilde{T}=k_0 T_0$. That is, each initial time series of length
$\tilde{T}$ has been divided into $k_0$ time series of length
$T_0$, assuming, of course, that $\tilde{T}$ is large enough to
allow for this in an acceptable way. The consistency of the
estimates obtained by Eq.~(\ref{eq:emcr4}) can be checked for
large values of $\xi$ by the observation that
$\textrm{Var}[N^+(\xi; 0,t)] = \nu^+(\xi)t$ since $N^+(\xi; 0,t)$
is then a Poisson random variable by assumption. This leads to the
equation
\begin{align}\label{eq:P}
\hat{s}(\xi)^2 =
\frac{1}{k}\textrm{Var}\left[\sum_{j=1}^k\frac{N_j^+(\xi;
	0,T_0)}{T_0}\right]=\frac{\nu^+(\xi)}{T_0} \,,
\end{align}
where $\{N^+_1(\xi;0,T_0),...,N^+_k(\xi;0,T_0)\}$ denotes a random
sample with a possible outcome
$\{n^+_1(\xi;0,T_0),...,n^+_k(\xi;0,T_0)\}$. Hence,
$\hat{s}(\xi)^2/k\approx\nu^+(\xi)/kT_0$. Since this last relation
is consistent with the adopted assumptions, it could have been
used as the empirical estimate of the sample variance in the first
place. It is also insensitive to the blocking of data discussed
above since $kT_0=\tilde{k}\tilde{T}$. However, the advantage of
Eq.~(\ref{eq:emcr4}) is that it applies whatever the value of
$\xi$, and it does not rely on any specific assumptions about the
statistical distributions involved.





\section{Numerical Examples for Gaussian Random Fields}\index{Gaussian random field!numerical examples}


Since a major motivation for developing the approach discussed in this chapter is
application to random ocean wave fields, the numerical examples presented will
be for zero mean, homogeneous Gaussian fields specified by spectral densities defined in terms of a so called JONSWAP spectrum
$S(\omega)$, which is a one-sided spectrum, and a directional spreading function $D(\theta)$ \citep{sais81}:
\begin{equation}\label{eq:ne1}
S(\omega)=\frac{\alpha g^2}{\omega^5} \exp\Bigg\{ - \frac{5}{4}
\Big(\frac{\omega_p}{\omega}\Big)^4 + \ln \gamma \, \exp\Big[ -
\frac{1}{2 \sigma^2} \Big(\frac{\omega}{\omega_p} - 1 \Big)^2
\Big] \Bigg\}, \: \omega > 0,
\end{equation}
where $g = 9.81$ m\,sec$^{-2}$, $\omega_p$ denotes the peak frequency
in rad/sec and $\alpha$,
$\gamma$ and $\sigma$ are parameters related to the spectral shape.
$\sigma = 0.07$ when $\omega \le \omega_p$, and $\sigma = 0.09$ when
$\omega > \omega_p$. The parameter $\gamma$ is chosen to be equal
to 3.0. The parameter $\alpha$ is determined from the following
empirical relationship \citep{nagate07aor},
%
\begin{equation}\label{eq:ne2}
\alpha = 5.06 \Big(\frac{H_s}{T_p^2}\Big)^{2} \big( 1 - 0.287 \,
\ln \gamma \, \big)
\end{equation}
%
$H_s$ = significant wave height and $T_p= 2 \pi/ \omega_p$ =
spectral peak wave period. $H_s$ = 14.0 m and $T_p$ = 16 sec are
chosen for all the following examples. The directional spreading function is expressed as,
%
\begin{equation}\label{eq:ne3}
D(\theta) = {2^{2s-1}\Gamma^2(s+1)\over\pi\Gamma(2s+1)}
\cos^{2s}{(\theta-\theta_0)\over 2} \,,
\end{equation}
%
where $\theta_0$ is the main wave direction, chosen so that
$\theta_0 = 0$, and then $- \pi < \theta \le \pi$. The choice of
spreading parameter $s = 8$, which is a typical value often used.

The extreme value distributions PA and PB may be considered as
having been derived under the Poisson assumption, i.e
when independence of individual upcrossings or exceedances is assumed,
supplemented by Rice-like formulas for the exceedance rates
$\nu^+_A(\xi)$ and $\nu^+_B(\xi)$. Thus PA and PB are liable to
work in case of a wide-band spectrum $\Psi(\mathbf{k}, \omega)$,
as they do in the one-dimensional case of a random process. 
It is therefore to be expected that for the narrow-band case, 
there will be some discrepancy between the exact values and the predictions 
provided by PA and PB. Also the simplified approach proposed in this chapter is based on the use of the mean upcrossing rate function. However, this approach can easily be amended to cope with statistical dependence between the peak values in the extracted time series by using the ACER method, cf. Chapter~\ref{chap:acer}.  

Another aspect that comes into question in the case of a random ocean wave field, is the effective dimension of the field. Due to the dispersion relation, which implies a strong coupling between frequency and wave number, the time and space dimensions of the wave field will likewise be strongly coupled, indicating a degeneracy. In fact, the dispersion relation implies that all the mass of the spectral density function $\Psi(\mathbf{k}, \omega)$ will be localized on a surface in $(\omega, {\bf k})$-space.

Numerical results obtained for three examples will now be presented for $d \le 2$. The simulations of the example random fields were carried out on a standard desktop computer, and the required CPU time was of the order of one hour for each case. For each of the examples, predictions of the 99.9\% quantile of a three hour extreme value distribution by the proposed method (marked by an asterix in the relevant figures) and the analytical method are also provided.






\subsection{1+1-dimensional Gaussian field}\index{Gaussian random field!1+1-dimensional Gaussian field}

In the first example, a wide-band 1+1-dimensional
Gaussian field is considered. It has a spectral density given as,
%
%
\begin{equation}
\Psi(\omega,k)=S(|\omega|)S(|k|)\,,
\end{equation}
%
where $S(\cdot)$ is defined by Eq.~\eqref{eq:ne1}. This is clearly an artificial example, but it serves the purpose of securing no interaction between the time and space dimensions of the random field.
It is therefore to be expected that the theoretical formulas with $d = 1$ should apply. A plot of the spectral density is shown in Figure~\ref{fig:stf1}, while part of a realization of the associated Gaussian random field is given in Figure~\ref{fig:stf2}. Plots of the empirical and analytical results are presented in Figures~\ref{fig:stf3} and \ref{fig:stf4} for one single point ($L_1 = 0$), which corresponds to an ordinary Gaussian temporal process, and for an interval ($L_1 = 100$), respectively. It is seen that the empirical and analytical results are in very good agreement in the tail for both cases, and that the predictions of the 99.9\% quantile in the three hour extreme value distributions are in good agreement. 
%

\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.5in]{Figures/1psd.pdf}
	\end{center}
	\caption{The spectral density $\Psi(\omega,k)$ of the 1+1-dimensional Gaussian field.} \label{fig:stf1}
\end{figure}


\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.5in]{Figures/1surf.pdf}
	\end{center}
	\caption{Part of a realization of the 1+1-dimensional Gaussian random field.} \label{fig:stf2}
\end{figure}



\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.5in]{Figures/1_x=0.pdf}
	\end{center}
	\caption{Plot of the empirical $\log_{10}(\nu^+(\xi))$ ($\bullet \, \bullet \, \bullet$), analytical $\log_{10}(\nu^+_B(\xi))$ (thick solid line), optimal fitted curve ($- \, \cdot \, -$) with both empirical (dotted lines) and fitted (dashed lines) 95\% confidence band for a point $L_1 = 0$ for Example~1. Predictions of 99.9\% quantile: 5.22 (analytical); 5.06 and 95\% CI = (4.93, 5.18) (proposed method).} \label{fig:stf3}
\end{figure}


\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.5in]{Figures/1_x=100.pdf}
	\end{center}
	\caption{Plot of the empirical $\log_{10}(\nu^+(\xi))$ ($\bullet \, \bullet \, \bullet$), analytical $\log_{10}(\nu^+_B(\xi))$ (thick solid line), optimal fitted curve ($- \, \cdot \, -$) with both empirical (dotted lines) and fitted (dashed lines) 95\% confidence band for an interval $L_1 = 100$\,m for Example~1. Predictions of 99.9\% quantile: 6.09 (analytical); 6.11 and 95\% CI = (6.03, 6.18) (proposed method).} \label{fig:stf4}\vspace{4mm}
\end{figure}



\vspace{3mm}


\subsection{1+1-dimensional Gaussian sea}\index{Gaussian random field!1+1-dimensional Gaussian sea}

The second example is a 1+1-dimensional Gaussian sea, which corresponds to the limiting case of long-crested waves. In this case the power spectral density assumes the form,
%
\begin{equation}
\Psi(\omega,k)={1\over4}S(|\omega|) \, \delta\left(k-{\omega|\omega|\over
	g}\right) \,,
\end{equation}
%
where the dispersion relation for deep water waves is implemented in the form of a delta function. The resulting degeneracy would expectedly have some influence on the extreme value distribution. In the present case, the wave field can be represented as,
%
\begin{equation}\label{eq:st140}
X (t,x)={\rm Re} \Big\{ \sum_{j=0}^{N}\sqrt{S(\omega_j)\Delta\omega} \cdot
C_j \  e^{\ii\omega_jt - \ii {\omega_j^2\over g}x} \Big\},
\end{equation}
%
where $C_j$ are complex $N(0,1)$-distributed random variables, that is, $C_j = R_j + \ii S_j$ with $R_j$ and $S_j$ two independent $N(0, 1/2)$-distributed variables. 

Part of a realization of the Gaussian random wave field generated by Eq.~\eqref{eq:st140}, is plotted in Figure~\ref{fig:stf5}. It is clearly seen how the strong coupling between frequency and wave number manifests itself. 


\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.5in]{Figures/4surf.pdf}
	\end{center}
	\caption{Part of a realization of the long-crested random wave field.} \label{fig:stf5}
\end{figure}


Plots of the empirical and analytical results are presented in Figures~\ref{fig:stf6} and \ref{fig:stf7} for a point ($L_1 = 0$), which again corresponds to an ordinary Gaussian temporal process, and for an interval ($L_1 = 100$\,m). It is seen that the empirical and analytical results are still in complete agreement in the tail for the one-point case, while there is a significant discrepancy between the results for $L_1 = 100$\,m. In fact, the results for $L_1 = 100$\,m are only slightly higher than the corresponding results for the one-point case, illustrating the effect of the degeneracy on the extreme values. It is interesting to observe that if the normalization constant $\ell_1^0$ is calibrated to the empirical results, very accurate predictions are obtained, and they remain accurate for other domains without recalibration. This will be illustrated in the next example. 


\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.5in]{Figures/4_x=0.pdf}
	\end{center}
	\caption{Plot of the empirical $\log_{10}(\nu^+(\xi))$ ($\bullet\, \bullet \, \bullet$), analytical $\log_{10}(\nu^+_B(\xi))$ (thick solid line), optimal fitted curve ($- \, \cdot \, -$) with both empirical (dotted lines) and fitted (dashed lines) 95\% confidence band for a point $L_1 = 0$ for Example~2. Predictions of 99.9\% quantile: 5.22 (analytical); 5.27 and 95\% CI = (5.14, 5.41) (proposed method).} \label{fig:stf6}
\end{figure}


\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.5in]{Figures/4_x=100.pdf}
	\end{center}
	\caption{Plot of the empirical $\log_{10}(\nu^+(\xi))$ ($\bullet\, \bullet \, \bullet$), analytical $\log_{10}(\nu^+_B(\xi))$ (thick solid line), optimal fitted curve ($- \, \cdot \, -$) with both empirical (dotted lines) and fitted (dashed lines) 95\% confidence band for an interval $L_1 = 100\,\mbox{m}$ for Example~2. Predictions of 99.9\% quantile: 5.70 (analytical); 5.47 and 95\% CI = (5.36, 5.58) (proposed method).} \label{fig:stf7}\vspace{4mm}
\end{figure}





\subsection{A short-crested Gaussian sea}\index{short-crested Gaussian sea}\index{Gaussian random field!short-crested Gaussian sea}


In this example a model of a short-crested Gaussian random sea way is simulated. It corresponds to a case where the main direction of propagation of the waves is the $x$-direction. Having assumed deep water, the relation between the wave number
vector and the spreading angle will be $\mathbf{k} = (k_x, k_y) =
(\omega^2 /g) \big( \cos \theta, \sin \theta  \big)$. The spectral
density of the random field $X$ can then be expressed as
\begin{equation}\label{Psi}
\Psi({\bf k},\omega) = \int{ S(\omega)\, D(\theta)\, \delta\!\!\left(k_x - {\omega^2\over g} \cos\theta\right)\, \delta\!\!\left(k_y - {\omega^2\over g} \sin\theta\right)} \, {\rm d} \theta\,.
\end{equation}
This leads to the expressions,
\begin{eqnarray}\label{n_txy}
&&\Lambda_{xx} = \int{ k_x^2\, \Psi({\bf k},\omega)\, {\rm d}\omega\, {\rm d}{\bf k} }
= D_c\, \omega_c^4\, \sigma^2 / g^2 \, , \label{Lxx} \nonumber \\
&&D_c= \int{  D(\theta)\, \cos^2(\theta)} \,{\rm d}\theta \,, \nonumber \\
&&\omega_c^4 = {1\over\sigma^2} \int{\omega^4 S(\omega) {\rm d}\omega} = (2\pi/t_c)^4\,, \nonumber \\
&&\Lambda_{yy} = \int{ k_y^2\, \Psi({\bf k},\omega)\, {\rm d}\omega\, {\rm d}{\bf k} }
= (1-D_c)\, \omega_c^4\, \sigma^2 / g^2 \,
, \label{Lyy} \nonumber \\
&&\Lambda_{tt}= \int{ \omega^2\, \Psi({\bf k},\omega)\, {\rm d}\omega\, {\rm d}{\bf k} }  = \omega_z^2 \sigma^2 \,, \label{Ltt} \\
&&\ell_1^0= 2\pi\sigma/\sqrt{\Lambda_{xx}} = {2\pi g\over \sqrt{D_c}\omega_c^2}
= {g\,  t_c^2\over 2\pi \sqrt{D_c}} \,
, \label{x_0} \nonumber \\
&&\ell_2^0= 2\pi\sigma/\sqrt{\Lambda_{xx}} = {2\pi g\over \sqrt{1-D_c}\omega_c^2}
= {g\,  t_c^2\over 2\pi \sqrt{1-D_c}}
, \label{y_0} \nonumber \\
&&\ell_3^0= 2\pi\sigma/\sqrt{\Lambda_{tt}} = 2\pi/\omega_z = T_z
, \label{t_0} \nonumber \\  
&& q_x = L_1/\ell_1^0,\quad q_y = L_2/\ell_2^0 \, \quad q_t = L_3/\ell_3^0 = T/T_z \,, \nonumber
\end{eqnarray}
where $(L_1, L_2, T)$ is the size of the rectangular domain of interest.


Calculating the values of the parameters using the formulas above gives the values $\ell_1^0=135$\,m and $\ell_2^0=280$\,m.
Plotting the empirical against the analytical results obtained for the calculated parameter values for $L_1 = L_2 = 0$,  $L_1 \times L_2 = 100\,\mbox{m} \times 0$, $L_1 \times L_2 = 0 \times 100\,\mbox{m}$ and $L_1 \times L_2 = 100\,\mbox{m} \times 100\,\mbox{m}$ leads to Figures~\ref{fig:stf8} - \ref{fig:stf11}. 


\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.5in]{Figures/Sea_1_XY=0x0.pdf}
	\end{center}
	\caption{Plot of the empirical $\log_{10}(\nu^+(\xi))$ ($\bullet\, \bullet \, \bullet$), analytical $\log_{10}(\nu^+_B(\xi))$ (thick solid line), optimal fitted curve ($- \, \cdot \, -$) with both empirical (dotted lines) and fitted (dashed lines) 95\% confidence band for $L_1 = L_2 = 0$ for Example~3. Predictions of 99.9\% quantile: 5.22 (analytical); 5.13 and 95\% CI = (5.01, 5.25) (proposed method).} \label{fig:stf8}
\end{figure}


\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.5in]{Figures/Sea_1_XY=100x0.pdf}
	\end{center}
	\caption{Plot of the empirical $\log_{10}(\nu^+(\xi))$ ($\bullet\, \bullet \, \bullet$), analytical $\log_{10}(\nu^+_B(\xi))$ (original: thick dashed line; tuned: thick solid line), optimal fitted curve ($- \, \cdot \, -$) with both empirical (dotted lines) and fitted (dashed lines) 95\% confidence band for the interval $L_1 \times L_2 = 100\,\mbox{m} \times 0$ for Example~3. Predictions of 99.9\% quantile: 5.47 (tuned analytical); 5.33 and 95\% CI = (5.26, 5.40) (proposed method).} \label{fig:stf9}
\end{figure}


\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.5in]{Figures/Sea_1_XY=0x100.pdf}
	\end{center}
	\caption{Plot of the empirical $\log_{10}(\nu^+(\xi))$ ($\bullet\, \bullet \, \bullet$), analytical $\log_{10}(\nu^+_B(\xi))$ (original: thick dashed line; tuned: thick solid line), optimal fitted curve ($- \, \cdot \, -$) with both empirical (dotted lines) and fitted (dashed lines) 95\% confidence band for the interval $L_1 \times L_2 = 0 \times 100\,\mbox{m}$ for Example~3. Predictions of 99.9\% quantile: 5.47 (tuned analytical); 5.31 and 95\% CI = (5.18, 5.41) (proposed method).} \label{fig:stf10}
\end{figure}


\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.5in]{Figures/Sea_1_XY=100x100.pdf}
	\end{center}
	\caption{Plot of the empirical $\log_{10}(\nu^+(\xi))$ ($\bullet\, \bullet \, \bullet$), analytical $\log_{10}(\nu^+_B(\xi))$ (original: thick dashed line; tuned: thick solid line), optimal fitted curve ($- \, \cdot \, -$) with both empirical (dotted lines) and fitted (dashed lines) 95\% confidence band for the square domain $L_1 \times L_2 = 100\,\mbox{m} \times 100\,\mbox{m}$ for Example~3. Predictions of 99.9\% quantile: 5.72 (tuned analytical); 5.74 and 95\% CI = (5.56, 5.89) (proposed method).} \label{fig:stf11}
\end{figure}

It is seen from these figures that there is also in this case a significant effect of dimensional degeneracy. Since the main wave direction is along the $x$-axis, the strongest effect of dimensional degeneracy is to be expected for the cases  $L_1 \times L_2 = 100\,\mbox{m} \times 0$ and $L_1 \times L_2 = 100\,\mbox{m} \times 100\,\mbox{m}$, which is also corroborated by the plots. 

As mentioned in the previous example, by tuning the parameters $\ell_1^0$ and $\ell_2^0$ to the empirical results, the tuned theoretical predictions become very accurate. The tuning can be done for each dimension separately. That is, $\ell_1^0$ is calibrated to the results for the case $L_1 \times L_2 = 100\,\mbox{m} \times 0$, while $\ell_2^0$ is calibrated to the results for the case $L_1 \times L_2 = 0 \times 100\,\mbox{m}$. It has been verified that the calibrated values apply also to other domain sizes. For the particular example at hand, it was found that the tuned parameters are $\ell_1^0 = 474$\,m, $\ell_2^0 = 477$\,m. The results obtained by using the tuned parameters have also been plotted in Figures~\ref{fig:stf8} - \ref{fig:stf11} together with the predictions of the 99.9\% quantiles. It is seen that fairly good agreement is now achieved for all cases. 




It has been mentioned already that the tuned version of PB seems to give accurate predictions for practically any size. This has been verified for a range of domain sizes providing substantial support for this assertion. It is also of interest to investigate to what extent PA can give accurate predictions. If tuned to each particular size, it could perhaps provide reasonable predictions, but that is impractical. Let us therefore illustrate how it performs using the same tuned parameters as previously established, that is, $\ell_1^0 = 474$\,m, $\ell_2^0 = 477$\,m. The results for the two square domains of size $100\,\mbox{m} \times 100\,\mbox{m}$ and $300\,\mbox{m} \times 300\,\mbox{m}$ have been plotted in Figures~\ref{fig:stf16} and \ref{fig:stf17}. It is seen that there is a significant discrepancy between the tuned PA results and the empirical results for the smaller domain, while there is good agreement for the large one, giving a good indication of the importance of the neglected terms in PA. If the domain is too small, the leading term is not sufficient to obtain accurate estimates.  



\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.5in]{Figures/Sea_1P_XY=100x100.pdf}
	\end{center}
	\caption{Plot of the empirical $\log_{10}(\nu^+(\xi))$ ($\bullet\, \bullet \, \bullet$), analytical $\log_{10}(\nu^+_A(\xi))$ (original: thick dash dot line; tuned: thick solid line) results for the square domain $L_1 \times L_2 = 100\,\mbox{m} \times 100\,\mbox{m}$ for Example 3. Empirical 95\% confidence band (thin dashed lines).} \label{fig:stf16}
\end{figure}


\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.5in]{Figures/Sea_1P_XY=300x300.pdf}
	\end{center}
	\caption{Plot of the empirical $\log_{10}(\nu^+(\xi))$ ($\bullet\, \bullet \, \bullet$), analytical $\log_{10}(\nu^+_A(\xi))$ (original: thick dash dot line; tuned: thick solid line) results for the square domain $L_1 \times L_2 = 300\,\mbox{m} \times 300\,\mbox{m}$ for Example 3. Empirical 95\% confidence band (thin dashed lines).} \label{fig:stf17}
\end{figure}





The value of the 99\% fractile of a three hour extreme value distribution, which corresponds to $\nu^+ = 10^{-6}$, as a function of the size of a square domain has been plotted in Figure~\ref{fig:stf20}. It is seen that the extreme values show a fairly strong dependence on the area of the domain. For prediction of extremes over e.g. the deck area of an offshore structure, which would typically be about $100\,\mbox{m} \times 100\,\mbox{m}$, it is seen that the area effect amounts to 10-15\% larger extremes than predicted by a one-point estimate.

\begin{figure}[!hbt]
	\begin{center}
		\includegraphics[width=2.2in]{Figures/Sea_1t_frac(-6)_X.pdf}
	\end{center}
	\caption{The 99\% fractile value of a three hour extreme value distribution as a function of the size of a square domain with $x = L_1 = L_2$ for the short-crested Gaussian sea.} \label{fig:stf20}
\end{figure}







\section{Numerical Examples for Non-Gaussian Random Fields}\label{sec:ngf}\index{non-Gaussian random field!numerical examples}


\subsection{A second-order wave field}\index{non-Gaussian random field!second-order wave field}

Since a major motivation for developing the approach discussed in this chapter is
application to random ocean wave fields, the numerical example presented will
be for a zero mean, homogeneous second-order short-crested random wave field. It has been shown, cf. e.g. \citep{tobionba08}, that such a wave field
can be built up in the following way. First, let $\eta(t,r,\theta)$ denote the wave elevation of a long-crested sea propagating in the direction specified by the angle $\theta$, where $r = x \cos \theta + y \sin \theta$. Then $\eta(t,r,\theta) = \eta_1(t,r,\theta) + \eta_2(t,r,\theta)$. Here the linear, first-order part $\eta_1(t,r,\theta)$ is given by the relation,
 \begin{equation}
  \eta_1(t,r,\theta)= \sum_{k=1}^{N} |C_k| \cos(\chi_k),
\end{equation}
where $C_k=|C_k|e^{\ii \varepsilon_k}$, $k = 1, \ldots, N$, is a set of independent complex random variables, whose real and imaginary parts are independent and normally distributed as $N\left(0, {S(\omega_k)\Delta\omega}\right)$. $\chi_k = r_k - \omega_k t + \varepsilon_k$, where $r_k = \omega^2_k r / g$. $S(\omega)$ is the JONSWAP spectrum, which was introduced in the previous section. The non-linear, second-order part $\eta_2(t,r)$ has the following representation,
\begin{align}
  \eta_2(t,&r,\theta) =  
  \sum_{k,l=1}^{N} \Big( m_{kl} |C_k C_l| \cos(\chi_k) \cos(\chi_l) \nonumber \\ &- M_{kl} |C_k C_l| \sin(\chi_k) \sin(\chi_l)\Big),
\end{align}
where $m_{kl}= (2g)^{-1} \min(\omega^2_k, \omega^2_l)$ and $M_{kl}= (2g)^{-1} \max(\omega^2_k, \omega^2_l)$.

The short-crested wave field $X(x,y,t)$ is then obtained by superposition of $m$ long-crested wave fields as follows,
\begin{align}
X(x,y,t) = \sum_{k=1}^m  \sqrt{D(\theta_k) \Delta \theta} 
\, \eta(t, x\cos \theta_k + y\sin \theta_k, \theta_k)
\end{align}
where $D(\theta)$ denotes the directional spreading function.


 

In the numerical example, a specific model of a short-crested, second-order random sea way is simulated. It corresponds to a case where the main direction of propagation of the waves is the $x$-direction. Having assumed deep water, the relation between the wave number
vector and the spreading angle will be $\mathbf{k} = (k_x, k_y) =
(\omega^2 /g) \big( \cos \theta, \sin \theta \big)$. The spectral
density in terms of wave numbers and frequency of the first-order part of the random field $X$ can then be expressed as
\begin{align}\label{Psi2}
&\Psi({\bf k},\omega) = \int_{- \pi}^\pi S(\omega)\, D(\theta)\, \, \nonumber \\ &\cdot \delta \left(k_x - \frac{\omega^2}{g} \cos\theta\right) \delta \left( k_y - \frac{\omega^2}{g} \sin\theta \right) \, {\rm d} \theta\,.
\end{align}



\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.5in]{Figures/2+1surf(3).pdf}
\end{center}
\caption{Part of a realization of the short-crested, second-order random wave field.} \label{fig:ngf1}
\end{figure}


Part of a realization of the short-crested, second-order random wave field is shown in Fig.~\ref{fig:ngf1}. The standard deviation of the wave field was estimated to be $\sigma = 3.53$m. From the simulated wave field, the four upcrossing rate functions of Eq.~\eqref{eq:st35} have been estimated by the linear regression approach described in Section~\ref{sec:stengrf} on space-time extremes. The results are shown in Figures~\ref{fig:ngf2} - \ref{fig:ngf5} together with the optimally fitted curves obtained by the point process procedure.


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.5in]{Figures/nu_0.pdf}
\end{center}
\caption{The one-point upcrossing rate $\nu_0^+(\xi)$, cf. Eq.~\eqref{eq:st35}.} \label{fig:ngf2}
\end{figure}

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.5in]{Figures/nu_x.pdf}
\end{center}
\caption{The upcrossing rate $\nu_x^+(\xi)$, cf. Eq.~\eqref{eq:st35}.} \label{fig:ngf3}
\end{figure}

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.5in]{Figures/nu_y.pdf}
\end{center}
\caption{The upcrossing rate $\nu_y^+(\xi)$, cf. Eq.~\eqref{eq:st35}.} \label{fig:ngf4}
\end{figure}


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.5in]{Figures/nu_xy.pdf}
\end{center}
\caption{The upcrossing rate $\nu_{xy}^+(\xi)$, cf. Eq.~\eqref{eq:st35}.} \label{fig:ngf5}
\end{figure}

The prediction results obtained by the simplified empirical approach in combination with the point process procedure have been plotted in Figures~\ref{fig:ngf6} - \ref{fig:ngf9}. Similarly, the results obtained by the parametric approach of Eq.~\eqref{eq:st35} are also plotted. In the figures the predicted 99.9\% fractile value of the three hour extreme value distribution obtained by the simplified approach has been indicated by an asterisk. This fractile value is achieved at $\nu^+ = 10^{-7}$. It is seen that the agreement between the two approaches proposed in this chapter is very good for the cases studied.


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.3in]{Figures/sea_0x0.pdf}
\end{center}
\caption{Case 1: Single point. Plot of the empirical upcrossing rate $\log_{10}(\nu^+(\xi))$ ($\bullet\, \bullet \, \bullet$), of the optimal fitted curve for the simplified approach ($- - -$), of $\log_{10}(\nu^+_P(\xi))$ (solid line) for the parametric approach. Both empirical (dotted lines) and fitted (dash-dotted lines) 95\% confidence bands are shown.  Predictions of 99.9\% quantile: 5.80 and 95\% CI = (5.75, 5.84) (simplified method); 5.77 (calibrated parametric).} \label{fig:ngf6}
\end{figure}


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.3in]{Figures/sea_99x99.pdf}
\end{center}
\caption{Case 2: $100 \times 100$m$^2$. Plot of the empirical upcrossing rate $\log_{10}(\nu^+(\xi))$ ($\bullet\, \bullet \, \bullet$), of the optimal fitted curve for the simplified approach ($- - -$), of $\log_{10}(\nu^+_P(\xi))$ (solid line) for the parametric approach. Both empirical (dotted lines) and fitted (dash-dotted lines) 95\% confidence bands are shown.  Predictions of 99.9\% quantile: 6.98 and 95\% CI = (6.73, 7.23) (simplified method); 6.89 (calibrated parametric).} \label{fig:ngf7}
\end{figure}



\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.5in]{Figures/sea_300x300.pdf}
\end{center}
\caption{Case 3: $300 \times 300$m$^2$. Plot of the empirical upcrossing rate $\log_{10}(\nu^+(\xi))$ ($\bullet\, \bullet \, \bullet$), of the optimal fitted curve for the simplified approach ($- - -$), of $\log_{10}(\nu^+_P(\xi))$ (solid line) for the parametric approach. Both empirical (dotted lines) and fitted (dash-dotted lines) 95\% confidence bands are shown.  Predictions of 99.9\% quantile: 7.35 and 95\% CI = (7.11, 7.62) (simplified method); 7.45 (calibrated parametric).} \label{fig:ngf8}
\end{figure}



\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.5in]{Figures/sea_501x501.pdf}
\end{center}
\caption{Case 4: $500 \times 500$m$^2$. Plot of the empirical upcrossing rate $\log_{10}(\nu^+(\xi))$ ($\bullet\, \bullet \, \bullet $), of the optimal fitted curve for the simplified approach ($- - -$), of $\log_{10}(\nu^+_P(\xi))$ (solid line) for the parametric approach. Both empirical (dotted lines) and fitted (dash-dotted lines) 95\% confidence bands are shown.  Predictions of 99.9\% quantile: 7.79 and 95\% CI = (7.50, 8.10) (simplified method); 7.72 (calibrated parametric).} \label{fig:ngf9}
\end{figure}


The value of the 99.9\% fractile of a three hour extreme value distribution as a function of the size of a square domain has been plotted in Figure~\ref{fig:ngf10}. It is seen that the extreme values show a fairly strong dependence on the area of the domain. For prediction of extremes over e.g. the deck area of an offshore structure, which would typically be about $100\,\mbox{m} \times 100\,\mbox{m}$, it is seen that the area effect may amount to as much as 15-20\% larger extreme crest heights than predicted by one-point estimates.

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.5in]{Figures/fractile(1e-7).pdf}
\end{center}
\caption{The 99.9\% fractile value of a three hour extreme value distribution as a function of the size of a square domain with $x = L_x = L_y$.} \label{fig:ngf10}
\end{figure}



\subsection{A Student's $t$ random field}\index{non-Gaussian random field!Student's $t$ random field}

To test the proposed methods on a strongly non-Gaussian random field, what shall be referred to as a Student's $t$ field has been constructed. This was done in the following manner. Let $X(x,t) = \eta_1(t,x,0)$, where $\eta_1(t,x,0)$ is as defined in the previous example. $X(x,t)$ is then a homogeneous Gaussian random field of zero mean with one spatial dimension. Let $X_j(x,t)$, $j = 1, \ldots, 4$ denote independent copies of $X(x,t)$. The random field $Z(x,t)$ is now constructed as follows,
\begin{equation}
Z(x,t) = \frac{X(x,t)}{\sqrt{\frac{1}{4}\sum_{k=1}^4 X_j(x,t)^2} + \frac{\sigma}{20} } \,,
\end{equation}
where $\sigma^2$ denotes the variance of $X(x,t)$. The added term $\sigma / 20$ in the denominator is introduced to avoid near singularities in the generated fields and make sure that it satisfies the conditions discussed after Eq.~\eqref{eq:st14}. Part of a realization of the resulting random field is shown in Figure~\ref{fig:ngf11}.

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.5in]{Figures/1+1surf.pdf}
\end{center}
\caption{Part of a realization of the Student's $t$ random field.} \label{fig:ngf11}
\end{figure}

In Figures~\ref{fig:ngf12} and \ref{fig:ngf13} are shown the obtained results for $\nu_0^+(\xi)$ and $\nu_x^+(\xi)$, respectively. These figures clearly display the strongly non-Gaussian characteristics of this random field.
 
\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.5in]{Figures/1nu_0.pdf}
\end{center}
\caption{The one-point upcrossing rate $\nu_0^+(\xi)$, cf. Eq.~\eqref{eq:st35}.} \label{fig:ngf12}
\end{figure} 

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.5in]{Figures/1nu_x.pdf}
\end{center}
\caption{The upcrossing rate $\nu_x^+(\xi)$, cf. Eq.~\eqref{eq:st35}.} \label{fig:ngf13}
\end{figure} 
 

Figure~\ref{fig:ngf14} shows the one-point prediction result, while Figure~\ref{fig:ngf15} shows the prediction results for an interval of size 100~m. It is seen that there is a strong extensional effect of about 50\% on the predicted extreme values . It is also clearly demonstrated that the predictions obtained by the two proposed methods are in excellent agreement. 

\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.5in]{Figures/1+1_0.pdf}
\end{center}
\caption{Single point. Plot of the empirical upcrossing rate $\log_{10}(\nu^+(\xi))$ ($\bullet\, \bullet \, \bullet$), of the optimal fitted curve for the simplified approach ($- - -$), of $\log_{10}(\nu^+_P(\xi))$ (solid line) for the parametric approach. Both empirical (dotted lines) and fitted (dash-dotted lines) 95\% confidence bands are shown.  Predictions of 99.9\% quantile: 32.33 and 95\% CI = (30.25, 34.69) (simplified method); 32.33 (calibrated parametric).} \label{fig:ngf14}
\end{figure}


\begin{figure}[!hbt]
\begin{center}
\includegraphics[width=2.5in]{Figures/1+1_100.pdf}
\end{center}
\caption{$x = 100$m. Plot of the empirical upcrossing rate $\log_{10}(\nu^+(\xi))$ ($\bullet\, \bullet \, \bullet$), of the optimal fitted curve for the simplified approach ($- - -$), of $\log_{10}(\nu^+_P(\xi))$ (solid line) for the parametric approach. Both empirical (dotted lines) and fitted (dash-dotted lines) 95\% confidence bands are shown.  Predictions of 99.9\% quantile: 48.92 and 95\% CI = (47.73, 49.78) (simplified method); 47.70 (calibrated parametric).} \label{fig:ngf15}
\end{figure}





\bigskip

\section{Comments}


Analytical formulas for the extreme value distribution of a homogeneous Gaussian random field has been discussed at some length, and some properties of these formulas have been highlighted. It has been shown that by proper calibration, the analytical formulas provide good results also for the case of a Gaussian random sea.

Provided there is no degeneracy in the dimensionality of the Gaussian random field, the extreme value distribution PA is asymptotically accurate with respect to the level $\xi$ as well as to the size of the domain, while PB is asymptotically accurate with respect to the level $\xi$ for any size. 

The region of applicability of the approximation PB is suggested to be $\nu^+_B(\xi) \le 10^{-2}$. It is then formulated in terms of the value of $\nu^+_B$ rather then for the level $\xi$, and it appears to be valid for a domain of arbitrary size $L_1\times L_2$. 

The accuracy of the proposed prediction procedure obtained by optimal fitting to the simulated data has also been amply demonstrated. This leaves the door open for accurate prediction of extremes of non-Gaussian random fields. While the analytical formulas are restricted to homogeneous Gaussian fields, no such restriction applies to the proposed prediction method.


A simplified, empirically based procedure for prediction of space-time extremes of homogeneous non-Gaussian random fields has been presented. The advantage of this method is its flexibility and general applicability, while it appears to be accurate and robust. A drawback of the method is that if the spatial domain is changed, the whole analysis has to be repeated. In an effort to ameliorate this situation, a quasi-parametric representation of the extreme value distribution for rectangular domains is proposed. When this approach applies, it has the advantage that it only needs to be calibrated for one rectangular domain. After that, it applies to any rectangular domain (with the same orientation). The case studies presented show that the two methods proposed seem to provide very good practical tools for the prediction of space-time extremes of random fields.     








% ------------------------------------------------------------------------
