Measurements of atmospheric turbulence have been studied and found to deviate from a Gaussian process, in particular regarding the velocity increments over small time steps, where the tails of the pdf are exponential rather than Gaussian. Principles forextreme event counting and the occurrence of cascading events are presented. Empirical extreme statistics agree with Rice’s exceedence theory, when it is assumed that the velocity and its time derivative are independent. Prediction based on the assumptionthat the velocity is a Gaussian process underpredicts the rate of occurrence of extreme events by many orders of magnitude, mainly because the measured pdf is non-Gaussian. Methods for simulation of turbulent signals have been developed and theircomputational efficiency are considered. The methods are applicable for multiple processes with individual spectra and probability distributions. Non-Gaussian processes are simulated by the correlation-distortion method. Non-stationary processes areobtained by Bezier interpolation between a set of stationary simulations with identical random seeds. Simulation of systems with some signals available is enabled by conditional statistics. A versatile method for simulation of extreme events has beendeveloped. This will generate gusts, velocity jumps, extreme velocity shears, and sudden changes of wind direction. Gusts may be prescribed with a specified ensemble average shape, and it is possible to detect the critical gust shape for a givenconstruction. The problem is formulated as the variational problem of finding the most probable adjustment of a standard simulation of a stationary Gaussian process subject to relevant event conditions, which are formulated as linear combination of pointsin the realization. The method is generalized for multiple correlated series, multiple simultaneous conditions, and 3D fields of all velocity components. Generalization are presented for a single non-Gaussian process subject to relatively simpleconditions, i.e. gusts and velocity jumps. Further generalizations for simulation of multiple correlated non-Gaussian processes are suggested.