Multipoint reconstruction of wind speeds

The most intermittent behaviour of atmospheric turbulence is found for very short timescales. Based on a concatenation of conditional probability density functions (cpdf's) of nested wind speed increments, inspired by a Markov process in scale, we derive a short-time predictor for wind speed fluctuations around a non-stationary mean value and with a corresponding non-stationary variance. As a new quality this short-time predictor enables a multipoint reconstruction of wind data. The used cpdf's are (1) directly estimated from historical data from the offshore research platform FINO1 and (2) obtained from numerical solutions of a family of Fokker–Planck equations in the scale domain. The explicit forms of the Fokker–Planck equations are estimated from the given wind data. A good agreement between the statistics of the generated and measured synthetic wind speed fluctuations is found even on timescales below 1 s. This shows that our approach captures the short-time dynamics of real wind speed fluctuations very well. Our method is extended by taking the non-stationarity of the mean wind speed and its non-stationary variance into account.

1 Introduction

The transition of our energy system, formerly strongly relying on gas and coal, to a decarbonised one, mainly based on wind, solar and hydropower, is still ongoing work, but great progress has been made. From 2005 to 2017 the share in installed capacity of wind (solar) has been increased from 6 % (0.3 %) to 18 % (11.5 %) in the European Union (WindEurope, 2018). The downside of this increasing share of fluctuating renewable energy sources is their integration into the power grid. By analysing measurements of fed-in wind and solar power, it could be shown that their fluctuations strongly deviate from Gaussian behaviour on timescales ranging from hours to seconds (Anvari et al., 2016) and for wind power even for scales below 1 s (Haehne et al., 2018). This survival of the atmospheric intermittency in the power grid poses the grid operators with the great challenge to ensure stable power supply, even under highly volatile conditions. Within this context the term intermittency is used in the spirit of Kolmogorov (1962) to describe the characteristic heavy-tailed shape of pdf's often found at small scales in time series of turbulent systems (Frisch, 2004).

To aid the design of our future energy systems, for example to size the needed energy storage or the power generation capacity of conventional power plants, much work has been done in the field of long-term wind speed and power modelling, utilising Markov chain models. Whereas simple first-order Markov chain models cannot grasp the characteristics of long-term correlations of wind speeds (Brokish and Kirtley, 2009), higher-order Markov chain models perform better, but will require more input data for estimating the transition matrices or some simplifications (Pesch et al., 2015; Brokish and Kirtley, 2009; Papaefthymiou and Klockl, 2008; Weber et al., 2018).

Despite their dramatic effect of long-range correlation and fluctuations of wind speeds on the power generation (and thus the grid stability), wind speed fluctuations are known to be most intermittent on short timescales (Boettcher et al., 2003). With short timescales we refer to timescales in the range of seconds to minutes. As can be seen in Boettcher et al. (2003), the effect of intermittency is most prominent at timescales <1 s, but as the timescales increase, the pdf's broaden. Models considering time steps ranging from minutes to seconds or even below are of course not suited for energy system analysis on national levels, but they are useful tools for wind turbine operators. For a time step of 10 min Carpinone et al. (2010) presented a higher-order Markov chain model for wind power, and Milan et al. (2013) showed a stochastic power model based on a conditional Langevin equation to work even in the range of seconds.

The knowledge of full three-dimensional wind fields for all three velocity components and the pressure in all details would be desirable. The lack of basic understanding, the impracticability of handling such huge data sets and the complexity of the wind energy conversion process often lead to the demand for simplified models for wind speed. Common approaches for the design of wind turbines are the so-called Mann uniform shear and the Kaimal spectral and exponential coherence model (IEC, 2005). Both models take spectral and coherence aspects of turbulent velocity fluctuations into account, thus handling the fluctuations as Gaussian distributed and stationary. Higher-order statistical effects like the prominent intermittency effect of turbulence and non-stationarities are not taken into account; see for example Mücke et al. (2011). Another approach is to use one-dimensional effective wind speed time series, representing for example the wind field together with the rotor aerodynamics as it impacts the drivetrain or can be used to model the above-discussed energy conversion process (Wächter et al., 2011).

Within this work we propose a novel stochastic generator of one-dimensional wind speed fluctuations with a sampling interval of 0.1 s. One main novelty is that we show how to grasp by this model multipoint statistics of wind structures in time. While commonly applied methods, like spectral analysis and two-point correlations, limit themselves to two-point statistics, here we extend the methodology to more than two points in time. We obtain generalised correlations between multiple points in time, in terms of probability density functions (pdf's), for the occurrence of a whole sequence of wind speeds. Those pdf's we denote multipoint pdf's, and they constitute the basic concept of our approach. Such a stochastic multipoint approach should in principle be able to grasp wind structures like gusts as well as clustering of wind fluctuations. The method was initially developed by Nawroth and Peinke (2006) in the context of homogeneous isotropic turbulence and later on applied to the modelling of log-return rates of current exchange rates (Nawroth et al., 2010) and velocity increments of idealised homogeneous isotropic turbulence (Stresing and Peinke, 2010). The successful application to ocean gravity waves (Hadjihosseini et al., 2016) showed that structures of monster waves can be grasped by this approach correctly (Hadjihoseini et al., 2018). For a recent review see (Peinke et al., 2019). Finally we want to point out that we also show how to handle the aspect of non-stationary wind conditions.

We will continue as follows. In the first part we discuss the method for a subset of wind data characterised by its mean wind speed and its standard deviation. For such data it is shown how, arising from a Langevin process in scale, a predictor for the upcoming wind speed fluctuation around a mean value can be derived by a nesting of conditional probability density functions. Afterwards we check for Markovian properties of the wind speed fluctuations in scale and set up a Fokker–Planck equation, corresponding to the Langevin process in scale, and we show how it contributes to the improvement of our stochastic prediction method. Finally, in the second part, the non-stationary mean wind speed and its non-stationary variance are incorporated into our approach to achieve more realistic wind speed time series.

2 Method

In this section we present the stochastic framework used for our multipoint reconstruction scheme. As a simplification we start this discussion for blocks of wind data U(t), with the time t, which share the same mean wind speed $\stackrel{\mathrm{‾}}{U}$ and the same standard deviation σ_U, as suggested in Morales et al. (2012). Fixing $\stackrel{\mathrm{‾}}{U}$ and σ_U one would generate quasi-stationary subsets of data. We follow this idea, but instead of creating quasi-stationary subsets, we aim at creating quasi-stationary time series. With U(t) we refer to the resulting wind speed from the horizontal components. The quasi-stationary wind speed u(t) is then obtained from U(t) by normalising it with the mean $\stackrel{\mathrm{‾}}{U}$ and standard deviation σ_U within blocks of 1 min length. We use measured data from the offshore research platform FINO1: the data were recorded at a sampling frequency of 1 Hz between calendar weeks 1 and 10 in 2007 with an ultrasonic anemometer, mounted at 80 m height, resulting in approximately 6×10⁶ samples. The use of our method for non-stationary time series is outlined in Sect. 3.

2.1 Multipoint statistics

Since we assume wind speeds to emerge from a turbulent cascade, increments will play a key role in our method. Having a time series of wind speeds u(t), the corresponding increment time series Δu(τ), depending on a certain scale τ, is given by

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }u\left(\mathit{\tau }\right)=u\left(t\right)-u(t-\mathit{\tau }).\end{array}$$

This is the definition of so-called right-sided increments. Note that the calculation of Δu(τ) after this definition depends on the current wind speed u(t) and a past value u(t−τ), whereas the use of left-sided increments $\mathrm{\Delta }u\left(\mathit{\tau }\right)=u(t+\mathit{\tau })-u\left(t\right)$ would premise the knowledge of future values u(t+τ), creating a contradiction as we aim at producing synthetic wind speed time series. To ease readability we use the shorthand notations $u_i:=u(t-{\mathit{\tau }}_i)$ and $\mathrm{\Delta }{u}_i:=u\left(t\right)-u(t-{\mathit{\tau }}_i)$ with ${\mathit{\tau }}_i=i\cdot \mathit{\tau }$ ($i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots }$) in the following.

As a further remark we note that although we consider in this work time series of wind speed, we often talk of multipoint statistics. Δu(τ) is considered to be a statistical two-point quantity, which more correctly could be denoted as a two-time quantity. Based on the commonly used hypotheses on frozen turbulence by Taylor, for short-time fluctuations time and space statistics are related linearly by the mean wind speed (see also Peinke et al., 2019).

Our idea is to predict a wind speed $u^{*}\left(t^{*}\right)$ only by knowledge of its N past values $\mathit{\left\{}{u}_{\mathrm{1}}\right(t^{*}-{\mathit{\tau }}_{\mathrm{1}}),\mathrm{\dots },u_N(t^{*}-{\mathit{\tau }}_N\left)\mathit{\right\}}$. The probability of an event u^* to happen at time t^* can then be expressed by the conditional probability density function (cpdf)

$$\begin{array}{}\text{(2)}& \begin{array}{rl}p(& u^{*},t^{*}|u_{\mathrm{1}},t^{*}-{\mathit{\tau }}_{\mathrm{1}};\mathrm{\dots };u_N,t^{*}-{\mathit{\tau }}_N)=\\ & {\displaystyle \frac{p(u^{*},t^{*};u_{\mathrm{1}},t^{*}-{\mathit{\tau }}_{\mathrm{1}};\mathrm{\dots };u_N,t^{*}-{\mathit{\tau }}_N)}{p(u_{\mathrm{1}},t^{*}-{\mathit{\tau }}_{\mathrm{1}};\mathrm{\dots };u_N,t^{*}-{\mathit{\tau }}_N)}},\end{array}\end{array}$$

using the definition of conditional probabilities. Note this cpdf is the key quantity to estimate a new wind speed value u^*, and it can be used iteratively to generate new time series, as we show below.

Now we link Eq. (2) to the idea of an underlying turbulent cascade, which can be described by means of a Markov process. Thus we assume that there exists a scale separation $\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{ME}}={\mathit{\tau }}_j-{\mathit{\tau }}_i$ (j>i), after which the Markov condition is fulfilled for all greater scales. This scale separation Δτ_ME is often called the Markov–Einstein length (Einstein, 1905). It enables us to express the multiple cpdf $p(u^{*},t^{*}|u_{\mathrm{1}},t^{*}-{\mathit{\tau }}_{\mathrm{1}};\mathrm{\dots };u_N,t^{*}-{\mathit{\tau }}_N)$ by a multiplication of a simpler double cpdf and marginal pdf's:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}p(& u^{*},t^{*}|u_{\mathrm{1}},t^{*}-{\mathit{\tau }}_{\mathrm{1}};\mathrm{\dots };u_N,t^{*}-{\mathit{\tau }}_N)=\\ & {\displaystyle \frac{p\left(u^{*}\right)}{p\left(u_{\mathrm{1}}\right)}}\cdot {\displaystyle \frac{{\mathrm{\Pi }}_{i=\mathrm{1}}^{N-\mathrm{1}}p\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_{i+\mathrm{1}};u^{*})}{{\mathrm{\Pi }}_{i=\mathrm{2}}^{N-\mathrm{1}}p\left({\stackrel{\mathrm{̃}}{\mathrm{\Delta }u}}_i\right|{\stackrel{\mathrm{̃}}{\mathrm{\Delta }u}}_{i+\mathrm{1}};u_{\mathrm{1}})}}\cdot {\displaystyle \frac{p\left(\mathrm{\Delta }{u}_N\right|u^{*})}{p\left({\stackrel{\mathrm{̃}}{\mathrm{\Delta }u}}_N\right|u_{\mathrm{1}})}},\end{array}\end{array}$$

with $\stackrel{\mathrm{̃}}{\mathrm{\Delta }{u}_i}:=u(t^{*}-{\mathit{\tau }}_{\mathrm{1}})-u(t^{*}-{\mathit{\tau }}_i)$ with the timescale τ_i−τ₁. For details we refer the reader to Appendix A.

It is known that the evolution of cpdf's of a Markov process can be described by the famous Kramers–Moyal expansion (Risken, 1996), which notes for our stochastic process in scale Δu_i

$$\begin{array}{}\text{(4)}& \begin{array}{rl}-{\mathit{\tau }}_i& {\displaystyle \frac{\partial }{\partial {\mathit{\tau }}_i}}p\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_j;u^{*})=\\ & \sum _{n=\mathrm{1}}^{\mathrm{\infty }}{\left(-{\displaystyle \frac{\partial }{\partial \mathrm{\Delta }{u}_i}}\right)}^n\left[D^{\left(n\right)}(\mathrm{\Delta }{u}_i,{\mathit{\tau }}_i,u^{*})p\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_j;u^{*})\right],\end{array}\end{array}$$

requiring ${\mathit{\tau }}_j-{\mathit{\tau }}_i\ge \mathrm{\Delta }{\mathit{\tau }}_{\mathrm{ME}}$ and with the Kramers–Moyal coefficients D⁽ⁿ⁾ being defined as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{D}^{\left(n\right)}& (\mathrm{\Delta }{u}_i,{\mathit{\tau }}_i,u^{*})=\\ & \underset{\mathrm{\Delta }\mathit{\tau }\to \mathrm{0}}{lim}{\displaystyle \frac{{\mathit{\tau }}_i}{n\mathrm{!}\mathrm{\Delta }\mathit{\tau }}}\langle {\left[\mathrm{\Delta }{u}_i^{\prime }({\mathit{\tau }}_i-\mathrm{\Delta }\mathit{\tau },u^{*})-\mathrm{\Delta }{u}_i({\mathit{\tau }}_i,u^{*})\right]}^n\rangle .\end{array}\end{array}$$

In contrast to the Kramers–Moyal expansion in time domain, a minus sign on the left-hand side (lhs) of Eq. (4) has to be added for scale processes, since during evolution of the process the scale τ is decreasing. According to the Pawula theorem, the Kramers–Moyal coefficients vanish for n≥3, if D⁽⁴⁾=0 (Risken, 1996); thus the Kramers–Moyal expansion reduces to the Fokker–Planck equation (FPE):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}-{\mathit{\tau }}_i& {\displaystyle \frac{\partial }{\partial {\mathit{\tau }}_i}}p\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_j;u^{*})=\\ & -{\displaystyle \frac{\partial }{\partial \mathrm{\Delta }{u}_i}}\left[D^{\left(\mathrm{1}\right)}(\mathrm{\Delta }{u}_i,{\mathit{\tau }}_i,u^{*})p\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_j;u^{*})\right]+\\ & -{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial \mathrm{\Delta }{u}_i^{\mathrm{2}}}}\left[D^{\left(\mathrm{2}\right)}(\mathrm{\Delta }{u}_i,{\mathit{\tau }}_i,u^{*})p\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_j;u^{*})\right],\end{array}\end{array}$$

with the drift function D⁽¹⁾ accounting for the deterministic evolution of the stochastic process, whereas the so-called diffusion function D⁽²⁾ scales the amplitude of the noise term of the corresponding Langevin equation

$$\begin{array}{}\text{(7)}& \begin{array}{rl}-{\displaystyle \frac{\partial }{\partial \mathit{\tau }}}& \mathrm{\Delta }u(\mathit{\tau },u^{*})=\\ & {\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}{D}^{\left(\mathrm{1}\right)}(\mathrm{\Delta }u,\mathit{\tau },u^{*})+\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}{D}^{\left(\mathrm{2}\right)}(\mathrm{\Delta }u,\mathit{\tau },u^{*})}\cdot \mathrm{\Gamma }\left(\mathit{\tau }\right)\end{array}\end{array}$$

with the Gaussian noise Γ(τ), fulfilling 〈Γ(τ)〉=0 and also $\langle \mathrm{\Gamma }\left(\mathit{\tau }\right)\mathrm{\Gamma }\left({\mathit{\tau }}^{\prime }\right)\rangle =\mathrm{2}{\mathit{\delta }}_{\mathit{\tau }{\mathit{\tau }}^{\prime }}$. This equation directly describes the evolution of a single trajectory along the scale τ.

As can be seen, the FPE can be used to determine the factorised pdf's from Eq. (3) if the process is Markovian and higher-order Kramers–Moyal coefficients are zero. This equivalence between a cpdf with N conditions and a nested chain of several cpdf's with only two conditions, stemming from the three-point closure of a cascade process, is extremely helpful if one aims at estimating the pdf's from measurements, since the high-dimensional pdf's in Eq. (2) would require a tremendous number of realisations in order to be estimated well. For a more detailed discussion of this multipoint approach we refer to Peinke et al. (2019).

2.2 Preliminary analysis of wind speed data

Next we check if wind speed data are suitable for the reconstruction method just described. According to the right-hand side (rhs) of Eq. (3) the estimation of the double cpdf's $p\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_{i+\mathrm{1}},u^{*})$ is necessary. However, to reduce computational costs or the number of data points, it would be much more convenient to use the single cpdf's $p\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_{i+\mathrm{1}})$ by excluding the condition on the wind speed u^* to be predicted (Nawroth et al., 2010; Peinke et al., 2019). Thus the equality

$$\begin{array}{}\text{(8)}& p\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_{i+\mathrm{1}};u^{*})=p(\mathrm{\Delta }{u}_i\left|\mathrm{\Delta }{u}_{i+\mathrm{1}}\right)\end{array}$$

must hold. As can be seen from Fig. 1, Eq. (8) holds for $u^{*}\approx \mathrm{0}$ but shows a significant shift for $u^{*}\approx \mathrm{2.5}$. From this we reason that the equality in Eq. (8) cannot generally be assumed, so we have to stick to the double cpdf's $p\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_{i+\mathrm{1}},u^{*})$. Similar results have been reported for idealised turbulence (Stresing and Peinke, 2010) and sea waves (Hadjihosseini et al., 2016). On derivation of our final formula for the reconstruction of time series (see Eq. 3) an essential step was to premise the underlying scale process to be Markovian. Thus it has to be shown that for $\mathrm{\Delta }\mathit{\tau }={\mathit{\tau }}_{i+\mathrm{1}}-{\mathit{\tau }}_i\ge \mathrm{\Delta }{\mathit{\tau }}_{\mathrm{ME}}$

$$\begin{array}{}\text{(9)}& p\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_{i+\mathrm{1}};\mathrm{\dots };\mathrm{\Delta }{u}_N;u^{*})=p(\mathrm{\Delta }{u}_i|\mathrm{\Delta }{u}_{i+\mathrm{1}};u^{*})\end{array}$$

is a valid assumption. As the verification of this expression is not feasible in its generality, we limit ourselves by reducing the number of dimensions involved and just check the equality of

$$\begin{array}{}\text{(10)}& {\displaystyle }p\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_{i+\mathrm{1}};\mathrm{\Delta }{u}_{i+\mathrm{2}};u^{*})=p(\mathrm{\Delta }{u}_i|\mathrm{\Delta }{u}_{i+\mathrm{1}};u^{*}).\end{array}$$

To check this we utilise the Chapman–Kolmogorov equation (CKE) (Friedrich et al., 2011). The cpdf $p\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_{i+\mathrm{2}};u^{*})$ is estimated directly from observational data and afterwards compared with the cpdf $\stackrel{\mathrm{̃}}{p}\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_{i+\mathrm{2}};u^{*})$ obtained numerically by use of the CKE

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\stackrel{\mathrm{̃}}{p}& \left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_{i+\mathrm{2}};u^{*})=\\ & \int p\left(\mathrm{\Delta }{u}_i\right|\mathrm{\Delta }{u}_{i+\mathrm{1}};u^{*}\left)p\right(\mathrm{\Delta }{u}_{i+\mathrm{1}}|\mathrm{\Delta }{u}_{i+\mathrm{2}};u^{*}),\end{array}\end{array}$$

whereas the two cpdf's within the integral on the rhs are also directly estimated from data. Figure 2 shows that Eq. (11) only holds for Δτ≥Δτ_ME. We find a Markov–Einstein length of Δτ_ME≤0.1 s, which we are going to use henceforth.

Figure 1 Comparison of single cpdf's p(Δu₁|Δu₂) (black) and double cpdf's $p\left(\mathrm{\Delta }{u}_{\mathrm{1}}\right|\mathrm{\Delta }{u}_{\mathrm{2}};u^{*})$ (red) for $u^{*}\approx \mathrm{0}$ (a, c), $u^{*}\approx \mathrm{2.5}$ (b, d), τ₁=1 s and τ₂=2τ₁. Dashed lines indicate cuts through the cpdf's at Δu₂≈1.2. The levels of the contour plots are given in natural logarithmic scale.

Figure 2 Comparison of the double cpdf's $p\left(\mathrm{\Delta }{u}_{\mathrm{1}}\right|\mathrm{\Delta }{u}_{\mathrm{3}};u^{*})$ (estimated from data) and $\stackrel{\mathrm{̃}}{p}\left(\mathrm{\Delta }{u}_{\mathrm{1}}\right|\mathrm{\Delta }{u}_{\mathrm{3}};u^{*})$ (obtained from CKE) for τ₁=0.1 s (a) or 1 s (b) and τ₂=2τ₁, τ₃=3τ₁, $\mathrm{\Delta }\mathit{\tau }={\mathit{\tau }}_{\mathrm{2}}-{\mathit{\tau }}_{\mathrm{1}}={\mathit{\tau }}_{\mathrm{3}}-{\mathit{\tau }}_{\mathrm{2}}$ and $u^{*}=\mathrm{0}$.

2.3 Parameterisation of the Fokker–Planck equation

As was mentioned in Sect. 2.1 the FPE may be used to generate solutions for the needed cpdf's in Eq. (1). Aiming for this, one needs a parameterisation of the FPE reflecting the scale process of the real-world wind speed data. No general, physical formulation for a FPE describing the scale process of wind speeds is known; thus we use the possibility to estimate a parameterisation directly from the given data (compare Eq. 5 and Reinke et al., 2018; Peinke et al., 2019).

This way we obtain estimations of the drift and diffusion functions $D^{\left(\mathrm{1}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$, $D^{\left(\mathrm{2}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$ along the scale τ_i and for every wind speed value u^*. From these estimations one then usually finds parameterisations of the FPE by fitting appropriate polynomial surfaces to the estimated functions, which then may be used to obtain numerical solutions of the FPE.

To match the functional shape of the estimated $D^{\left(\mathrm{1}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$ and $D^{\left(\mathrm{2}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$ (see Fig. 3), we require the polynomials to be of the order of 3 and 2 respectively. We find a significant shift, denoted with ${\mathit{\gamma }}_{D^{\left(\mathrm{1}\right)}}({\mathit{\tau }}_i,u^{*})$, depending on u^* and τ_i for the drift along all scales τ_i which has to be taken into account for a parameterisation suitable to the given data.

Figure 3 Exemplary estimations of the drift and diffusion functions $D^{\left(\mathrm{1}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$ and $D^{\left(\mathrm{2}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$ for τ=65 s and $u^{*}\approx -\mathrm{1.58}$ (red), ≈0 (black) and ≈1.58 (blue). Fits for D⁽¹⁾ and D⁽²⁾ were achieved by applying Eqs. (12) and (13).

We set for $D^{\left(\mathrm{1}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{D}^{\left(\mathrm{1}\right)}& (\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})=\\ & d_{\mathrm{10}}({\mathit{\tau }}_i,u^{*})+d_{\mathrm{11}}({\mathit{\tau }}_i,u^{*})\left[u-{\mathit{\gamma }}_{D^{\left(\mathrm{1}\right)}}({\mathit{\tau }}_i,u^{*})\right]+\\ & d_{\mathrm{13}}({\mathit{\tau }}_i,u^{*}){\left[u-{\mathit{\gamma }}_{D^{\left(\mathrm{1}\right)}}({\mathit{\tau }}_i,u^{*})\right]}^{\mathrm{3}}\end{array}\end{array}$$

and for $D^{\left(\mathrm{2}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{D}^{\left(\mathrm{2}\right)}& (\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})=\\ & d_{\mathrm{20}}({\mathit{\tau }}_i,u^{*})+d_{\mathrm{21}}({\mathit{\tau }}_i,u^{*})u+d_{\mathrm{22}}({\mathit{\tau }}_i,u^{*})u^{\mathrm{2}}.\end{array}\end{array}$$

Similar findings were obtained for other systems (Hadjihosseini et al., 2016; Stresing and Peinke, 2010), where the dependence on u^* was limited only to the drift function $D^{\left(\mathrm{1}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$. For wind speed data the contribution of u^* is clearly more complex.

Furthermore we check the validity of the Pawula theorem, requiring D⁽⁴⁾=0. As one can see in Fig. 4, the fourth Kramers–Moyal coefficient is slightly larger than zero, but negligible compared to the magnitude of the diffusion function D⁽²⁾.

Figure 4 Exemplary estimations of the second and fourth Kramers–Moyal coefficients $D^{\left(\mathrm{2}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$ and $D^{\left(\mathrm{4}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$ for τ=65 s.

Next, we check if the parameterisation of the FPE, given by Eqs. (12) and (13), performs well in describing the underlying scale process, before one uses it for the reconstruction scheme presented in Sect. 2.1. This can easily be done by comparing cpdf's estimated directly from the data and the ones obtained from numerical solutions of the FPE. The latter one is not carried out by common finite-difference schemes but by an iterative approach (for first works see Renner et al., 2001; Wächter et al., 2003). For a (very) small step size in scale Δτ the functions $D^{\left(\mathrm{1}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$ and $D^{\left(\mathrm{2}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$ can be assumed to be constant in τ, leading to an exact solution (see Risken, 1996) for the cpdf

$$\begin{array}{}\text{(14)}& \begin{array}{rl}p& (\mathrm{\Delta }{u}_j,{\mathit{\tau }}_i-\mathrm{\Delta }\mathit{\tau }|\mathrm{\Delta }{u}_i,{\mathit{\tau }}_i,u^{*})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\sqrt{\mathit{\pi }{D}^{\left(\mathrm{2}\right)}({\mathit{\tau }}_i,u^{*})\mathrm{\Delta }\mathit{\tau }}}}\cdot \\ & \exp \left[-{\displaystyle \frac{{\left(\mathrm{\Delta }{u}_j-\mathrm{\Delta }{u}_i-D^{\left(\mathrm{1}\right)}({\mathit{\tau }}_i,u^{*})\mathrm{\Delta }\mathit{\tau }\right)}^{\mathrm{2}}}{\mathrm{4}{D}^{\left(\mathrm{2}\right)}({\mathit{\tau }}_i,u^{*})\mathrm{\Delta }\mathit{\tau }}}\right],\end{array}\end{array}$$

describing the transition between wind speed increments of a larger scale τ_i and a smaller scale τ_i−Δτ. By iteratively combining this so-called short-time propagator (STP) with the CKE (see Eq. 11), one is able to obtain cpdf's $p(\mathrm{\Delta }{u}_j,{\mathit{\tau }}_j|\mathrm{\Delta }{u}_i,{\mathit{\tau }}_i,u^{*})$ for arbitrary large-scale differences ${\mathit{\tau }}_i-{\mathit{\tau }}_j\gg \mathrm{\Delta }\mathit{\tau }$. In theory this procedure would lead to exact solutions of the FPE, but since one is limited to finite step sizes Δτ, this method of course only provides numerical approximations of the true cpdf's.

Comparing cpdf's estimated directly from the data and from the numerical solution, we see (Fig. 5) that our proposed estimation in terms of $D^{\left(\mathrm{1}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$ and $D^{\left(\mathrm{2}\right)}(\mathrm{\Delta }u,{\mathit{\tau }}_i,u^{*})$ is well suited to describe the underlying scale process. Thus we use this parameterisation for the reconstruction scheme.

Figure 5 Isoline plots of a cpdf $p(\mathrm{\Delta }{u}_i,{\mathit{\tau }}_i|\mathrm{\Delta }{u}_j,{\mathit{\tau }}_j;u^{*})$ estimated from the data (black) and from the numerical solution of the FPE (red) for τ_i=1.6 s and τ_j=3.2 s and for $u^{*}\approx \mathrm{0}$ and $u^{*}\approx \mathrm{2}$.

2.4 Results of the multiscale reconstruction

Alternatively to the presented approach to obtain the cpdf's from numerical solutions of the FPE, it is of course possible and much less cumbersome to estimate them directly from observational data. (Note that due to the use of the FPE, the obtained pdf's are less noisy and extend to large values as seen in Fig. 5.) In this section we will present the results for the multipoint reconstruction achieved for u^*, yielded from both approaches.

Figure 6 Comparison of original wind speed time series u^* (black, a) with the reconstructed one using Eq. (8) (red, b). The vertical line marks the transition from the N=128 starting values to the simulated wind speeds.

To start the reconstruction scheme, we provided a short piece of the original time series of N wind speeds to the algorithm to compute the pdf $p\left(u^{*}\right|u_{\mathrm{1}},\mathrm{\dots },u_N)$, from which the first point of our simulation can be drawn. The reconstruction scheme is then shifted by one time step τ and applied again. By iteratively applying our method a new artificial time series of arbitrary length can be generated. After N iterations of the reconstruction scheme no data from measurements are required any more to generate new values of u^*. From visual comparison in Fig. 6 one finds a realistic looking simulated time series of u^*, retaining the characteristics of dynamics on the smaller scales as well as the ones of the larger scales.

Figure 7 Comparison of the marginal increment pdf's computed from the empirical data (black), the simulated data using the directly estimated pdf's (red) and the simulated data using the cpdf's obtained from numerical solution of the FPE (blue). The scales range from Δτ_ME to N⋅τ_ME; they are explicitly: 2ⁱ⋅Δτ_ME with $(i=\mathrm{0},\mathrm{1},\mathrm{\dots },\mathrm{7})$ and Δτ_ME=0.1 s. For better visualisation the pdf's were shifted along the vertical axis.

To confirm this visual impression in a quantitative way, we compare the increment pdf's p(Δu_i,τ_i) obtained from the reconstructed and the measured time series. The synthetic time series were generated by using the cpdf's estimated directly both from data and from numerical solutions of the FPE.

As shown in Fig. 7 the increment pdf's yielded from both stochastic simulations nicely coincide with the pdf's from observations. This attests that the presented reconstruction scheme is able to capture the complex dynamics of wind speeds, characterised by a gradual shift of increment pdf's of a Gaussian-like shape (larger scales) to increment pdf's of a heavy-tailed shape (smaller scales).

A striking difference between the increment pdf's from the stochastic simulation can be noted. Whereas the tails of the original pdf's are systematically underestimated by the reconstruction using the directly estimated pdf's, the reconstruction from the numerically obtained cpdf's is able to keep track of the tails of the original pdf. This stems from the fact, mentioned above, that by estimating a pdf from observational data one in general underestimates the outer tails, since there are only a few measured points available. Considering the estimation of cpdf's, the estimation error of course worsens, which is even more severe when additional conditioning is applied, like in our case in terms of the additional condition on u^*.

The tails of cpdf's $p(\mathrm{\Delta }{u}_i,{\mathit{\tau }}_i|\mathrm{\Delta }{u}_j,{\mathit{\tau }}_j;u^{*})$ computed from the family of FPEs can be extrapolated to areas where no measurement points are available. This effect is a direct consequence of the CKE (Eq. 11), as the tails are the product of quite well-estimated less probable but not rare events.

Certainly this approach indirectly suffers from the limited number of observations as well, as the estimation of the functions $D^{\left(\mathrm{1}\right)}(u,r,u^{*})$ and $D^{\left(\mathrm{2}\right)}(u,r,u^{*})$ is based on observational data, too (see Eq. 5). Furthermore the parameterisation of these functions is always only a approximation of the real drift and diffusion functions, introducing deviations from the real-world system.

But we see from Fig. 7 that the majority of increments are well grasped; only the occurrence of a few rare events are underpredicted. A more detailed investigation of the pdf's shows the pdf's obtained by the FPE deviate from the original shape for the largest scale $\mathit{\tau }=N\cdot \mathrm{\Delta }{\mathit{\tau }}_{\mathrm{ME}}$ but performs better in the outer wings of the pdf.

From Fig. 8 a better understanding of the reconstruction method can be gained. The pdf used for the simulation $p\left(u^{*}\right|u_{\mathrm{1}},\mathrm{\dots },u_N)$ is not stationary, even though it is computed from completely stationary cpdf's (see Eq. 8) and changes sensitively with respect to the N past values. While the snapshot pdf (shown as black circles) at the time marked by the black vertical line has a rather clear shape, it undergoes a change, becoming broader (red line and red circles). Due to the spreading of the wings of the pdf, values of u^* of larger magnitude become more likely to be drawn. This may lead to a distinct shift of the wind speed values to $u^{*}\le \mathrm{0}$ as seen in this example. After this transition the broadness decreases (blue line and circles) again and a tendency to relax back to $u^{*}=\mathrm{0}$ can be seen. Furthermore the increased broadness of the red pdf (red line) can be seen as an early warning signal that the wind speed is prone to fluctuate in a stronger way.

Figure 8 Evolution of $p\left(u^{*}\right|u_{\mathrm{1}},\mathrm{\dots },u_N)$ during reconstruction. Horizontal lines indicate snapshots of the pdf used for drawing the next sample of u^*. The colours of the horizontal lines correspond to the snapshot pdf's.

3 Extension to non-stationary wind speeds

In the preceding part for the reconstruction scheme, block-wise normalised wind speeds with a window length of 1 min were used. These blocks were defined by common mean wind speed $\stackrel{\mathrm{‾}}{U}$ and standard deviation σ_U. For the normalised wind speed u, we showed how to generate new time series; see Fig. 9.

There are different methods to generate more general non-stationary wind data. Knowing the slow variation in $\stackrel{\mathrm{‾}}{U}\left(t\right)$ and σ_U(t), the drift and diffusion coefficients D⁽ⁱ⁾ are taken as a slowly changing function of $D^{\left(i\right)}(\mathrm{\Delta }u,\mathit{\tau },u^{*},\stackrel{\mathrm{‾}}{U},{\mathit{\sigma }}_U)$. If due to the normalisation of U to u the coefficients D⁽ⁱ⁾ are in a good approximation independent of $\stackrel{\mathrm{‾}}{U}$ and σ_U, the slow variation in the real wind conditions can be incorporated over the back-transformation of the newly estimated values $U^{*}=({\mathit{\sigma }}_U\cdot u^{*})+\stackrel{\mathrm{‾}}{U}$. The slow dynamics of $\stackrel{\mathrm{‾}}{U}\left(t\right)$ and σ_U(t) may be given from measured data, meteorological simulations or other modelling.

Figure 9 (a, b) Comparison between measured non-stationary speeds U and reconstructed non-stationary wind speeds. (c, d) Increment and one point pdf's of original and reconstructed U.

A third possibility is a self-adaptive procedure which we show here. Instead of using given values of $\stackrel{\mathrm{‾}}{U}\left(t\right)$ and σ_U(t), the intrinsic fluctuation of these quantities are used. Given an initial pair $\left(\stackrel{\mathrm{‾}}{U}\right(t),{\mathit{\sigma }}_U(t\left)\right)$ estimated over 1 min from measured data, a time series of non-stationary wind speeds U^*(t) is obtained by applying the above-mentioned back-transformation to the generated wind speeds u^* from our algorithm. For the upcoming simulation window of 1 min length, a new pair $\left(\stackrel{\mathrm{‾}}{U}\right(t^{\prime }),{\mathit{\sigma }}_U(t^{\prime }\left)\right)$ is estimated from the just generated block of wind speed data U^*(t). This procedure is carried out until a time series of non-stationary wind speeds of desired length is obtained. In Fig. 9 such a time series is shown together with the statistical analysis of the increment pdf's and the marginal pdf of the non-stationary wind speed U. We observe that the pdf's of the reconstructed time series match the shape of the empirical pdf's for both the increments and the wind speeds. At this point we would like to emphasise that we do not aim to create copies of historical wind speeds, but rather to be able to generate stochastically equivalent time series.

4 Conclusions

We presented a stochastic approach based on multipoint statistics to generate surrogate short-time wind speed fluctuations with stochastic processes. Note these stochastic processes can be estimated self-consistently from given data. By using the normalised wind speeds u^* and wind speed increments Δu(τ_i), Δu(τ_j) from two separate scales τ_i and τ_j, a three-point closure to the complex systems of wind speeds was achieved.

It was shown that our method works well in describing the dynamics of block-wise normalised wind speeds u^* along scales τ_i in terms of a stochastic scale process, governed by a family of Fokker–Planck equations. This separation of the fluctuations from the mean values is similar to the Reynolds-averaged Navier–Stokes (RANS) approach widely used in fluid dynamics (Frisch, 2004), with the difference that we have a description of the underlying stochastic process of the fluctuation and “only” lack the mean values. With the modified reconstruction (see Sect. 3) we are able to generate mean values on the basis of past values of the reconstructed time series, yielding realistic non-stationary wind speed time series U^*. As the typical response times of wind turbines and their control systems have duration of seconds to minutes, our reconstructed wind data are suitable for investigation of many dynamical effects of the wind energy conversion process. Note for these times the wind energy system may be driven into non-stationary response dynamics. Thus for these timescales it is important to have to maximal information on all statistical details of the driving wind source. The stochastic model presented here is based on multipoint statistics and is able to capture small-scale intermittency effects, extreme events and clustering of fluctuations, up to now not addressed in wind energy research.

For timescales larger than the response times of wind turbines, the turbines operate with fully adapted control systems in a stationary state. To estimate effects, like loads, of such stationary states, the temporal order of the states becomes unimportant. It is sufficient to know how often which wind situation emerges. Thus the knowledge of the valid Weibull distribution $p\left(\stackrel{\mathrm{‾}}{U}\right)$ should be sufficient. Note our result here indicates that it would be better to extend the Weibull distribution to the joint probability $p(\stackrel{\mathrm{‾}}{U},{\mathit{\sigma }}_U)$.

Finally we emphasise that the presented stochastic multipoint approach to small-scale wind speed fluctuations should encompass automatically extreme short-term wind fluctuations, commonly added to wind investigations in terms of standard or multi-year gusts. These methods can be applied easily to other wind quantities like the temporal behaviour of shears, or wind veers, eventually combined in higher-dimensional stochastic processes (Siefert and Peinke, 2006). The results reported in Ali et al. (2019) show that such a stochastic modelling can also be used for wake flows.

Appendix A: Derivation of simpler cpdf's

To link Eq. (2) to the idea of an underlying turbulent cascade, we identify the conditional pdf on the left-hand side (lhs) with the cpdf $p(u^{*},t^{*}|\mathrm{\Delta }{u}_{\mathrm{1}},t^{*}-{\mathit{\tau }}_{\mathrm{1}};\mathrm{\dots };\mathrm{\Delta }{u}_N,t^{*}-{\mathit{\tau }}_N)$. Thus the numerator of the right-hand side (rhs of Eq. 2) can be rewritten as

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}p& (u^{*},t^{*};u_{\mathrm{1}},t^{*}-{\mathit{\tau }}_{\mathrm{1}};\mathrm{\dots };u_N,t^{*}-{\mathit{\tau }}_N)=\\ & p(u^{*},u^{*}-u_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{1}};\mathrm{\dots };u^{*}-u_N,{\mathit{\tau }}_N)\end{array}\end{array}$$

and the nominator as

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}p& (u_{\mathrm{1}},t^{*}-{\mathit{\tau }}_{\mathrm{1}};\mathrm{\dots };u_N,t^{*}-{\mathit{\tau }}_N)=\\ & p(u_{\mathrm{1}};u_{\mathrm{1}}-u_{\mathrm{2}},{\mathit{\tau }}_{\mathrm{2}}-{\mathit{\tau }}_{\mathrm{1}};\mathrm{\dots };u_{\mathrm{1}}-u_N,{\mathit{\tau }}_N-{\mathit{\tau }}_{\mathrm{1}}).\end{array}\end{array}$$

The identity of the expressions in Eqs. (A1) and (A2) can mathematically rigorously be shown, as done in Nawroth et al. (2010), but intuitively speaking the sequences on the lhs and rhs must yield the same joint pdf, since the increments on the rhs have a common reference point u^* or u₁. Next we factorise the joint pdf's from Eqs. (A1) and (A2) by iteratively using cpdf's

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}p& (u^{*};\mathrm{\Delta }{u}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{1}};\mathrm{\dots };\mathrm{\Delta }{u}_N,{\mathit{\tau }}_N)=\\ & p(\mathrm{\Delta }{u}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{1}}|\mathrm{\Delta }{u}_{\mathrm{2}},{\mathit{\tau }}_{\mathrm{2}};\mathrm{\dots };\mathrm{\Delta }{u}_N,{\mathit{\tau }}_N;u^{*})\cdot \\ & p(\mathrm{\Delta }{u}_{\mathrm{2}},{\mathit{\tau }}_{\mathrm{2}}|\mathrm{\Delta }{u}_{\mathrm{3}},{\mathit{\tau }}_{\mathrm{3}};\mathrm{\dots };\mathrm{\Delta }{u}_N,{\mathit{\tau }}_N;u^{*})\mathrm{\cdots }\\ & p(\mathrm{\Delta }{u}_{N-\mathrm{1}},{\mathit{\tau }}_{N-\mathrm{1}}|\mathrm{\Delta }{u}_N,{\mathit{\tau }}_N;u^{*})\cdot \\ & p(\mathrm{\Delta }{u}_N,{\mathit{\tau }}_N|u^{*})\cdot p(u^{*})\end{array}\end{array}$$

and with $\stackrel{\mathrm{̃}}{\mathrm{\Delta }{u}_i}:=u(t^{*}-{\mathit{\tau }}_{\mathrm{1}})-u(t^{*}-{\mathit{\tau }}_i)$ with the timescale τ_i−τ₁:

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}& p(u_{\mathrm{1}};\stackrel{\mathrm{̃}}{\mathrm{\Delta }{u}_{\mathrm{2}}},{\mathit{\tau }}_{\mathrm{2}}-{\mathit{\tau }}_{\mathrm{1}};\mathrm{\dots };\stackrel{\mathrm{̃}}{\mathrm{\Delta }{u}_N},{\mathit{\tau }}_N-{\mathit{\tau }}_{\mathrm{1}})=\\ & p(\stackrel{\mathrm{̃}}{\mathrm{\Delta }{u}_{\mathrm{2}}},{\mathit{\tau }}_{\mathrm{2}}-{\mathit{\tau }}_{\mathrm{1}}|\stackrel{\mathrm{̃}}{\mathrm{\Delta }{u}_{\mathrm{3}}},{\mathit{\tau }}_{\mathrm{3}}-{\mathit{\tau }}_{\mathrm{1}};\mathrm{\dots };\stackrel{\mathrm{̃}}{\mathrm{\Delta }{u}_N},{\mathit{\tau }}_N-{\mathit{\tau }}_{\mathrm{1}};u_{\mathrm{1}})\cdot \\ & p(\stackrel{\mathrm{̃}}{\mathrm{\Delta }{u}_{\mathrm{3}}},{\mathit{\tau }}_{\mathrm{3}}-{\mathit{\tau }}_{\mathrm{1}}|\stackrel{\mathrm{̃}}{\mathrm{\Delta }{u}_{\mathrm{4}}},{\mathit{\tau }}_{\mathrm{4}}-{\mathit{\tau }}_{\mathrm{1}};\mathrm{\dots };\stackrel{\mathrm{̃}}{\mathrm{\Delta }{u}_N},{\mathit{\tau }}_N-{\mathit{\tau }}_{\mathrm{1}};u_{\mathrm{1}})\\ & \mathrm{\cdots }p(\stackrel{\mathrm{̃}}{\mathrm{\Delta }{u}_{N_{\mathrm{1}}}},{\mathit{\tau }}_{N-\mathrm{1}}-{\mathit{\tau }}_{\mathrm{1}}|\stackrel{\mathrm{̃}}{\mathrm{\Delta }{u}_N},{\mathit{\tau }}_N-{\mathit{\tau }}_{\mathrm{1}};u_{\mathrm{1}})\cdot \\ & p(\stackrel{\mathrm{̃}}{\mathrm{\Delta }{u}_N},{\mathit{\tau }}_N-{\mathit{\tau }}_{\mathrm{1}}|u_{\mathrm{1}})\cdot p(u_{\mathrm{1}}).\end{array}\end{array}$$

A further step in reducing the dimensionality of the involved pdf's can be performed upon assuming the scale process to be Markovian; i.e. there exists a timescale separation $\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{ME}}={\mathit{\tau }}_j-{\mathit{\tau }}_i$ (j>i), where

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}& p(\mathrm{\Delta }{u}_i,{\mathit{\tau }}_i|\mathrm{\Delta }{u}_j,{\mathit{\tau }}_i+\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{ME}};\mathrm{\dots };\mathrm{\Delta }{u}_n,{\mathit{\tau }}_i+n\cdot \mathrm{\Delta }{\mathit{\tau }}_{\mathrm{ME}};u^{*})\\ & =p(\mathrm{\Delta }{u}_i,{\mathit{\tau }}_i|\mathrm{\Delta }{u}_j,{\mathit{\tau }}_i+\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{ME}};u^{*})\end{array}\end{array}$$

holds. The timescale separation Δτ_ME is often called the Markov–Einstein length (Einstein, 1905), and for various systems its existence could be shown empirically, ranging from jet streams in laboratory experiments (Renner et al., 2001), (Reinke et al., 2018) to large geophysical systems such as ocean gravity waves (Hadjihosseini et al., 2016).

Appendix B: Parameterisation of D⁽¹⁾ and D⁽²⁾

Here we present the polynomial coefficients used to parameterise the first and second Kramers–Moyal coefficients $D^{\left(\mathrm{1}\right)}(\mathrm{\Delta }{u}_i,{\mathit{\tau }}_i,u^{*})$ and $D^{\left(\mathrm{2}\right)}(\mathrm{\Delta }{u}_i,{\mathit{\tau }}_i,u^{*})$ (see Eqs. 12 and 13).

Coefficients for the drift function $D^{\left(\mathrm{1}\right)}(\mathrm{\Delta }{u}_i,{\mathit{\tau }}_i,u^{*})$:

$$\begin{array}{}\text{(B1)}& {\displaystyle }\begin{array}{rl}{d}_{\mathrm{10}}& =c_{\mathrm{0},d_{\mathrm{10}}}\cdot {\mathit{\tau }}_i+c_{\mathrm{1},d_{\mathrm{10}}}\cdot u^{*}\cdot {\mathit{\tau }}_i^{{\stackrel{\mathrm{̃}}{c}}_{\mathrm{1},d_{\mathrm{10}}}}+c_{\mathrm{2},d_{\mathrm{10}}}\cdot {\mathit{\tau }}_i^{{\stackrel{\mathrm{̃}}{c}}_{\mathrm{2},d_{\mathrm{10}}}}\\ & \cdot u^{*\mathrm{2}}+c_{\mathrm{3},d_{\mathrm{10}}}\cdot u^{*\mathrm{3}},\end{array}\text{(B2)}& {\displaystyle }{d}_{\mathrm{11}}=c_{\mathrm{0},d_{\mathrm{11}}}\cdot {\mathit{\tau }}_i+c_{\mathrm{1},d_{\mathrm{11}}}\cdot {\mathit{\tau }}_i^{{\stackrel{\mathrm{̃}}{c}}_{\mathrm{1},d_{\mathrm{11}}}}+c_{\mathrm{2},d_{\mathrm{11}}}\cdot {\mathit{\tau }}_i^{{\stackrel{\mathrm{̃}}{c}}_{\mathrm{2},d_{\mathrm{11}}}}\cdot u^{*\mathrm{2}},\text{(B3)}& {\displaystyle }{d}_{\mathrm{13}}=c_{\mathrm{0},d_{\mathrm{13}}}\cdot {\mathit{\tau }}_i^{{\stackrel{\mathrm{̃}}{c}}_{\mathrm{0},d_{\mathrm{13}}}}+c_{\mathrm{1},d_{\mathrm{13}}}\cdot u^{*},\text{(B4)}& {\displaystyle }{\mathit{\gamma }}_{D^{\left(\mathrm{1}\right)}}=c_{\mathrm{1},{\mathit{\gamma }}_{D^{\left(\mathrm{1}\right)}}}\cdot {\mathit{\tau }}_i^{{\stackrel{\mathrm{̃}}{c}}_{\mathrm{1},{\mathit{\gamma }}_{D^{\left(\mathrm{1}\right)}}}}\cdot u^{*},\end{array}$$

with $c_{\mathrm{0},d_{\mathrm{10}}}=-\mathrm{0.006}$, $c_{\mathrm{1},d_{\mathrm{10}}}=-\mathrm{0.888}$, ${\stackrel{\mathrm{̃}}{c}}_{\mathrm{1},d_{\mathrm{10}}}=\mathrm{0.098}$, $c_{\mathrm{2},d_{\mathrm{10}}}=\mathrm{0.137}$, ${\stackrel{\mathrm{̃}}{c}}_{\mathrm{2},d_{\mathrm{10}}}=\mathrm{0.019}$, $c_{\mathrm{3},d_{\mathrm{10}}}=-\mathrm{10.566}$, $c_{\mathrm{0},d_{\mathrm{11}}}=-\mathrm{1.656}$, $c_{\mathrm{1},d_{\mathrm{11}}}=-\mathrm{0.018}$, ${\stackrel{\mathrm{̃}}{c}}_{\mathrm{1},d_{\mathrm{11}}}=-\mathrm{8.853}e-\mathrm{05}$, $c_{\mathrm{2},d_{\mathrm{11}}}=-\mathrm{0.268}$, ${\stackrel{\mathrm{̃}}{c}}_{\mathrm{2},d_{\mathrm{11}}}=\mathrm{1.671}$, $c_{\mathrm{0},d_{\mathrm{13}}}=-\mathrm{0.005}$, ${\stackrel{\mathrm{̃}}{c}}_{\mathrm{0},d_{\mathrm{13}}}=\mathrm{0.012}$, $c_{\mathrm{1},d_{\mathrm{13}}}=\mathrm{1.023}$, $c_{\mathrm{1},{\mathit{\gamma }}_{D^{\left(\mathrm{1}\right)}}}=\mathrm{0.341}$, ${\stackrel{\mathrm{̃}}{c}}_{\mathrm{1},{\mathit{\gamma }}_{D^{\left(\mathrm{1}\right)}}}=\mathrm{0.247}$.

And for the diffusion function $D^{\left(\mathrm{2}\right)}(\mathrm{\Delta }{u}_i,{\mathit{\tau }}_i,u^{*})$

$$\begin{array}{}\text{(B5)}& {\displaystyle }\begin{array}{rl}{d}_{\mathrm{20}}& =c_{\mathrm{0},d_{\mathrm{20}}}\cdot {\mathit{\tau }}_i^{{\stackrel{\mathrm{̃}}{c}}_{\mathrm{0},d_{\mathrm{20}}}}+c_{\mathrm{1},d_{\mathrm{20}}}\cdot {\mathit{\tau }}_i^{{\stackrel{\mathrm{̃}}{c}}_{\mathrm{1},d_{\mathrm{20}}}}\cdot u^{*}+c_{\mathrm{2},d_{\mathrm{20}}}\\ & \cdot {\mathit{\tau }}_i^{{\stackrel{\mathrm{̃}}{c}}_{\mathrm{2},d_{\mathrm{20}}}}\cdot u^{*\mathrm{2}},\end{array}\text{(B6)}& {\displaystyle }{d}_{\mathrm{21}}=c_{\mathrm{0},d_{\mathrm{21}}}\cdot {\mathit{\tau }}_i^{{\stackrel{\mathrm{̃}}{c}}_{\mathrm{0},d_{\mathrm{21}}}}+c_{\mathrm{1},d_{\mathrm{21}}}\cdot {\mathit{\tau }}_i^{{\stackrel{\mathrm{̃}}{c}}_{\mathrm{1},d_{\mathrm{21}}}}\cdot u^{*},\text{(B7)}& {\displaystyle }{d}_{\mathrm{22}}=c_{\mathrm{0},d_{\mathrm{22}}}\cdot {\mathit{\tau }}_i^{{\stackrel{\mathrm{̃}}{c}}_{\mathrm{0},d_{\mathrm{22}}}}+c_{\mathrm{1},d_{\mathrm{22}}}\cdot {\mathit{\tau }}_i^{{\stackrel{\mathrm{̃}}{c}}_{\mathrm{1},d_{\mathrm{22}}}}\cdot u^{*},\end{array}$$

with $c_{\mathrm{0},d_{\mathrm{20}}}=\mathrm{0.024}$, ${\stackrel{\mathrm{̃}}{c}}_{\mathrm{0},d_{\mathrm{20}}}=-\mathrm{0.0001}$, $c_{\mathrm{1},d_{\mathrm{20}}}=\mathrm{0.0002}$, ${\stackrel{\mathrm{̃}}{c}}_{\mathrm{1},d_{\mathrm{20}}}=\mathrm{1.076}$, $c_{\mathrm{2},d_{\mathrm{20}}}=\mathrm{1.573}$, ${\stackrel{\mathrm{̃}}{c}}_{\mathrm{2},d_{\mathrm{20}}}=\mathrm{1.622}$, $c_{\mathrm{0},d_{\mathrm{21}}}=\mathrm{0.002}$, ${\stackrel{\mathrm{̃}}{c}}_{\mathrm{0},d_{\mathrm{21}}}=-\mathrm{0.001}$, $c_{\mathrm{1},d_{\mathrm{21}}}=\mathrm{1.104}$, ${\stackrel{\mathrm{̃}}{c}}_{\mathrm{1},d_{\mathrm{21}}}=\mathrm{1.395}$, $c_{\mathrm{0},d_{\mathrm{22}}}=\mathrm{0.042}$, ${\stackrel{\mathrm{̃}}{c}}_{\mathrm{0},d_{\mathrm{22}}}=\mathrm{0.002}$, $c_{\mathrm{1},d_{\mathrm{22}}}=\mathrm{0.555}$, ${\stackrel{\mathrm{̃}}{c}}_{\mathrm{1},d_{\mathrm{22}}}=\mathrm{0.364}$.