Home / Archiv / Realistic simulation of ocean surface using wave spectra
Document Actions

No section

Realistic simulation of ocean surface using wave spectra

  1. Jocelyn Fréchot LaBRI - Laboratoire bordelais de recherche en informatique


We present a method to simulate ocean surfaces away from the coast, with correct statistical wave height and direction distributions. By using classical oceanographic parametric wave spectra, our results fit real world measurements, without depending on them. Since wave spectra are independent of the ocean model, Gerstner parametric equations and Fourier transform method can be used with them. Moreover, since they are simple to use and need very few parameters, they allow easy production of ocean surface animations usable in movies and games. We explain how to accurately sample them, to achieve oceanic waves in deep water according to given wind parameters, in a realistic way.

  1. published: 2007-11-14


1.  Introduction

The reproductions of many natural phenomena is still an open problem in Computer Graphics. In particular, the realistic simulation of oceanic surfaces continues to be of great interest. Many fields of activity rely on it: virtual reality, movies and games, but also sailing [ CGG01 ] and oceanographic simulations, among others.

While beautiful pictures and animations such as wave refraction [ GlS00 ] [ GM02 ] are produced, there is still a need for oceanic scenes with realistic wave features. Models based on fluid mechanics equations [ EMF02 ] [ MMS04 ] should give the more realistic results, but are inadequate for large oceanic surfaces.

For twenty years, procedural models for wave representation in Computer Graphics have been developed and improved [ Igl04 ]. The Gerstner parametric equations and their Fourier transform rewriting are well known by the oceanographic community. Based on the linear wave theory, they allow representation of natural wave shapes and moves in deep water above the capillary size. As they not describe any net mass transport, they are limited to non-breaking waves and to scenes without violent storms.

Wave height measurement data allow to simulate a given sea state. In the more general case, oceanic surface statistics are taken in account by the use of parametric wave spectra. For given wind parameters, they indicate the distribution of wave energy as a function of frequency. However, none of the existing methods are able to correctly reproduce sea states described by these spectra.

In this paper, we explain how wave spectra are related to oceanic surface energy and we detail a method to deduce wave heights from it. We propose a solution for optimal spectrum sampling, suitable for parametric and Fourier transform models. This method enables the accurate simulation of oceanic surface. Users can easily produce realistic wave animations, with respect to supplied weather conditions.

The structure of this paper is as follows. Section 2 summarizes ocean models for Computer Graphics. In section 3 , we detail our method and give classical parametric spectrum examples. We handle wave statistics in section 4 . We discuss our results in section 5 , andconclude in section 6 .

2.  Ocean models

2.1.  Definitions and relationships

The surface elevation, η, is the vertical displacement from the water level of the oceanic surface at rest at given point and time. The wave amplitude, a, is the maximum surface elevation due to this wave. The wave height, H, is the vertical distance between a trough and an adjacent crest. For a single wave, H = 2a. The wavelength, λ, is the distance between two successive crests. The period, T, is the time interval between the passage of successive crests through a fixed point. The frequency, f = 1/T, is the number of crests that pass a fixed point in 1 second. The wave speed or phase speed, c = λ/T, is the speed of wave crests or troughs. It can also be expressed as c = ω/κ, where ω = 2πf is the angular velocity or angular frequency, and κ = 2π/λ is the wavenumber. The water depth, d, is the vertical distance between the floor and the water level at rest.

The concept of deep water is wave dependent and is related to the ratio of the water depth to the wavelength of the wave. This comes from the transcendental expression


where g ≈ 9.807 m⋅s -2 is the standard acceleration of gravity. Since tanh(x) ≈ 1 when x > π, and tanh(x) ≈ x when x < π/10, water is considered deep if d/λ > 1/2 and shallow if d/λ < 1/20. This leads to some simplifications in the expressions of the terms given above (see table 1 ).

Table 1. Relations between oceanographic terms. The subscript ∞ indicates a deep water term.

2.2.  Parametric equations

The base model for ocean waves simulation comes from the linear (or small-amplitude) wave theory by [ Air45 ], widely used in oceanic engineering, and in Computer Graphics by [ Pea86 ]. It describes waves with a sinusoidal shape, which corresponds to calm weather conditions. Considering a location x0 lying on a 1D surface at rest, the elevation z = η at time t due to a wave with wavenumber κ, amplitude a and angular velocity ω = is


When the steepness of the wave increases, its crests become sharper and its troughs flatter. Therefore, [ FR86 ] used a more realistic description based on trochoids, made by [ vG04 ] and by [ Ran63 ]. The surface equation is now


where z0 is the surface elevation at rest (typically 0). Note that the important point here is the phase term difference of -π/2 between x(x0,t) and z(x0,t), not the use of a particular sine or cosine function. This Lagrangian model describes the trajectory in a vertical plane of a particle (x,z) around its position at rest (x0,z0). Summing several waves and extending equation 3 to a 2D surface gives


where = (x,y) is the horizontal particle position at time t, 0 = (x0,y0) its position at rest, = (κxy) a wave vector, i.e. a vector with magnitude || || = κ and direction of the considered wave propagation and = /|||| is the unit vector of . Because of the presence of uniformly distributed random phase terms, () ∈ [0,2π[, this model is also known as random waves. Any set of 0 can be taken, which allows the use of adaptive mesh, for optimal surface sampling [ HNC02 ].

2.3.  Inverse Fourier transform

Another method for computing equation 4 is the 2D inverse discrete Fourier transform, used by Mastin, Watterberg and Mareda [ MWM87 ], Premože and Ashikhmin [ PA01 ], and Tessendorf [ Tes01 ]. The set of N x M complex numbers Fn,m is transformed into a set of complex numbers fp,q by


where p ∈ {1,..,N} and q ∈{1,..,M} . This is of particular interest because the fast Fourier transform (FFT) algorithm is an efficient way to compute these sum. However, in comparison with the previous method, its usage is more restrictive. The set of 0 is defined as a regular grid of N x M points with dimensions Lx x Ly and origin at the center, such that 0 = (nx0 Lx/N,my0 Ly/M), where nx0 ∈ [-N/2,n/2[ and mx0 ∈ [-M/2,M/2[ are integers. We rewrite equation 4 as


where () and () correspond, respectively, to the imaginary part and real part of the expression, and


Rather than the random phase term , authors commonly use a normally distributed random amplitude (see section 4 ). To comply with equation 5 , the set of wave vectors must be


where nκ and mκ are defined as nx0 and mx0 above. The resulting surface is periodic in all directions and can be used as a tile [ Tes01 ].

To decrease the FFT computation time by a factor two, real numbers can be used instead of complex ones. This is achieved by adding two waves with same amplitude and travelling in opposite directions [ Tes01 ]. Equation 7 becomes


This results in a single stationary wave with time dependent amplitude. This phenomenon is known as standing wave. The drawback of this method is that the generated surface is only half the size of the expected one. The second half is obtained with a rotation of the first one by 180 degrees around the origin, leading to pattern duplications.

3.  Wave spectra

To be able to use the models discussed so far, we need a function a(). For the parametric model, we also need to know which set of wave vectors characterizes the surface we want to generate. The more natural method is to use records of surface elevation measurement data used by oceanic engineering [ TG02 ]. Measurements are made withaccelerometers mounted on buoys, satellite altimeters, high frequency radars and others techniques. A spectral analysis of the collected data is done using the Fourier transform from time domain to frequency domain. This decomposition of the surface elevation gives a histogram of wave frequencies ω, with corresponding amplitudes and possibly directions, known as frequency amplitude spectrum. The same analysis method can be used with oceanic surface photographs [ Tes01 ] [ TG02 ].

An amplitude spectrum can be used to generate a surface with the same characteristics as the related one. However, to get several oceanic surfaces with different sea states, as many spectra are needed. Instead, we rely on parametric wave spectra, which can fit any environment conditions, like wind speed and fetch.

3.1.  Wave energy and amplitudes

Oceanographers are mostly concerned by wave energy rather than amplitude. For this purpose, they use another type of spectrum, S(ω), called wave frequency energy spectrum or wave spectrum, which gives information about wave energy distribution as a continuous function of frequency. This is a smoothed version of a modified amplitude spectrum (see below). Directional information are taken into account with a directional wave spectrum:


where θ is the direction of the wave and D(ω,θ) is a directional spreading function, defined such that


for all values of ω. The Pierson-Moskowitz parametric wave spectrum (see section 3.2 ) has been regularly used in the field of Computer Graphics since [ MWM87 ] as a direct evaluation of wave amplitudes. But this method cannot reflect statistical characteristics of the oceanic surface. First, this is a continuous function, so it needs to be integrated to be used with a discrete model. Second, wave spectrum does not give direct information about wave amplitude but energy.

The mean energy per unit surface area is E = ρg var(η), where ρ is the water density and var(η) the variance of the surface elevation. The wave spectrum is defined to be related to the surface energy by


For a zero mean process, like the sinusoidal wave modelfrom equation 2 , the variance is


where xrms is the root mean square of the values. For a general process, this becomes x2 rms- 2 . The mean amplitude of the trochoidal wave model from equation 3 is = -(κa2/2, which can be reasonably neglected in most cases for practical purposes. So, we apply equation 13 to this model as well for simplicity. The rms of a sinusoidal or trochoidal wave with amplitude a is


As we see the oceanic surface as the sum of many waves with independent random phases, the total variance of the surface elevation is the sum of the variance of each wave. So, taking a finite range of frequencies minmax] such that it is representative of the whole energy, i.e.


we divide it in N samples of width Δωn and tag the mean frequency of each of them as ωn to get the relation


Regarding equations 12 and 15 , this is rewritten as


and so we found the amplitude of each wave with


To take into account the direction of the wave, we sample the whole directional spectrum E(ω,θ):


So, for each wave with angular frequency ω and direction θ, we can now compute the corresponding amplitude a(ω,θ).

3.2.  Parametric spectra

Parametric spectrum models are empirical expressions with user defined parameters that have been found to fit ocean surface elevation measurements. Most common parameters include wind speed and direction and fetch. They are useful to simulate realistic ocean waves without measurement data, but can also be altered to fit specific data.

The base of modern parametric wave spectrum was found by [ Phi57 ] and has the form


where α = 0.0081 is called the Phillips constant. Taking into account wind speed, [ PM64 ] found an expression that describes a theoretical fully developed sea with infinite fetch (see figure 1 ). The wind has blown with no change over a large area for several hours, and waves growth is almost null. The spectrum is


where ωp is the frequency of the peak of the spectrum, i.e. dSPM/dω = 0 for ωp . Usual value for ωp is

where U10 is the wind speed at a height of 10 meters above the sea surface.

Figure 1. Pierson-Moskowitz frequency spectrum SPM(ω) with different wind speed values U10 .

Pierson-Moskowitz frequency spectrum SPM(ω) with different wind speed values U10.

For a fetch limited sea, where waves continue to grow, the JONSWAP spectrum has been found by [ HBB73 ] to be more appropriate (figure 2 ). This is the previous spectrum with an additional term:



Here, γ is a peak enhancement factor, i.e. when same values are used for α and ωp in both spectra, we get γ = SJp/SPMp). The parameter σ controls the width of the peak. Usual parameter values for this spectrum are

where F is the fetch in meters,

γ = 3.3, but can vary between 1 and 7, and

As the fetch is an important element of the sea state description, we use the JONSWAP spectrum, rather than the Pierson-Moskowitz one.

Figure 2. JONSWAP frequency spectrum, SJ(ω), with a wind speed of U10 = 15m⋅s-1 , and different fetch values F.

JONSWAP frequency spectrum, SJ(ω), with a wind speed of U10 = 15m⋅s-1, and different fetch values F.

Other parametric spectrum models can be found in oceanographic literature, including double peaked spectrum for swell coming from distant storm. These spectra are often expressed as a function of frequency f = ω/(2π). To preserve integral equality, conversion between spectra is done using substitution:


For the dispersion relation, we use the expression found by [ LH62 ]:


where θw is the main direction of the spectrum, usually the direction of the wind, and N(s(ω)) is a normalization factor defined as

and Γ() is the gamma function (figure 3 ). The function s(ω) controls the sharpness of the directional spreading and has been defined by [ MTS75 ] as


Figure 3. 3D plot of the JONSWAP spectrum, EJ(ω,θ), with U10 = 25m⋅s-1 , F = 100 km and θw = 0 rad.

3D plot of the JONSWAP spectrum, EJ(ω,θ), with U10 = 25m⋅s-1, F = 100 km and θw = 0 rad.

For a simple use, only wind speed, direction and fetch need to be given to the spectrum. In order to make the spectrum curve fit particular data, like oceanic measurements, the default parameter values can be changed, until the desired curve shape is reached (figure 4 ).

Figure 4. JONSWAP frequency spectrum, SJ(ω), showing how different parameter values alter the curve shape (U10=15m⋅s-1,F=50km).

JONSWAP frequency spectrum, SJ(ω), showing how different parameter values alter the curve shape (U10=15m⋅s-1,F=50km).

3.3.  Spectrum sampling

We now look for a selection of wave frequencies and directions. For the FFT based model, the size and sample number of the rendered surface determine the set of wave vectors to be used (see equation 8 ). For the parametric one, there is not such restriction. We sample a spectrum E() in the Cartesian wavenumber domain (figure 5). To preserve integral equality, spectrum conversion is done according to substitution rule. First, we convert the frequency spectrum Sω(ω) to a wavenumber spectrum Sκ(κ):


Figure 5. The same spectrum as in figure 3 , but in the wave vector domain κx x κy (i.e. EJ( The same spectrum as in figure 3 , but in the wave vector domain κx x κy (i.e. EJ()). )).

The same spectrum as in figure 3

Second, we convert the wavenumber directional spectrum Eκ,θ(κ,θ) = Sκ(κ)D(ω = θ) from polar coordinates to Cartesian ones:


Finally, we get


we select waves taht are the most representative of the spectrum energy dispersion. For this purpose, we sample the wave spectrum over a limited domain xmin, κxmax] x [κymin, κymax]. We want this domain to be the most representative of the whole spectrum, i.e. we want it to be as dense as possible. So, by iteration, we found one of the smallest domains that contains roughly 95% of the total energy.

For simplicity, a regular sampling could be used by taking a constant size

Figure 6. The same spectrum as in figure 3 , but in the wave vector domain κx x κy (i.e. EJ( The same spectrum as in figure 3 , but in the wave vector domain κx x κy (i.e. EJ()). )).

The same spectrum as in figure 3

4.  Statistical properties of the ocean surface

Measurements of oceanic wave statistics have shown that, to a reasonable accuracy, the surface elevation η follows a normal distribution. Attributes of oceanic waves we simulate are consistent with these observations. If we look at the oceanic surface as an infinite sum of waves with infinitely small amplitude, we see that the wave spectrum is a probability density function of wave frequencies and directions. Furthermore, the variance of the surface elevation is finite (equation 12 ). So, by virtue of the central limit theorem, we know that the more waves we sum, the closer to the normal distribution is the simulated surface elevation.

Another observed characteristic of oceanic waves is that their heights follow a Rayleigh distribution, as noted by [ FR86 ]. This distribution is used with a parameter σ = Hs/2, where Hs ≈ 4 is known as the significant height. From this, it can be shown that the probability that a wave has a larger height than Hs is exp(-2) ≈ 0.1353.

Since they use defined amplitudes, the ocean models presented in section 2 are known as deterministic methods. [ TCC84 ] showed that uniformly distributed random phase terms should be replaced with random amplitudes. An expression like a cos(κx0 - ωt + ) becomes


where r1 and r2 are random numbers from the standard normal distribution, i.e. with a mean of zero and a standard deviation of one, R = follows a Rayleigh distribution and

However, since the values of x0 and t depend on arbitrary origins, we can use the uniformly distributed term in place of .

This non-deterministic method has better statistical properties than the previous one. In order to use it with parametric and FFT methods, we rewrite equations 4 and 6 according to equation 28 .

5.  Results and discussion

We made a simple real-time implementation of ocean models. We focused on wave shapes and animations, without considering effects like Fresnel reflectivity and transmissivity, or foam and spray (figure 7 ). Of course, rendering could have been achieved with classical high quality techniques as well. Interactions of objects with ocean surfaces are not specially handled by the Lagrangian model we use, and are beyond the scope of this paper.

Simulation of the interactions of objects with the ocean surface is limited, due to the Langrangian nature of the model. Velocity vectors can easily be computed at each point of the surface grid, and can be used to simulate small object motion. More complex interactions will require additional information and resort to fluid dynamics.

Since our method preserves the main spectrum energy, the global structure of the rendered surface is independent of the number of waves we sum. For both ocean models, the sea state we get is always consistent with the provided wind parameters, which cannot be achieved with other existing methods. Since these parameters are the only ones needed to have full control of the method, there is no need for the user to adjust the resulting surface by trial and error.

As previously noted, the FFT model can be particularly tedious to use. To catch a reasonable part of the spectrum, the grid length and the number of samples have to be carefully chosen, whatever are the desired surface characteristics. For example, taking 512 samples up to a frequency of 25 rad⋅s-1 requires a grid length of no more than 25 m. Furthermore, as this leads to a lowest frequency of about 1.5 rad⋅s-1, the spectrum peak may be under-sample.

We have tested our implementation with a 3 GHz Pentium 4 PC and a Radeon 9200 graphic board. With the FFT model, we got about 65 fps and 13 fps with, respectively, a grid of 128 x 128 and 256 x 256 samples. With the parametric equations, we got 8 fps with a regular grid of 128 x 128 samples and 50 waves. And taking less than 100 waves leads to poor detailed results. Clearly, only the FFT can reach interactive rate. Although we do not implement an adaptive surface mesh, it seems obvious this could not compete with the FFT speed. Parametric equations should be kept for non-interactive rendering, since they are easiest to use.

Figure 7.  FFT with a grid of 512 x 512 samples. Top: grid length is 50 m and wind speed is 5 m⋅s-1. Middle: grid length is 450 m and wind speed is 15 m⋅s-1. Bottom: grid length is 1500 m and wind speed is 25 m⋅s-1. White square is 1m wide.

FFT with a grid of 512 x 512 samples. Top: grid length is 50 m and wind speed is 5 m⋅s-1. Middle: grid length is 450 m and wind speed is 15 m⋅s-1. Bottom: grid length is 1500 m and wind speed is 25 m⋅s-1. White square is 1m wide. FFT with a grid of 512 x 512 samples. Top: grid length is 50 m and wind speed is 5 m⋅s-1. Middle: grid length is 450 m and wind speed is 15 m⋅s-1. Bottom: grid length is 1500 m and wind speed is 25 m⋅s-1. White square is 1m wide.

FFT with a grid of 512 x 512 samples. Top: grid length is 50 m and wind speed is 5 m⋅s-1. Middle: grid length is 450 m and wind speed is 15 m⋅s-1. Bottom: grid length is 1500 m and wind speed is 25 m⋅s-1. White square is 1m wide. FFT with a grid of 512 x 512 samples. Top: grid length is 50 m and wind speed is 5 m⋅s-1. Middle: grid length is 450 m and wind speed is 15 m⋅s-1. Bottom: grid length is 1500 m and wind speed is 25 m⋅s-1. White square is 1m wide.

FFT with a grid of 512 x 512 samples. Top: grid length is 50 m and wind speed is 5 m⋅s-1. Middle: grid length is 450 m and wind speed is 15 m⋅s-1. Bottom: grid length is 1500 m and wind speed is 25 m⋅s-1. White square is 1m wide. FFT with a grid of 512 x 512 samples. Top: grid length is 50 m and wind speed is 5 m⋅s-1. Middle: grid length is 450 m and wind speed is 15 m⋅s-1. Bottom: grid length is 1500 m and wind speed is 25 m⋅s-1. White square is 1m wide.

6.  Conclusion

We have presented a method for accurate wave energy spectrum sampling that allows realistic ocean surface synthesis and animation. For given wind parameters, the wave heights and directions are computed such that statistical properties of the resulting surface are correct. Since it does not rely on any ocean model, this method is suitable for Gerstner equations and Fourier transforms.

7.  Acknowledgments

The author wish to thank Bertrand le Saëc and Jean-Christophe Gonzato for rereading this paper.


[Air45] George Biddell Airy Tides and waves Encyclopaedia Metropolitana,  1845Vol. 5,  no. 192,  241—396 B. Fellowes London.

[CGG01] Jean-Marc Cieutat Jean-Christophe Gonzato, and Pascal Guitton A new efficient wave model for maritime training simulator Proceedings of the 17th Spring conference on Computer graphics,  p. 2022001isbn 0-7695-1215-1.

[EMF02] Douglas Enright Stephen Marschner, and Ronald Fedkiw Animation and rendering of complex water surfaces ACM Transactions on Graphics,  21 (2002)no. 3736—744 ACM Press issn 0730-0301.

[FR86] Alain Fournier and William T. Reeves A simple model of ocean waves SIGGRAPH Computer Graphics,  20 (1986)no. 475—84issn 0097-8930.

[GlS00] Jean-Christophe Gonzato and Bertrand le Saëc On modeling and rendering ocean scenes Journal of visualisation and computer simulation,  11 (2000)no. 127—37Januaryissn 1049-8907.

[GM02] Manuel N. Gamito and F. Kenton Musgrave An accurate model of wave refraction over shallow water Computers & Graphics,  26 (2002)no. 2291—307Aprilissn 0097-8493.

[HBB73] Klaus Hasselmann T. P. Barnett E. Bouws H. Carlson D. E. Cartwright K. Enke J. A. Ewing H. Gienapp D. E. Hasselmann P. Kruseman A. Meerburg P. Müller D. J. Olbers K. Richter W. Sell, and H. Walden Measurements of wind-wave growth and swell decay during the Joint North Sea Project (JONSWAP) Ergänzungsheft zur Deutschen Hydrographischen Zeitschrift,  1973no. 12A8.

[HNC02] Damien Hinsinger Fabrice Neyret, and Marie-Paule Cani Interactive Animation of Ocean Waves Symposium on Computer Animation,  pp. 161-1662002isbn 1-58113-573-4.

[Igl04] A. Iglesias Computer graphics for water modeling and rendering: a survey Future generation computer systems,  20 (2004)no. 81355—1374Novemberissn 0167-739X.

[LH62] Michael S. Longuet-Higgins The distribution of intervals between zeros of a stationary random function Philosophical Transactions for the Royal Society of London,  254 (1962)no. 1047557—599Series A, Mathematical and Physical Sciencesissn 0080-4614.

[MMS04] Viorel Mihalef Dimitris Metaxas, and Mark Sussman Animation and control of breaking waves Symposium on Computer Animation,  2004pp. 315—324isbn 3-905673-14-2.

[MTS75] Hisashi Mitsuyasu Fukuzo Tasai Toshiko Suhara Shinjiro Mizuno Makoto Ohkusu Tadao Honda, and Kunio Rikiishi Observations of the directional spectrum of ocean waves using a cloverleaf buoy Jour. of physical oceanography,  5 (1975)no. 4750—760issn 0022-3670.

[MWM87] Gary A. Mastin Peter A. Watterberg, and John F. Mareda Fourier synthesis of ocean scenes IEEE Computer Graphics and Applications,  7 (1987)no. 316—23Marchissn 0272-1716.

[PA01] Simon Premože and Michael Ashikhmin Rendering Natural Waters Computer Graphics Forum,  20 (2001)no. 4189—199issn 0167-7055.

[Pea86] Darwyn R. Peachey Modeling waves and surf SIGGRAPH Computer Graphics,  20 (1986)no. 465—74Augustissn 0097-8930.

[Phi57] Owen M. Phillips On the generation of waves by turbulent wind Journal of fluid mechanics,  2 (1957)417—445issn 0022-1120.

[PM64] Willard J. Pierson and Lionel Moskowitz A Proposed Spectral Form for Fully Developed Wind Seas Based on the Similarity Theory of S. A. Kitaigorodskii Journal of geophysical research,  69 (1964)5181—5203issn 0022-1406.

[Ran63] William John Macquorn Rankine On the exact form of waves near the surface of deep water Philosophical transactions of the Royal society of London,  153 (1863)127—138, issn 1471-2946.

[TCC84] M. J. Tucker P. G. Challenor, and D. J. T. Carter Numerical simulation of a random sea Applied ocean research,  6 (1984)no. 2118—122issn 0141-1187.

[Tes01] Jerry Tessendorf Simulating ocean water ACM SIGGRAPH course notes,  http://graphics.ucsd.edu/courses/rendering/2005/jdewall/tessendorf.pdf2001Updated in 2004.

[TG02] Sébastien Thon and Djamchid Ghazanfarpour Ocean waves synthesis and animation using real world information Computers and Graphics,  26 (2002)no. 199—108Februaryissn 0097-8493.

[vG04] Franz Joseph von Gerstner Theorie der Wellen Abhandlungen der Königlichen Böhmischen Gesellschaft der Wissenschaften,  18041—65.


Any party may pass on this Work by electronic means and make it available for download under the terms and conditions of the Digital Peer Publishing License. The text of the license may be accessed and retrieved at http://www.dipp.nrw.de/lizenzen/dppl/dppl/DPPL_v2_en_06-2004.html.