Bernoulli-Langevin Wind Speed Model for Simulation of Storm Events

The motivation for this work is derived from the goal of supporting air traffic controllers by automatic assistance systems that, in addition to supporting traffic control under normal conditions, are able to also derive suggestions from real-time data during extreme weather (rare or “Xevents” [1], [2]) for optimal anticipative air traffic control and management actions at airports. A good prediction of, e.g. wind/gust speed profiles during extreme low-pressure (storm) events by real-time updating of weather forecast should help minimise performance disruption. It would be realised, e.g. by proposing appropriate arrival/departure rescheduling [3], [4] and re-direction to alternative runways or airports as modification of the original plan. Kantz et al. [5], [6], provide an example for a Markov chain based stochastic time series analysis of turbulent surface wind and gust for prediction of disruptive disturbance as precondition for anticipatory wind turbine rotor blade control on the time frame of seconds. The main modeling focus here is on the average behavior with time scale of less than an hour up to the order of a day or so, i.e. the duration of a typical autumn storm event in Middle Europe that in addition to (averaged) wind speed increase includes the increase of dangerous gust as speed fluctuations. Recent advanced wind speed models that include extended time frames of a couple of hours and more are based, e.g. on sophisticated combinations of artificial intelligence (AI) methods like neural networks (i.e. model-free estimators), classical autoregressive moving average methods, Kalman filters, and other approaches for optimal control (see [7], [8] and references therein). Wind speed models using stochastic differential equations similar to our approach, however based on the (stationary) Ornstein-Uhlenbeck (O-U) process (ds​/​dt= s˙=−a s+ε(t) [9]) with linear drift term and Gaussian noise ε(τ) are derived in [10]. By means of fitting empirical storm data with a heuristic approximation to our nonlinear Bernoulli model, we provide evidence that the considered class of extreme storm events in fact require a non-stationary dynamical description allowing for state (or phase) transition between two (stationary) states (i.e. normal wind speed and storm with increased speed fluctuations). This situation is not covered by the O-U dynamics that could, however, complement the storm model before and after the Xevent.

In this paper a simple nonlinear dynamics stochastic wind speed model ds/dt=B(k(t); s, ε(t)) (s=modulus of velocity vector) for extreme storm events is derived and compared with historical data. It is based on the second-degree Bernoulli equation [11] (s˙+p(t)s=q(t)s2), however with a single time-dependent rate parameter k=–p=–q and a Langevin random force term ε(t). The model formalises state transition times (pre-storm to storm ≈1/|k|) below 1 h including fluctuations, with typical data sample times of 30 min. A related (second-degree) approach was used as a simplified model of Laser action that describes the external pump energy induced non-equilibrium (disorder-to-order) phase transition at the laser threshold from normal light (random photon emission) to coherent photon emission [12]. The equation in this case describes the time variation of photon number in the laser cavity n˙= κn−κ1n2  with control parameters κ, κ₁ where increasing pump power induces a sign change κ(t)<0→>0 at the onset of coherent laser light emission. In terms of the Landau (equilibrium) theory of phase transitions the full semiclassical laser theory with third-degree nonlinearity corresponds to a second-order (continuous) transition [12]. In what follows we will discuss our analytically treatable Bernoulli-Langevin (B-L) model with a single rate parameter k(t) controlling the wind speed variation and compare model-based simulations with logistic function fits to empirical data.

After providing the dynamical Langevin model, we will first derive the heuristic analytic approximation to the general solution of the former for fitting averaged wind speed profiles of historical storm events. The results provide evidence for our underlying two-state transition hypothesis formalised by the Bernoulli equation. Empirical stochastic speed fluctuation parameters are derived from the fit residuals of two example wind profiles measured for historic low-pressure events. These data also yield probability densities from the empirical histograms of acceleration fluctuations which are used as random external force input into Monte Carlo simulations with the B-L equation for reproducing the empirical non-stationary stochastic time series.

The Langevin equation of our dynamical wind speed model is given by the combination of the quadratic Bernoulli force function as deterministic drift term B(v,t)) with a single time-dependent rate parameter k(t) as external control and the (normalised) randomly fluctuating external force ε(t):

with normalised speed v=s/s_a, acceleration s˙​/​sa=v˙=dv​/​dt , and stationary solutions v₁*=0, v₂*=1=upper asymptotic speed v_a (s=s_a)=1. B(v, t) is assumed to fulfill the potential condition, i.e. the potential function U(v)=k(t) v² (v/3–1/2) is obtained via integration of the (conservative force) relationship B(v)=–dU/dv. With 0≤v≤1 the system state moves within rigid potential boundaries at fixed points v_1,2*. The stability of v* at the potential minima is controlled by the sign of k(t), e.g. v₁*/v₂* unstable/stable for k>0. In terms of equilibrium thermodynamics (with U equivalent to free energy) the state change v₁*↔v₂* with sign change of k(t) corresponds to a first-order phase transition due to the discontinuous change of entropy S=–dU/dk at k=0: S₁(v₁*)=0, S₂(v₂*)=1/6. Following [12], ε(t) is defined by random pushes modelled by δ-correlated Gaussian noise, represented by the autocorrelation function K_tt′=<ε(t) ε(t′)>=Qδ(t–t′). Via the Wiener-Khintchine theorem the power spectral density is related to K_tt′ through Q=J_ω, with variance σ²=J_ω/t_S, t_S=sampling time. Equation (1) allows for simulating empirical time series with time-dependent rate (control) parameter k(t) shaping the wind speed profile with state transition time constant τ(t) (=1/|k(t)|).

In order to provide evidence for our nonlinear dynamics two-state approach for modeling the transition between normal and storm state, we will first analyse two empirical wind profiles of historical storm events. The speed profiles are fitted by a heuristic approximation to the general solution of the Bernoulli drift model [see below (6) and (7)] using constant rate parameter k>0 that turns (1) into the logistic (or Verhulst) equation for population growth under limited resources or carrying capacity (here replaced by asymptotic maximum speed s_a>s_max during the storm). Correspondingly, we may approximate the wind speed increase dv/dt as proportional to the momentary value v(t) (with initial value v₀) yielding exponential time dependence of speed. For the upper asymptotic limit v_a=1 (≈v_max for long storm duration T>>τ) we have 0<v₀≤v(t)≤v_max<v_a. The reduction of the v(t) increase rate when approaching v_a is realised by the factor (1–v(t)) yielding a sigmoid function as solution

Different rate parameters k_i>0, k_d<0 are required for fitting the increasing and the decreasing speed section of the storm time series, with the two time shift values μ_i, μ_d for maximum and minimum of averaged speed increase and decrease, respectively. μ=half value time of the state transition: v(t=μ)=1/2. According to the sign of k the (inverse) sign of the exponent represents increasing (k=+1/τ>0) or decreasing (k=–1/τ>0) state transition (sigmoid) with time constant τ=1/|k|. It is easily verified that, e.g. for the negative exponent (k>0) the stationary states are lim v(t→–∞)=0=v₁* and lim v(t→+∞)=1=v_a=v₂*. Moreover, v(t=0)=v₀=1/(1+exp(μ/τ) (≈v₁* for large μ), yielding as definition for the transition time μ of the considered time series:

The analytical approximation to the general solution of (1) [see (6) and (7)] to be used for time series fitting is constructed by considering the two discrete instants of time only (μ_i, μ_d). Written in dimensional magnitudes via s=vs_a, the wind speed profile is approximated by the sum of the two corresponding sigmoids, to be used as heuristic model for time series fitting and as approximation for the solution of (1):

By inspecting the asymptotic values, one confirms that s(lim t→–∞)=s_d and s(lim t→+∞)=s_i. Transition with time constant τ_i (k_i>0) from the (stationary) pre-storm state s_d to the (quasi) stationary storm state s_max<s_a is followed, after a delay Δμ=μ_d–μ_i=storm duration by wind speed decrease with τ_d (k_d<0) from v_max back to the (stationary) post-storm wind speed s_i≈s_d. An asymptotic (stationary) speed s_max≈s_a is achieved only if Δμ>>τ_i, τ_d (see the “Xaver” event below). Otherwise, the speed starts decreasing before s_a is approached (see the “Christian” event below).

Equation (4) allows to attribute characteristic times, time constants, stationary wind speed levels, and speed maximum to the Xevent. An estimate of the (maximum) quasi stationary speed value s_c in the center between increase and decrease at time t_c=(μ_i+μ_d)/2=μ_i+Δμ/2 is obtained as

For large storm duration Δμ it approaches s_max=s_i+s_d (=2s_a if s_i≈s_d). This sigmoid approximation to the solution of (1) [see (6) and (7)] provides a fit function for speed/gust measurements with estimates of two triplets of constant parameters {s_i, μ_i, τ_i=1/|k_i|}, {s_d, μ_d, τ_d=1/|k_d|}, with (quasi) stationary pre-/post-storm speed values s_d/s_i. We will see below that within the B-L model with continuous control k_c=k(t=t_c)=0 defines a critical point of the hypothesised phase transition (Xevent (storm)→post storm wind).

Figures 1 and 2 depict two wind speed and gust profiles of extreme storm (X-)events “Christian” and “Xaver” measured at Hamburg airport in October and December 2013, respectively. Wind time series were taken from the aeronautical standard METAR (Meteorological Terminal Air Report) weather reports of the German Weather Service [http://www.skybrary.aero/index.php/wind_velocity_reporting, http://www.dwd.de/DE/fachnutzer/luftfahrt/luftfahrt_node.html], issued in 30 min intervals for aviation purposes, measured according to ICAO Annex 3 (International Civil Aviation Organisation) at standard height of 10 m (in knots, 1 kt=1 nm/h=1.852 km/h=0.514 m/s). Speed values are the last 10 min averages of the intervals, and gust values are the maximum speed values (exceeding 2 min averages) within the 10 min intervals. The first Xevent occurred around noon to afternoon of October 28, 2013 (Δμ≈5 h) with maximum wind/gust speed (>30/50 kt) around t_c=(μ_inc+μ_dec)/2≈1:00 pm. The second one extended over 2 days from December 5 to 6, 2013 (Δμ≈40 h) with maximum wind/gust speed around 23/45 kt. Both figures include nonlinear least squares fits to wind and gust profiles using model (4) that exhibit significant parameter estimates within 95% confidence intervals.

Figure 1:

Empirical data and model fit. Wind speed data (bottom, 10 min averages, every 30 min) and gust (top) for low-pressure storm event “Christian” on October 28, 2013. Logistic model fits with parameter estimates for wind (inserted) and gust speed profiles according to (4) (smooth blue lines), with 95% confidence intervals (dotted lines). inc=i=increasing, dec=d=decreasing. For details, see text and Table 1.

Figure 2:

Empirical data and model fit. Wind speed data (bottom, 10 min averages, every 30 min) and gust (top) for low-pressure storm event “Xaver” on December 5–6, 2013. Logistic model fits with parameter estimates for wind (inserted) and gust speed profile according to (4) (smooth blue lines), with 95% confidence intervals (dotted lines). inc=i=increasing, dec=d=decreasing. For details, see text and Table 1.

Table 1 collects the fit parameters for wind and gust speed profiles of the two examples. As expected, the maximum wind speeds of approximately 23–24 kt (“Xaver”) and 31–32 kt (“Christian”) as derived from the empirical data correspond precisely to the theoretical model based estimates at the time near speed maximum t_c≈(μ_dec+μ_inc)/2. s_max is calculated with (5) from the wind/gust parameter estimates of Table 1 based on the time constants τ, delays μ, and stationary speed before and after the Xevent: s_max(Xaver)=23.3 kt, s_max(Christian)=32.5 kt. For Δμ>>τ (as in the “Xaver” case) the maximum speed is approximated reasonably well by the sum of the sigmoid asymptotes s_max=s_i+s_d≈2s_a. An estimate of the empirical wind speed fluctuations is obtained through the fit residuals. Gaussian density fits N(μ, σ) of speed histogram residuals yield significant results (parameter estimates±SD) at χ² test error probability level α=5%, i.e. χ²<χ²_α or quantile α_S(χ²_sample)>0.05 both for “Christian” (μ, σ, α_S)=(−0.1±0.3, 1.8±0.2, 0.53), and “Xaver” (μ, σ, α_S)=(0±0.5, 3.5±0.3, 0.32). The statistics of the random speed fluctuations will be used below for generating the stochastic force input ε(t) into simulations with the dynamical model (1).

While σ(residuals) averages over the complete dataset, one interesting observation is the difference (s_max(gust)–s_max(wind average))>(s_i,d(gust)–s_i,d(wind average)). This is reflected by the apparently instationary speed fluctuations around the fitted wind profiles which seem to increase with speed, i.e. fluctuation increase during transition to the upper speed state with sign change of the control parameter at k(t=t_c). Increase of fluctuations is typical for non-equilibrium disorder-to-order phase transitions like in the laser with incoherent to coherent light emission, and it may be an indication for a corresponding dynamical background for the considered class of extreme storm events. The promising results of data fitting with the sigmoid approximation (4) to the general solution of (1) is taken as even stronger evidence for a dynamic non-stationary two-state approach and non-equilibrium phase transition. The sigmoid fit function is, however, only of limited value for real-time prediction of wind/gust speed because it relies on a priori knowledge of the two triples of logistic function parameters, derived from previously measured data. For anticipative management actions the stochastic nonlinear dynamics (1) has to be solved allowing for real-time predictive filtering and extrapolation of the observed time series.

The deterministic (drift) part of (1) allows for modelling the transition between the two (quasi) stationary states (normal wind speed and storm) via continuous sign change k(t)>0↔k(t)<0. The well-known general solution [11] of our second-degree Bernoulli equation is given by

with integration constant c determined by the initial speed v₀. An explicit analytical model may be derived by integrating (6), if a not too complex time dependence of the rate parameter k(t) can be determined. k(t) should provide a change of sign from positive to negative, corresponding to the transition from increasing to decreasing speed during the time interval≈Δμ of high speed state≈v_max. As a simple ansatz we try a sinusoidal rate parameter k(t)=k_max sin{2π(t–t₁)/T} between start and end times of the storm profile (t₁, t₂; T=t₂–t₁), with sign change of k(t) (change of potential minima from v₁* to v₂*) at critical times t₁ (normal wind becomes instable: start of storm transition), and t_c=t₁+T/2 (storm state becomes instable). Constant speed v₁=v₀ and k=0 is assumed for simplicity before (t≤t₁) and after the storm (v₂=v₀ for t>t₂). Introducing this function into (6) yields for the normalised wind speed v(t) during the storm event

It fulfills the initial condition v₁=v₀ at t₁=t₀ if c=1/v₀–1=s_a/s₀–1 with s(t)=vs_a in dimensional units (kt). This yields also for t₂=t₁+T: v₂=v(t₂)=v₀. Within the potential picture [see text following (1)] the low speed boundary of U(k,v) at v₁*=0 is shifted to v₀. Maximum (normalised) speed v_max<v_a for this symmetric model is obtained as

For T>>τ_min=1/k_max the asymptotic speed converges to v_max≈v_a=1, as expected. Figure 3 compares the solution (7) of (1) for “Christian” (T=10 h) and “Xaver (T=42 h) after multiplication with a scaling factor s_max/v_max=s_a for absolute speed dimension (kt). s_max is the estimated maximum of the empirical speed value according to (5), and v_max is the maximum normalised speed according to (8) at time t=t_c=(t₁+t₂)/2.

Figure 3:

Deterministic wind speed profiles as obtained with (7), dimensionalised in (kt) through scaling factor s(t)=v(t) s_max/v_max, for comparison with empirical data and logistic fit of “Christian” (k-period T=10 h) and “Xaver” (k-period T=42 h) events (Figs. 1 and 2).

The storm duration Δμ of the heuristic approximation (4) and (5) is replaced in the Bernoulli model (7) and (8) by the single period T of the sinusoidal rate function. The B-L parameter set {v₀, k_max(s_max), T} is matched to reproduce the empirical data of “Christian” (Fig. 1) and “Xaver” (Fig. 2) with regard to the deterministic part of the 10 min averaged speed profile. While the parameters of the sigmoid approximation (4) are estimated from historical time series fitting, the initial B-L parameter estimates of (7) including σ(ε) may be derived in practical applications from weather forecast and real-time speed measurement which will also be the basis for parameter updates and predictive filtering like in [8].

The simulation of the fluctuating speed, in addition to the averaged speed, may be achieved by numerically solving (1) with inclusion of the random external force ε(t) as additive noise, i.e. interpreting ε(t) as random speed increments. A noise distribution for numerical computer experiments is obtained from the empirical time series (Figs. 1 and 2) via the speed increments Δs=s(t+Δt)–s(t) corresponding to accelerations or random forces Δs/Δt acting on the system state s in intervals Δt=30 min and resulting in random deviations from the deterministic two-state speed profile (7). The corresponding ε(t) histograms, together with Gaussian N(μ, σ) fits for the two events, are depicted in Figure 4a and b.

Figure 4:

Histograms with Gaussian fit N(μ, σ) of the 10 min average-speed increments (kt/h) measured in 30 min intervals (from Figs. 1 and 2) for (a) “Christian” and (b) “Xaver”. Fit significant for p(χ²)>0.05. Std. dev. σ quantifying external random force for numerical integration of Langevin (1). For details, see text.

The large deviation at zero mean in Figure 4b (Xaver) with χ²-test rejection of N(μ, σ) hypothesis is due to missing data points (see Fig. 2). As mentioned above the residuals of the fits in Figures 1 and 2 provide some indication of non-stationarity of the noise, with speed fluctuation amplitudes increasing during transition to the storm state. A theoretical estimate of the expected speed dependence of the noise distribution may be obtained from the stationary solution of the Fokker-Planck equation [12] with drift (B(v)) and diffusion term (Gaussian δ-correlated noise ε(t) of strength Q) corresponding to the B-L dynamics (1). For our one-dimensional case the stationary solution yields for the speed distribution the probability density [12]

with potential U=–dB/dv [see text following (1)], natural boundary conditions (f vanishes for v→±∞), normalisation constant N, and arbitrary constant U_c=0. After transition through v=v_max the speed decreases corresponding to dv/dt=B(v)≈–(v/τ)(1–v) with maximum acceleration at half value time μ if τ=constant (logistic approximation). With the exponent in (9) as ratio of a quadratic potential and variance, modified by a cubic term (U(v)=(v²/2τ)(1–2/3v)), f(v) may be interpreted as Gaussian density with variance ~Q(v). This results in an increase of the effective variance that is illustrated for three speed values (v<<1, v=1/2, v=1): exp{–v²/τQ}, exp{–v²/1.5τQ}, exp{−v²/3τQ}. The increase of the standard deviation is obtained as σ(v) = σ0/1−2v/3 (with σ₀=√(τQ/2)) corresponding to speed fluctuations growing by a factor 1.7 with v, which compares reasonably well with the empirical results of Figures 1 and 2. Figure 5 depicts example results for “Christian” and “Xaver” of numerical integration of the B-L equation (1) (using Matlab ode45) with sinusoidal k(t) rate parameter ansatz and δ-correlated Gaussian noise ε(t) which exhibit fluctuating time series for the two simulated quite different wind speed profiles, also in reasonable agreement with the empirical results.

Figure 5:

Example of stochastic wind speed profiles of “Christian” and “Xaver” events (Figs. 1 and 2) corresponding to analytic solution in Figure 3, as obtained from Monte Carlo simulation with the stochastic (1). Stochastic Langevin force derived from empirical speed increment distributions (Fig. 4, see text for details).

Major features of the two events are reproduced quite well although for “Xaver” the speed decrease starts too early, probably due to the early sign change of the sinusoidal k(t) in combination with large fluctuations. Detailed results of course vary between simulation runs with different seed values of the random number generator. For achieving even better agreement, one option is the consideration of asymmetry ({s_i, τ_i}≠{s_d, τ_d}), as apparent in the collection of fit parameters of Table 1. This could be achieved, e.g. through a weak time dependence of k_max and different pre-/post-storm speed levels. Like the Laser equation an improved B-L model should allow for different control functions k(t), k₁(t) of the linear and quadratic term of the Bernoulli equation which also would open the possibility of smooth transitions between the B-L (k₁=1) and the O-U model (k=constant, k₁=0 for pre-/post-storm phase) and for relating the control parameters k, k₁ to their physical/meteorological background.

In conclusion, a nonlinear dynamics non-stationary wind speed model for extreme storm situations (Xevents) is derived. It is based on a two-state phase transition hypothesis formalised by the deterministic second-degree Bernoulli dynamics with a single (time-dependent) control parameter k(t) and stochastic δ-correlated Langevin force. Initial evidence for a non-stationary nonequilibrium phase transition underlying the considered Xevents is obtained by fitting of empirical speed increase-decrease time series with a sum of two logistic functions (sigmoids, k=constant) as heuristic approximation to the general solution of the dynamical model. Speed and gust data were taken from standard aeronautical METAR weather reports [http://www.dwd.de/DE/fachnutzer/luftfahrt/luftfahrt_node.html] at Hamburg airport of two storm events in 2013. Fit estimates are obtained for characteristic wind-profile parameters: increase and decrease rate of speed change (|k|=inverse time constants τ), pre- and post-storm stationary speed values and corresponding event start/stop times with storm duration Δμ. Maximum speed within the sigmoid fit model is determined by the asymptotic, i.e. (quasi) stationary values before (and after) the storm and duration Δμ, i.e. the different time shifts μ of the two sigmoids. An analytic solution of the Bernoulli equation formalising the dynamical two-state hypothesis is provided with a sinusoidal time dependence of the rate parameter k(t) (with period T≈Δμ defining the storm duration) as most simple ansatz for the transition between the two stationary states. Despite simplifications (symmetric time constants and equal pre-/post-storm speed), the deterministic part of the wind profile is reproduced reasonably well for the two quite different event histories. Numerical simulations with inclusion of the stochastic Langevin force based on empirical fluctuations of measured speed increments reproduced the essential features and differences of the two events including noise. As a phase transition approach for extreme storm (X-)events with time scales of a few hours up to days the B-L dynamics may be seen as complimentary to the linear and stationary O-U process [10] and other approaches aiming at short-term prediction and (automatic) control application [6], [7], [8]. Our results appear supportive for a nonequilibrium phase transition during certain extreme storm events, in analogy to the laser transition from random to coherent photon emission [12]. Of course, this hypothesis requires the present initial evidence to be reproduced with more empirical data and a more detailed theoretical treatment, e.g. following [12]. With further evidence our results are expected to contribute to more resilient anticipative air traffic management through predictive filtering of real-time weather data (e.g. using extended Kalman filter) also under extreme storm situations.