# Multifractal analysis of wind turbine power and rainfall from an operational wind farm -Part 1: Wind turbine power and the associated biases

Jerry Jose, Auguste Gires, Yelva Roustan, Ernani Schnorenberger, Ioulia Tchiguirinskaia, Daniel Schertzer

► To cite this version:

Jerry Jose, Auguste Gires, Yelva Roustan, Ernani Schnorenberger, Ioulia Tchiguirinskaia, et al. Multifractal analysis of wind turbine power and rainfall from an operational wind farm -Part 1: Wind turbine power and the associated biases. Nonlinear Processes in Geophysics, 2024, 31, pp.587-602. 10.5194/npg-31-587-2024. hal-04828797

HAL Id: hal-04828797

https://enpc.hal.science/hal-04828797v1

Submitted on 10 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Multifractal analysis of wind turbine power and rainfall from an operational wind farm – Part 1: Wind turbine power and the associated biases

Jerry Jose¹, Auguste Gires¹, Yelva Roustan³, Ernani Schnorenberger², Ioulia Tchiguirinskaia¹, and Daniel Schertzer¹,⁴,⁵

¹HM&Co, École nationale des ponts et chaussées, Institut Polytechnique de Paris, Champs-sur-Marne, France

²Boralex, Lyon, France

³CEREA, École nationale des ponts et chaussées, Institut Polytechnique de Paris, EDF R&D, Île-de-France, France

⁴Department of Civil and Environmental Engineering, Imperial College London, London, UK

⁵Department of Complexity Science, Potsdam Institute for Climate Impact Research, Potsdam, Germany

**Correspondence:** Jerry Jose (jerry.jose@enpc.fr)

Received: 31 January 2024 – Discussion started: 2 February 2024

Revised: 20 September 2024 – Accepted: 10 October 2024 – Published: 10 December 2024

**Abstract.** The inherent variability in atmospheric fields, which extends over a wide range of temporal and spatial scales, is also transferred to energy fields extracted from them. In the specific case of wind power generation, this can be seen in the theoretical power available for extraction and the empirical power produced by turbines. To model and analyse them, it is important to quantify their variability, intermittency, and correlations with other interacting fields across scales. To understand the uncertainties involved in power production, power outputs from four 2 MW turbines are analysed (from an operational wind farm at Pay d'Othe, 110 km south-east of Paris, France) using the scale-invariant framework of universal multifractals (UM). Their scaling properties were compared with power available at the same location from simultaneously measured wind velocity.

While statistically analysing the turbine output, the rated power acts like an upper threshold that results in biased estimators. This is identified and quantified here using the theoretical framework of UM and validated using numerical simulations. Understanding the effect of instrumental thresholds in statistical analysis is important in retrieving actual fields and modelling them, more so in wind power production, where the uncertainties due to turbulence are already a leading challenge. This is expanded in Part 2, where the influence of rainfall on power production is studied across scales using UM and joint multifractals.

## 1 Introduction

In the increasing global transition towards renewable and carbon-neutral energy, wind power is extremely attractive, as it has some of the lowest carbon emissions in life cycle assessment (Li et al., 2020; Guezuraga et al., 2012; Wiser et al., 2011). The levelized cost of energy (LCOE, including building and operation) has also decreased drastically in recent decades for both offshore and onshore wind power (80% since early 1980 and a further 30% in the past 5 years), giv-

ing it better economic value (Beiter et al., 2021). However, wind is a fluctuating field and mainly owes its generation to uneven heating of Earth's surface by solar radiation and the pressure gradients generated from it. Further, atmospheric turbulence makes the characterization of the field a difficult task (with the governing Navier-Stokes equations still remaining unsolvable). The small-scale fluctuations and intermittency in wind are transferred to the power produced; this is further complicated by the fact that wind turbine hubs are located in the atmospheric boundary layer. In addition, an im-
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Log-log plot of the exceedance probability, $Pr(\Delta v > s)$, of only positive horizontal velocity increments, $\Delta v(\tau) = v(t + \tau) - v(t)$, from December 2020 to July 2021, at 1 Hz from location 1 of the RW-Turb meteorological mast along with a Gaussian distribution to illustrate the latter's inadequacy. $s$ is a threshold of intensity, and $\tau$ here is 15 s.</figcaption></figure>

proved understanding of turbulence is identified as one of the leading challenges in the field of wind power by experts (van Kuik et al., 2016). When it comes to the working of modern turbines, one way to account for wind variations is through variable speed turbines and adaptive torque control enabling maximum power capture. However, the commonly used parameter for control, “turbulence intensity” (standard deviation of wind speed divided by mean wind speed over 10 min), cannot fully capture the behaviour (see the non-Gaussian behaviour of wind velocity in Fig. 1) and is too coarse to represent the variability (active torque controls should be responsive down to a few seconds). Further, this does not consider any effect of rain that could get transferred to loads on turbines (Johnson, 2004).

To understand the complex effect of turbulence on power production, along with access to high-resolution data, an appropriate theoretical framework is required to characterize intermittency at all scales of measurement. The scale-invariant multifractal framework of universal multifractals (UM), which is widely used to study variability in geophysical fields, can be used to characterize this complexity (Schertzer and Lovejoy, 1987, 1997). Using the framework of UM, Fitton et al. (2011, 2014) studied scaling behaviour and multifractal properties of wind velocity and torque fluctuations at wind farm test sites (in Germany and Corsica) and made a case for multifractal modelling of atmospheric turbulence. Multi-fractality of wind speed and aggregate wind farm power was illustrated in Calif and Schmitt (2014), where the coupling between both fields was examined.

In light of the scientific perspectives mentioned so far, here we try to characterize the small-scale fluctuations in wind power production using data from an operational wind farm at Pays d'Othe, 110 km south-east of Paris, France. However,

while analysing the variability of the field using statistical methods, the presence of instrumental limits in the data can introduce biases. For example, the effect of an instrumental lower threshold is discussed in Jose et al. (2021) in the framework of UM analysis using an atmospheric extinction coefficient ($\sigma_e$) as the field. Similarly, there is also bias from the presence of zeroes in the data (Gires et al., 2014). Both of these biases are present in statistical analyses of empirical power from wind turbines since turbines are designed to work at a rated power (here, 2 MW) and can provide zero or negative power (more consumption than production). The major aim of this paper is to highlight these biases encountered during multifractal analysis and their influence on the direct statistical analysis of turbine power. For the theoretical aspect, only the effect of the upper threshold, which has not yet been addressed in the literature, is considered here to avoid complexity. More analysis is intended for a follow-up paper, where the influence of rain on wind power production is examined along with the coupling of power as a field with other atmospheric parameters. The details of the data collection and quality are presented in the following section on data and methods; the second part of this section briefly recapitulates the framework of UM. The biases encountered in the analysis of turbine power are presented in the section that follows, along with numerical simulations where it is identified and reproduced in the framework of UM. Acknowledging these biases, some efforts were made to characterize the effect of rainfall and wind velocity on turbine power. The final section concludes the study and summarizes the results.

# 2 Data and methods

## 2.1 Data and instrumentation

The Rainfall Wind Turbine or Turbulence project (RW-Turb; https://hmco.enpc.fr/portfolio-archive/rw-turb/, last access: 26 November 2024), supported by the Agence Nationale de la Recherche (ANR, French National Research Agency in English), is designed to understand the long- and short-term effects of rainfall on wind power production, with simultaneous high-resolution measurements on an operational wind farm. Interested readers are directed to Gires et al. (2022) for an overview of the campaign. To briefly summarize this, the RW-Turb measurement campaign (at Pay d'Othe, 110 km south-east of Paris, France) consists of a meteorological mast (Fig. 2b on the right side) on a wind farm (jointly operated by Boralex, https://www.boralex.com/our-projects-and-sites/, last access: 26 November 2024), and JP Énergie Environnement, https://pays-othe-89.parc-eolien-jpee.fr/, last access: 26 November 2024). Figure 2a shows the location of the project; the nine wind turbines of the Pays d'Othe wind farm (aligned south-east and within a 4 km radius) are marked as black vertical crosses and the meteorological mast as a star (in the middle). Data from four Vestas V-90 (marked 1, 2, 8,
and 9 in Fig. 2a) instruments are available: two are closer to the meteorological mast and two are farther from it (≈ 3.5 km from the mast). The five turbines of the Molinons wind farm in the north are also visible within a 5 km radius (grey vertical crosses). It should also be noted that a small groove is located just south of the mast at roughly 160 m; a larger one is to the east at roughly 100 m. Nearby the mast (i.e. within the 1 km radius), there is a small slope in the north-south direction. The meteorological mast consists of two sets of optical disdrometers (OTT Parsivel², 30 s, not used in current study), 3D sonic anemometers (ThiesCLIMA, 100 Hz), and a mini meteorological station (1 Hz) at heights of roughly 45 and 80 m (managed by the Hydrology Meteorology and Complexity laboratory of the École des Ponts ParisTech – HM&Co, ENPC).

Technical and working information on the turbine can be found in Vestas Wind Systems A/S (2023). The turbines have a rated power of 2.0 MW, which is pitch-regulated with variable speed. The hub height of the turbines is 80 m, which is closer to the vertical height of the upper set of devices on the mast (≈ 78 m). The turbines have a cut-in wind speed of 4 m s⁻¹ and a rated wind speed of 12 m s⁻¹. This can be seen in the power curves in Fig. 3 (last row), where the turbine register power at the cut-in speed maintains the rated power of 2000 kW after the rated wind speed. The cut-out speed of the turbine is 25 m s⁻¹ (the extreme x-axis point of the power curves): this is the speed at which a turbine stops registering power. Generally, the turbines register positive values of wind power. However, when the power retrieved from wind is less than what is required to work a turbine, it registers negative power. This can be seen in the power curves as clusters around 0. Along with the wind power, the turbine also provides information on the local velocity (at a sampling measurement rate of 15 s), which is used for internal regulation; this is used to plot the power curves in Fig. 3. The wind power data used for the studies come from the four turbines of Boralex – 1 and 2 are located closest to the mast (as can be seen in Fig. 2b), and 8 and 9 are located at the farthest end (at a sampling frequency of 15 s).

The temporal evolution and power curves (power vs. velocity, with the expected curve provided by the manufacturer in red) for the turbines are shown in Fig. 3 for 3 months (from 1 January to 1 April 2021). There are instances where the turbine failed to produce any power and had to consume energy for its basic operation. This results in negative data values, and for realistic analysis they were considered to be zero. This is why there is clustering of points at zero in the power curve (Fig. 3e-h).

In addition to the empirical power provided by turbines, the theoretical power available for extraction can be obtained by

$$P_a = \frac{1}{2} \rho A v^3 C_p, \qquad (1)$$

where ρ is the air density at the wind turbine height (h hub), A the swept area of the turbine rotor, v the wind velocity (m s⁻¹) approximated at the turbine height, and C p the power coefficient or Betz coefficient (for the Vestas V-90 instruments examined here, h hub = 80 m, A = 6362 m², the rated power is 2 MW, and C p is 0.593). Here, P a is estimated at the same sampling rate as that of P t (15 s) despite the 3D sonic anemometer and mini meteorological station registering data at finer sampling rates.

The value of air density is often approximated as 1.255 kg m⁻³ (standard value at sea level 15 °C). However, it is known to show fluctuations and is reported to have an effect on power generation at varying levels (Jung and Schindler, 2019; Ulazia et al., 2018). For the purpose of this analysis, air density was considered a varying quantity and estimated using the current official formula of the International Committee for Weights and Measures (CIPM), referred to as the CIPM-2007 equation, which accounts for humidity (Picard et al., 2008):

$$\rho(T, P, H_r) = \frac{P M_a}{Z(T, P, H_r) R T(K)} \times \left\{1 - x_v(T, P, H_r)\left[1 - \frac{M_v}{M_a}\right]\right\}, \qquad (2)$$

where T (°C), P (Pa), and H r (0 ≤ H r ≤ 1) are the temperature, pressure, and humidity from the meteorological station at h hub. The other derived parameters are the following:

- T (K) the air temperature (K; from T);

- Z the compressibility factor (a function of T and P);

- R the molar gas constant (J mol⁻¹ K⁻¹);

- x v the mole fraction of water vapour;

- M a the molar mass of dry air (g mol⁻¹); and

- M v the molar mass of water (g mol⁻¹).

## 2.2 Scaling analysis and UM framework

Spectral analysis is widely used for characterizing scaling properties; here, the second-order statistics of rain in the frequency domain were examined for power law scaling as follows (Mandelbrot, 1982; Schertzer and Lovejoy, 1985):

$$E(k) \approx k^{-\beta}, \qquad (3)$$

where k corresponds to the wave number and β is the spectral exponent.

However, to fully characterize the complexity of the process across its intensities and spatio-temporal variation, information on higher- and lower-order statistics is required. For this, we use UM, which rely on the assumption of the field being generated by an underlying cascade process with
<figure><img src="image_5.png" /><figcaption>Figure 2. (a) Map of the Pays d'Othe wind farm (inset: location in France): the meteorological mast is at the centre, and the available turbines are numbered – 1, 2, 8, and 9. (b) Turbine 1 and the mast. (c) Turbines 1 and 2 as seen from the bottom of the mast. The figures are adapted from Gires et al. (2022).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3. Time series of the empirical turbine power (first row) and the power vs. velocity plot with the theoretical state curve of the turbine (second row) for the wind turbines (1, 2, 8, and 9) at Pays d'Othe.</figcaption></figure>

conserved statistical properties at each scale while inher-
iting the scale-invariant properties of Navier–Stokes equa-
tions (Schertzer and Lovejoy, 1987, 1989; Schertzer and
Tchiguirinskaia, 2020). In this framework, the probability of
a field exceeding a particular threshold across all the scales is
captured using the scale-invariant notion of singularity (γ),
and for a multifractal field this probability scales accord-

ing to the resolution (λ: the ratio of L, the outer scale, to
l, the observational scale), with a corresponding fractal co-
dimension as the scaling exponent c(γ) (Schertzer and Love-
joy, 1987, 1988):

$$
p(\epsilon_{\lambda} \geq \lambda^{\gamma}) \approx \lambda^{-c(\gamma)}. \quad (4)
$$
This relation implies that the statistical moments $q$ of the field also scale with resolution and the moment scaling function $K(q)$ as (Schertzer and Lovejoy, 1987, 1988)

$$ \langle \epsilon_{\lambda}^{q} \rangle \approx \lambda^{K(q)}. \qquad (5) $$

Both functions are related by Legendre transform (Parisi and Frisch, 1985). For a conservative field in the UM framework, $K_c(q)$ can be fully determined with only two parameters, the multifractality index $\alpha$ and the mean intermittency co-dimension $C_1$.

$$ K_c(q) = \begin{cases} \frac{C_1}{\alpha - 1} (q^{\alpha} - q) & \alpha \neq 1 \\ C_1 q \ln q & \alpha = 1 \end{cases} \qquad (6) $$

$C_1$ measures the clustering of average intensity across scales ($C_1 \in [0, 1]$ for one-dimensional fields); when $C_1 = 0$, the field is homogeneous with little variability. $\alpha$ measures how this clustering changes with respect to intensity levels ($\alpha \in [0, 2]$): the higher the value of $\alpha$, the higher the variability, with $\alpha = 0$ being a mono-fractal field where intermittency of the extreme is the same as that of the mean. If the UM parameters are known, the co-dimension function of the conservative multifractal field $c_c(\gamma)$ can also be obtained using Lengendre transform (Schertzer and Lovejoy, 1987, 1988; Parisi and Frisch, 1985):

$$ c_c(\gamma) = \begin{cases} C_1 \left( \frac{\gamma}{C_1 \alpha'} + \frac{1}{\alpha} \right)^{\alpha'} & \alpha \neq 1, \\ C_1 \exp \left( \frac{\gamma}{C_1} - 1 \right) & \alpha = 1, \end{cases} \qquad (7) $$

where $\frac{1}{\alpha} + \frac{1}{\alpha'} = 1$.

For a non-conservative field $\phi_{\lambda}$, i.e. a field whose average $\langle (\phi_{\lambda}) \rangle$ changes with scales, a non-conservative parameter $H$ (not to be confused with the Hurst exponent: though for values $> 0$ both quantify long-range correlations, the latter does not have a simple general expression for the multifractal process; see the Appendix A for more details) is used in the expression of scaling (Schertzer and Lovejoy, 1987, 1988; Lovejoy and Schertzer, 2013):

$$ \phi_{\lambda} = \epsilon^{d} \lambda^{H}, \qquad (8) $$

where $\epsilon = d$ denotes equality in the distribution ($X = dY \Leftrightarrow \forall x$: $\Pr(X > x) = \Pr(Y > x)$) and $\epsilon$ is a conservative field characterized by $C_1$ and $\alpha$ for a conservative field $H = 0$. For a non-conservative field with positive $H$, fractional differentiation is required to retrieve a coarser field. Similarly, from a non-conservative field with a negative value of $H$, the conservative field is retrieved through fractional integration. $H$ is related to the spectral slope $\beta$ (Eq. 3) via the relationship (Tessier et al., 1993)

$$ \beta = 1 + 2H - K_c(2). \qquad (9) $$

The scaling behaviour of conservative multifractal fields can be examined using a trace moment (TM), where a log-log plot of upscaled fields against resolution $\lambda$ is taken for each moment $q$ (Eq. 5). The quality of the scaling is given by the estimate $r^2$ of the linear regression; the value for $q = 1.5$ is used as a reference. A double trace moment (DTM) is a more robust version of a TM tailored for UM fields, where the moment scaling function $K(q, \eta)$ of the field $\epsilon_{\lambda}^{(\eta)}$ (the initial field raised to power $\eta$ at maximum resolution and renormalized) is expressed as a function of the multifractal index $\alpha$ (Lavallée et al., 1993).

$$ \langle (\epsilon_{\lambda}^{(\eta)})^{q} \rangle \approx \lambda^{K(q, \eta)} = \lambda^{\eta \alpha} K(q) \qquad (10) $$

From the above equation, a value of $\alpha$ can be obtained as the slope of the linear part, where $K(q, \eta)$ is represented for a given $q$ as a function of $\eta$ in the log-log plot. Both the TM and DTM techniques give reliable estimates as long as $H < 0.5$ for the studied field.

Since multifractal processes are generated by cascade processes, the average values can become too concentrated over a certain area, leading to spurious estimates of moments above a particular value of $q$ (at $q_D$, $q$ is above $K(q) \approx +\infty$). This effect is called divergence of moments. The convex na-tures of the functions $K(q)$ and $c(\gamma)$ are also limited by the sample size of the data or the maximum value of the scale-invariant threshold or singularity ($\gamma_s$) and the corresponding moment ($q_s$). The details of their computation can be found in Schertzer and Lovejoy (1992), Schertzer and Lovejoy (1989), and Lovejoy and Schertzer (2007). For reliable statistical estimates of the moment scaling function and hence the UM parameters, the moment orders used should not exceed $q_s$ (moment order corresponding to the maximum singularity) or $q_D$ (moment order where divergence happens).

An in-depth discussion of the methodological choices we have made, along with an overview of other multifractal formalisms and their strengths and weaknesses, is provided in Appendix A.

# 3 Turbine power, biases, and associated issues in data analysis

## 3.1 Turbine power and biases

For the Vestas V-90 instruments, the rated power is 2 MW; this means that the maximum power the turbine can produce is 2000 kW. However, if we calculate the available power as per Eq. (1), there are many instances where it can go beyond the rated value (see Fig. 4c). When analysing the variability of a field using statistical methods, the presence of instrumental limits in the data (here an upper limit) can introduce biases. As briefly mentioned in the Introduction section, an instrumental lower threshold in the data can increase $\alpha$ and decrease $C_1$ (Jose et al., 2021). In addition to this, there is also bias from the presence of zeroes in the data (Gires et al., 2014) (underestimation of $\alpha$, deterioration of scaling), which
replaced negative values of the turbine power. Figure 4 shows the real and theoretical turbine power state curves along with the bias this poses in statistical analysis. Along with the power produced, the turbine data provide wind velocity at the hub (from a basic sensor installed on the hub), which is used for its internal monitoring. For research purposes, the wind velocity from the 3D anemometer at the mast is more desirable as it offers a more reliable measurement (on almost the same horizontal plane as the turbine hub). Figure 4a shows the empirical state curve of the turbine with this velocity, and Fig. 4b shows the same curve with the velocity from the anemometer. There is considerably more scatter with the latter. It should be stated that the turbines are not in the exact location of the mast (Turbines 8 and 9 are ≈ 3.5 km away), and hence approximation of the wind velocity from the mast (for computing the theoretical power in Eq. 1) comes with some biases.

From Fig. 4c and d, it can be clearly seen that the rated power imposes an upper threshold on the turbine power ($P_t$), while the power available ($P_a$) is the actual underlying field. For this 1-week-long series of $P_t$, 21.7% of the data were at an upper threshold and 2.9% were either zero or negative (taken as zeroes in the analysis); this percentage was found to change according to the data selected. The effect of these limits on UM analysis is shown in Fig. 5a, where the data in Fig. 4 are treated as an ensemble of 32 min. UM analysis was performed on direct fields as the values of $H$ were within acceptable limits ($H < 0.3$). A unique scaling regime from 15 s to 32 min was considered. The presence of the rated power clips the values of the field and results in a reduced value of $\alpha$ for $P_t$ (Fig. 5a: $\alpha = 1.36$, $C_1 = 0.00715$) from that of $P_a$ (Fig. 5b: $\alpha = 1.93$, $C_1 = 0.01753$). Imposition of a similar threshold ($P_a \le 2000 = 2000$) on $P_a$ was found to artificially reduce the estimates ($\alpha$ from 1.93 to 1.39, $C_1$ from 0.10753 to 0.0076) in Fig. 5c, bringing them closer to that of biased turbine power, $P_t$ (Fig. 5a). Even closer values of $\alpha$ were obtained when a lower threshold was also imposed (replacing $P_a$ values with zeroes at positions where $P_t$ was negative), giving $\alpha$ a value of 1.35 and $C_1$ a value of 0.0077 (Fig. 5d). The results are compiled in Table 1. The scaling quality remained similar for all the cases mentioned here, with the $r^2$ value (of the TM curve at $q = 1.5$) remaining around 0.99 (second column of Fig. 5).

It should be noted that the effect of a threshold could be different according to the size of the sample and the scaling regimes studied. In the same spirit as $\sigma_e$ in Jose et al. (2021), the effect of rated power as an upper threshold on $P_t$ is explored here in the theoretical framework of UM. The effect on different scaling regimes as well as the additional complexity from the known effect of zeroes (Gires et al., 2012), though identified here, are not considered to avoid complexity.

<table><caption>**Table 1.** Values of UM parameters for (a) the power produced by the turbines ($P_t$), which has intrinsic thresholds (21.7% as the upper threshold and 2.9% as the lower threshold); (b) the power available ($P_a$) without any thresholds; (c) $P_a$ with an imposed upper threshold; and (d) $P_a$ with an imposed upper threshold and a lower threshold. The graphs can be seen in Fig. 4.</caption><thead><tr><th rowspan="2">Field</th><th rowspan="2">Threshold</th><th colspan="4">UM parameters</th></tr><tr><th>α<sub>DTM</sub></th><th>C<sub>1,DTM</sub></th><th>β</th><th>H</th></tr></thead><tbody><tr><td>P<sub>t</sub></td><td>Upper + lower</td><td>1.36</td><td>0.00715</td><td>1.6</td><td>0.31</td></tr><tr><td>P<sub>a</sub></td><td>Upper</td><td>1.93</td><td>0.00753</td><td>1.39</td><td>0.21</td></tr><tr><td>P<sub>a</sub></td><td>Lower</td><td>1.35</td><td>0.0077</td><td>1.44</td><td>0.23</td></tr><tr><td>P<sub>a</sub></td><td>Upper + lower</td><td>1.35</td><td>0.0076</td><td>1.44</td><td>0.23</td></tr></tbody></table>

## 3.2 Understanding the effect of upper threshold in the UM framework

Let us take the upper threshold (rated power in this case) at the largest possible scale ratio as

$$T = \Lambda^{\gamma_T}, \qquad (11)$$

where $\gamma_T$ is the singularity corresponding to threshold $T$ and $\Lambda$ is the maximum resolution (length of the time series). For multifractal fields, the probabilities of exceeding scale-independent thresholds, $\lambda^\gamma$, scale with resolution $\lambda$ (see Eq. 4). At the upper threshold $T$, this yields

$$\mathrm{Pr}(\epsilon_{\lambda} \ge T) \approx \lambda^{-c(\gamma_T)}. \qquad (12)$$

If we set the upper threshold, i.e. set all the values of the field greater than $T$ equal to $T$ (represented here by this expression: $\{\epsilon_{\lambda} \ge T\} = T$), the probability of having values greater than $T$, $\mathrm{Pr}(\epsilon_{\lambda} > T)$, becomes 0, reducing the above relation into $\mathrm{Pr}(\epsilon_{\lambda} = T) \approx \lambda^{-c(\gamma_T)}$. This leaves the value of $c(\gamma)$ equal to $+\infty$ for singularities above $\gamma_T$ (for $\gamma > \gamma_T$, $c(\gamma) = +\infty$). Here, $c(\gamma_T)$ is the limiting non-zero value above which $c(\gamma)$ becomes $+\infty$.

This effect of the upper threshold ($c(\gamma) \to +\infty$ for $\gamma > \gamma_T$) is similar to the effect of the sampling dimension ($D_s$) in the UM framework. The maximum observable singularity can be defined by taking the probability at the corresponding threshold as in Eq. (12):

$$\mathrm{Pr}(\epsilon_{\lambda} \ge \lambda^{\gamma_s}) \approx \frac{1}{N_s \lambda^D}, \qquad (13)$$

where $N_s$ is the number of samples and $\lambda^D$ is the number of values per sample. $N_s = \lambda^{D_s}$ ($D_s$ being the sampling dimension, which quantifies the number of independent samples with resolution $\lambda$; for one sample, $D_s = 0$). Using the notions of $D_s$ and $D$, with $\gamma$ corresponding to the sampling resolution, $\gamma_s$ can be estimated from $c(\gamma_s) = D + D_s$ (Hubert and Carbonnel, 1989; Lovejoy and Schertzer, 2007).

The moment scaling function $K(q)$ and co-dimension function $c(\gamma)$ were discussed earlier in terms of UM parameters in Eqs. (6) and (7). Both are equivalent functions, and for
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 4.** Illustration of the upper threshold (by virtue of a rated power of 2000 kW) in power produced by turbines. (a) Empirical and theoretical power state curve of Turbine 1 with wind velocity from the turbine. (b) Wind velocity from location 1 on the mast. (c) Power produced by the turbine ($P_t$) and actual wind power available $P_a$. (d) Effect of the rated power as a threshold of the time series and effect of negative values on $P_t$ for 1-week-long data – 20 to 26 May 2021.</figcaption></figure>

multifractals they are related by a simple Legendre transform (Parisì and Frisch, 1985; Schertzer and Lovejoy, 1993):

$$ 
\begin{aligned}
K(q) &= \max_{\gamma}[q\gamma - c(\gamma)], \\
c(\gamma) &= \max_{q}[q\gamma - K(q)].
\end{aligned} 
\quad (14) $$

Hence, for every singularity $\gamma$, there is a corresponding order of moment $q$ associated with it and vice versa: $q = c'(\gamma_q)$ and $\gamma = K'(q_\gamma)$.

When $\gamma > \gamma_s$, $c(\gamma) = +\infty$; by Legendre transform $K(q)$ becomes linear from $q > q_s = c'(\gamma_s)$:

$$ 
\begin{aligned}
\gamma_s &= \alpha' C_1 \left( \frac{D + D_s}{C_1} \right)^{\frac{1}{\alpha'}} - \frac{C_1}{\alpha - 1}, \\
q_s &= \left( \frac{D + D_s}{C_1} \right)^{\frac{1}{\alpha}}.
\end{aligned} 
\quad (15) $$

In the case of the sampling dimension, $c(\gamma)$ varies as follows:

$$ c(\gamma) = \begin{cases} +\infty & \text{for } \gamma > \gamma_s \\ D + D_s & \text{for } \gamma = \gamma_s \\ c(\gamma) & \text{for } \gamma < \gamma_s \end{cases}. \quad (16) $$

Similarly, in the presence of the upper threshold here ($\{\epsilon_\lambda \ge T\} = T$), $c(\gamma)$ reaches $+\infty$ at an earlier limiting value $c(\gamma_T)$, where $\gamma_T < \gamma_s$ (Fig. 6a).

$$ c_T(\gamma) = \begin{cases} +\infty & \text{for } \gamma > \gamma_T \\ c(\gamma_T) & \text{for } \gamma = \gamma_T \\ c(\gamma) & \text{for } \gamma < \gamma_T \end{cases} \quad (17) $$

Here, $\gamma_T$ is defined from the threshold as initially stated in Eq. (11), and $c_T(\gamma)$ is estimated as above. From this, the corresponding limit moment $q_T$ can be obtained as in Eq. (15).

$$ q_T = \left( \frac{c(\gamma_T)}{C_1} \right)^{\frac{1}{\alpha}} \quad (18) $$

To summarize, in standard data analysis with sampling limitations, $c(\gamma)$ is bounded by a maximum value $c(\gamma_s)$, above which it becomes infinite. $K(q)$, which is connected to $c(\gamma)$ through Legendre transform (Eq. 14), becomes linear beyond this $q$ ($q \ge q_s$) value ($K(q) = (q - q_s)\gamma_s + K(q_s)$). Similarly, in this specific case, when an upper threshold is imposed ($\{\epsilon_\lambda \ge T\} = T$), $K(q)$ becomes linear at an earlier value of
<figure><img src="image_11.png" /><figcaption>Figure 5. Spectral analysis (Eq. 3), TM analysis (Eq. 5), DTM curve (Eq. 10), and K(q) for (a) power produced by a turbine ($P_t$) with intrinsic thresholds (upper: due to rated power; lower: due to negative values that are treated as zeroes); (b) power available ($P_a$), which is the unbiased actual field; (c) $P_a$, where an upper threshold is imposed at the rated power of the turbine (i.e. all values of $P_a$ above the rated power of the turbine are artificially replaced by 2 MW; values > rated power = rated power); and (d) $P_a$, where an upper threshold and a lower threshold (values set to zero, where $P_t < 0$) are imposed based on the turbine values. Data used: time series from 20 to 26 May 2021, with the lowest time step of 15 s and a sample size of 32 min.</figcaption></figure>

$q$ (at $q_T < q_s$) as defined by $\gamma_T$ (Fig. 6b).

$$K_T(q) = \begin{cases} \gamma_T(q - q_T) + K(q_T) & \text{for } q > q_T \\ K(q_T) = q_T\gamma_T - c(\gamma_T) & \text{for } q = q_T \\ K(q) & \text{for } q < q_T \end{cases} \quad (19)$$

Using the DTM technique, for a given $q$, $K(q, \eta) = K(q\eta) - qK(\eta)$, which for UM fields is $\eta^\alpha K(q)$. When no thresholds are applied, $K(q, \eta)$ varies as

$$K(q, \eta) = \begin{cases} (q-1)(D + D_s) & \text{for } \eta \ge \eta_+(q) \\ \eta^\alpha K(q) & \text{for } \eta < \eta_-(q) \end{cases}, \quad (20)$$

where $\eta_+(q)$ corresponds to the maximum values of $\eta$ above which $K(q, \eta)$ becomes a plateau due to sampling limitations (Eq. 15). To elaborate, $K(q, \eta)$ consists of two parts

$K(q\eta)$ and $K(\eta)$, and $\eta_+(q)$ corresponds to the value of $\eta$ above which both are linear (which is $q_s$). The transition to the plateau starts at the lower value $\eta_-(q)$ (which is $q_s/q$), above which only $K(q\eta)$ is linear. In the presence of an upper threshold ($\{\epsilon_\lambda \ge T\} = T$), the DTM curve will be (Fig. 6c)

$$K_T(q, \eta) = \begin{cases} (q-1)c(\gamma_T) & \text{for } \eta \ge \eta_+(q) \\ \eta^\alpha K(q) & \text{for } \eta < \eta_-(q) \end{cases}, \quad (21)$$

where $\eta_+(q) = q_T$ and $\eta_-(q) = q_T/q$.

It is important to note here that the value of $K(q, \eta)$ does not reach the upper plateau abruptly at $q_T$ or $q_s$. Rather, it flattens gradually, starting from a value of $\eta = q_s/q$ or $q_T/q$ (as per the value of $q\eta$ in $K(q, \eta)$). The presence of the upper threshold shifts this starting point and decreases the range of
<figure><img src="image_12.png" /><figcaption>**Figure 6.** Influence of the threshold on (a) the $c(\gamma)-\gamma$ curve, with $c(\gamma)$ reaching $+\infty$ at $\gamma_T$ and then $\gamma_S$; (b) the DTM curve, with $K(q)$ becoming linear at $q_T$ and then $q_S$; and (c) the $K(q)-q$ curve, with $K(q, \eta)$ reaching the upper plateau early. Arbitrary values were used for $\gamma_S$ and $\gamma_T$; UM parameter values of the fields were taken to be $\alpha = 1.8$ and $C_1 = 0.2$.</figcaption></figure>

possible values for estimation of $\alpha$ (slope of the DTM curve). Hence, the presence of a plateau will result in biased (reduced) estimates.

## 3.3 Numerical simulations

Underestimation of values of $\alpha$ due to application of the upper threshold was already observed in Fig. 5c using real data. To understand this further, conservative multifractal fields ($H = 0$) were simulated using discrete cascades with values of UM parameters close to those observed for empirical powers ($P_t$). Discrete cascade simulation here involves division of a parent structure into sub-structures (retaining the value of the parent structure multiplied by a random factor; Chambers et al., 1976), iteratively following a non-infinitesimal scale ratio while maintaining the validity of Eqs. (5) and (6).

For ease of contrast with the simulations examining lower thresholds in Jose et al. (2021), values of $\alpha = 1.8$ and $C_1 = 0.2$ were used. UM analysis was implemented in ensembles (sample size 128 and number of samples 100) by progressively applying the upper threshold till the percentages were observed in $P_t$ (~30%). $K(q)$ becomes linear at earlier and earlier values of $q$ (after $q_T$), and an increasing percentage of values at the threshold can be seen in the third column (like in Fig. 6b). The DTM curve in the second column shows that both $\alpha$ and $C_1$ decrease with progressive application of the thresholds (from 0% to 30%: $\alpha$ decreased from 1.8 to 1.56, while $C_1$ decreased from 0.17 to 0.05).

While discussing this bias in the framework before, it was mentioned that the upper threshold was introduced at the maximum resolution (Eq. 11, $\Lambda^{\gamma_T}$). Since in practice the lower scales in UM are obtained from averaging the outer scale (at the maximum resolution), the threshold values (and hence $\gamma_T$) at each stage do not correspond exactly to the originally defined one. This is the reason for an increased “transition part” (part of the curve from the straight line to the upper plateau) of the DTM curves (second column) in the simulations here (more than that in Fig. 6c). When the esti-

<table><caption>**Table 2.** Values of the UM parameters for the simulated fields with an artificial imposition of the upper thresholds.</caption><thead><tr><th rowspan="2">Percentage at the threshold</th><th colspan="4">UM parameters</th></tr><tr><th>α<sub>DTM</sub></th><th>C<sub>1,DTM</sub></th><th>β</th><th>H</th></tr></thead><tbody><tr><td>0%</td><td>1.81</td><td>0.17</td><td>0.95</td><td>0.108</td></tr><tr><td>5%</td><td>1.72</td><td>0.13</td><td>1.15</td><td>0.186</td></tr><tr><td>15%</td><td>1.64</td><td>0.10</td><td>1.18</td><td>0.178</td></tr><tr><td>30%</td><td>1.56</td><td>0.07</td><td>1.2</td><td>0.163</td></tr></tbody></table>

mation of $\alpha$ was forced at $\eta = 1$ (so that the TM and DTM estimates were the same), the bias in the values of $\alpha$ increased as the slope estimation moved to the transition part. For example, the already biased value of $\alpha$ at the 30% threshold, 1.58 (slope at $\log \eta$ between -0.1 and -0.5), was further reduced to 0.95 (slope around $\log \eta = 0$). In this estimation, $C_1$ moreover remained similar to previous estimates at all the thresholds.

It is interesting to note that the trend here (only for $\alpha$) is not exactly the opposite of what was observed during numerical simulations with lower thresholds in the section (Jose et al., 2021). While the imposition of a lower threshold increased $\alpha$ and decreased $C_1$, the upper threshold here reduces both UM parameters. In the specific case of turbine power, $P_t$ (Fig. 5a) has a combination of the upper threshold from the rated power and the lower threshold (zeroes) from a negative power (the latter is not considered here). This, in practice, further reduces the range of available $\eta$ for estimation of $\alpha$ by imposing a lower plateau as well (see Fig. 5, third column). Also, the effect of this bias could be different when fluctuations of the fields are selected to retrieve conservative fields since the simulations were performed directly on conservative multifractal fields here. Since two consecutive power values can be the same, thanks to the rated power, taking fluctuations will yield zeroes in the field and add to the zero bias. The effect of both biases could be different when
<figure><img src="image_13.png" /><figcaption>(b) % of values at threshold = 5 %</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c) % of values at threshold = 15 %</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d) % of values at threshold = 30 %</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 7.** Effect of the upper thresholds illustrated using numerical simulations – discrete cascades of size 128 with 100 samples with $\alpha = 1.8$ and $C_1 = 0.2$ as input. The thresholds were applied progressively to the simulated field: (a) no threshold, (b) 5%, (c) 15%, and (d) 30%. A decrease in $\alpha$ and an increase in $C_1$ with a threshold can be seen from the DTM curves (second column) of panels (a) to (d).</figcaption></figure>

the aggregate power of the wind farm is considered as well: this is not explored here either.

## 3.4 Data analysis reducing the biases

So far, the effect of thresholds in UM analysis has been illustrated in the framework of UM. To get a better idea of their effect on scaling, UM analysis was performed of a longer series of powers. For this, a 3-month-long continuous series was taken (from 1 January to 1 April 2021) and UM analysis was performed on $P_a$ and $P_t$ as ensembles of sample size 128 (32 min). Figure 8a and b show the curves for $P_a$ ($\alpha$: 1.93; $C_1$: 0.01) and $P_t$ ($\alpha$: 1.11; $C_1$: 0.0042) respectively. Considering $P_a$ to be the underlying field, the effects of the thresholds (upper and due to zeroes) in $P_t$ can be seen in the DTM curve (Fig. 8b, second row). The lower plateau corresponds

to the presence of a negative power in the data (which was replaced by zeroes), and the low value of $\alpha$ is due to $\eta$ being in the transition range (as already seen in Fig. 6c). Figure 8c and d show the UM analysis for the same data but remove the columns with thresholds. In the $P_t$ ensemble data, the columns with thresholds, zeroes, and repetition of data were removed by using a limit of 0.01%. For example, columns with more than 0.01% of their values $\ge$ 1600 were removed to be on the safe side when analysing data without the influence of a threshold. For a more accurate comparison, the same columns were removed from $P_a$ as well; the results are shown in Fig. 8c ($\alpha$: 1.76; $C_1$: 0.0095) and Fig. 8d. It can be seen that the lower plateau disappeared for $P_t$ and that the values of the UM parameters ($\alpha$: 1.5749 from 1.11; $C_1$: 0.00554 from 0.0042) improved; this also increased the value
of β (1.8 from 1.6) and consequently that of H (0.41 from 0.3). It should be noted that the values of the UM parameters get closer to those of P_a but are not the exact values. This suggests that there are slight differences in the properties of both fields, even though they appear comparable when biases are removed.

## 4 Conclusion

Wind turbines are designed to work at a rated power for optimal production of power as well as their safe functioning. This inadvertently creates an upper threshold in the output data of power production, and such an effect induces biases in statistical analysis, especially when the small-scale nonlinear variability and intermittency are to be studied. The backdrop of this study followed the campaign in Gires et al. (2022), where the main objectives were to analyse the turbine power P_t as a temporal field and to gain insights into its correlation with rainfall, which is poorly understood, and also with other meteorological fields, across scales using the data averaged to these reliable frequencies. However, direct analysis of empirical turbine power (using the framework of UM) was found to be difficult since the output from wind turbines is limited by a maximum or rated power. In time series analysis this acts as an upper threshold, resulting in reduced estimates of UM parameters compared to those of theoretically available wind power (P_a) for extraction. The reason for this decrease was identified in the framework of UM and is illustrated using theoretical formulations. The same has been confirmed through simulations of conservative multifractal fields as well. Basically, the presence of an upper threshold introduces an upper plateau into the DTM curve, which is similar to the one due to the sampling dimension, but it begins at a lower value of η. This reduces the range of available η values for estimation of the slope and hence results in a biased value of α (reduced α and C_1). Also, when it comes to the empirical power produced by the turbines, the biases are two-fold since a lower threshold (albeit to a much lesser extent) is involved, because the turbine does not necessarily always produce power and has moments that involve only operational consumption (leaving power production values negative). Since UM in their usual forms are not designed to handle negative values, based on how these values are managed (taken to be zero here), the values of α will be biased further due to the effect of a lower threshold.

Though these biases have been identified, as of now, no solutions are available to account for them, and more methodological developments are required for this solution. The same is required to understand the combined effect of both biases. It is also worth mentioning that such an upper threshold is very likely to affect other statistical analyses relying on scale invariance as well. However, this is beyond the scope of the current paper and would require separate investigations. Due to the presence of the above-mentioned biases in P_t, the

actual wind power available at the turbine hub for extraction (P_a = f(v, ρ); Eq. 1) was used as a proxy to understand the small-scale variability in the follow-up UM analysis. Since the presence of thresholds in the data – imposed by limitations of operations as well as measurement – exists in many geophysical situations, understanding them is important for retrieving the actual characteristic parameters and modelling them. In the case explored here, since the characterization of power production already suffers from various influences that are poorly understood and accounted for, understanding the biases in the data treatment will help avoid more uncertainties.

## Appendix A: Diversity of multifractal formalisms and theoretical choices

In response to questions from referees and therefore potential questions from readers, we felt it necessary to discuss in greater depth the methodological choices we have made. The purpose of this Appendix is therefore to provide, in a fairly autonomous way to ease reading, a better overview of the formalisms, highlighting their strengths and weaknesses, their common features, and their diversity. This Appendix is also valid for Part 2 of this publication.

### A1 Deterministic and stochastic multifractals as well as dimension and co-dimension formalisms

In this paper and its companion paper (Jose et al., 2024), we use a stochastic multifractal framework (Schertzer and Lovejoy, 1984, 1983, 1989, 1992), for the fundamental reason that it is much more general than a deterministic multifractal framework (Parisì and Frisch, 1985; Halsey et al., 1986). This large difference is fundamentally due to the number of samples required to get reliable information. While a unique sample is sufficient for a purely deterministic process, the determination of extremes of a stochastic process may require a very large number of samples. This difference is also illustrated by the fact that a stochastic event has a finite occurrence frequency, while it may occur on an infinite number of samples. This results from the fact that a probability frequency can be understood as the (finite) limit of the ratio of two diverging numerations: “favourable cases” and “all cases”. If both numerations scale with a dimension exponent, then the frequency scales with the difference of these dimensions, which is usually called co-dimension. This can be written as follows for a process ε_λ at resolution λ = L/ℓ (outer scale L, observation scale ℓ) when assessing its probability of diverging more quickly than λ^γ, i.e. with a singularity γ:

$$ \begin{aligned} \mathrm{Pr}(\varepsilon_{\lambda} \geq \lambda^{\gamma}) &\approx \frac{N_{\lambda}(\varepsilon_{\lambda} \geq \lambda^{\gamma})}{N_{\lambda}} \approx \frac{\lambda^{D(\gamma)}}{\lambda^{D}} = \lambda^{D(\gamma)-D} \\ &= \lambda^{-c(\gamma)}, \end{aligned} \quad (A1) $$
<figure><img src="image_17.png" /><figcaption>Figure 8. UM analysis of data from 1 January to 1 April 2021 as ensembles of sample size 128 (32 min): (a) $P_a$ (as the data are), (b) $P_t$ (as the data are), (c) $P_a$ (columns with the threshold removed on the basis of the $P_t$ data below), and (d) $P_t$ with the columns of the upper threshold and zero data removed.</figcaption></figure>

where the co-dimension $c(\gamma)$ and the dimension $D(\gamma)$ satisfy

$$D = c(\gamma) + D(\gamma), \qquad (A2)$$

where $D$ is the embedding dimension of the process and thus generalizes the definition of the co-dimension $C(A)$ of a $D(A)$-dimensional sub-space $A$ in a $D$-dimensional vector space:

$$D = D(A) + C(A). \qquad (A3)$$

Equation (A1) provides the first insight into the fundamental fact that scaling exponents of probabilities are co-dimensions, while those of numerations are dimensions.

It seems paradoxical that the dimension or deterministic multifractal formalism was introduced to explain the observed non-linearity of statistical scaling exponents, specifically that of the velocity structure functions (Anselmet et al., 1984), which are the statistical moments of the velocity increments. This was done with the help of strong assumptions:

the singularities of the velocity increments, defined as local Holder exponents, were assumed to be geometrically and rather deterministically distributed over embedded fractals. In many respects, the $f(\alpha)$ formalism (Halsey et al., 1986), which deals with multifractal strange attractors, further emphasized this implicit non-random and geometric framework. A dimension formalism such as the $f(\alpha)$ formalism is formally related to the co-dimensions $(\gamma, c(\gamma))$ according to

$$\alpha_D = D - \gamma; f_D(\alpha_D) = D - c(\gamma). \qquad (A4)$$

The sub-index $D$ is introduced to $\alpha$ and $f$ for two reasons:

- Both $\alpha$ and $f$ depend on the embedding dimension $D$, e.g. by taking cuts of dimensions smaller than $D$, while $(\gamma, c(\gamma))$ do not depend on it.

- Another “historical” $\alpha$ has a quite different meaning, as briefly mentioned below.
The same dependence on $D$ occurs for the scaling exponent $\tau(q)$ of the partition function (Hentschel and Procaccia, 1983; Jiang et al., 2019), which is strongly related to that of the statistical moment scaling exponent $K(q)$,

$$<\varepsilon_{\lambda}>\approx\lambda^{K(q)}, \qquad (A5)$$

as follows:

$$\tau_D(q) = (q-1)D - K(q). \qquad (A6)$$

## A2 Partial equivalence between dimension and co-dimension formalisms

In fact, we implicitly used the partial equivalence between both formalisms for the introduction of the co-dimension (Eq. A1). Before insisting on its partiality, let us stress that it is merely defined by Eq. (A4), a very broad but rather straightforward generalization of Eq. (A3). Equation (A4) thus defines the framework transformation of going from one formalism to another.

Unfortunately, this equivalence is only partial because the dimension framework is much more limited than the co-dimension one. A major difference is that the numeration dimension $D(\gamma)$ is bounded below and above ($0 \le D(\gamma) \le D$), while the co-dimension $c(\gamma)$) being statistically defined as the opposite of the scaling exponent of the probability exponent generally has no upper bound ($0 \le c(\gamma)$). For instance, $c(\gamma) = \infty$ merely corresponds to a singularity $\gamma$ that almost never happens, including at a finite resolution. According to Eq. (A4), $c(\gamma) > D$ would correspond to a negative dimension $D(\gamma)$.

## A3 Universality and stochasticity

One important feature of the stochastic framework is that it provides universal behaviour, i.e. processes that are attractive, stable, and determined by only a limited number of parameters. We briefly recall that this is the case for the UM (Schertzer and Lovejoy, 1987, 1997) that satisfy a broad generalization of the central limit theorem: their generators are attractive and stable through re-normalized summations. We recall that their statistics are defined by the three parameters which follow and which are physically meaningful:

- For the scaling exponent $H^1$ of the mean field, when $H=0$, the mean field is strictly scale-invariant and the field is said to be conservative. $H \neq 0$ often results from a fractional integration of this order of a conservative field.

- The co-dimension $C_1 \ge 0$ of the mean field measures the mean intermittency, i.e. how the mean fluctuations

are increasingly concentrated scale by scale. When $C_1=0$, there is no intermittency and the field is statistically homogeneous.

- The multifractal index $0 \le \alpha \le 2$ measures the variability of the intermittency when departing from the mean field. When $\alpha=0$, the field is uni-fractal or mono-fractal: $\alpha=2$ corresponds to a maximal intermittency and to the so-called log-normal model. $\alpha$ is also the Levy stability index (Lévy, 1937) of the cascade generator.

The corresponding UM scaling moment function is therefore

$$K(q) = -qH + K_c(q); \quad K_c(q) = C_1 \frac{q^{\alpha} - q}{\alpha - 1}, \qquad (A7)$$

where $K_c(q)$ denotes the scaling moment function of a conservative UM field. We believe that the case studies of the textual body confirm that Eq. (A7) allows a much richer data analysis than the frequently used deterministic indicators, such as $\Delta\alpha_D = \alpha_{D,max} - \alpha_{D,min}$.

## A4 Hurst exponent and its multifractal generalizations

Scaling time series analysis has been strongly focused on the estimation of the historical Hurst exponent $H$ (Hurst, 1951), in particular with respect to its critical value $H=0.5$ that is supposed to discriminate long-range dependency and persistence ($H > 1/2$) from short-range dependency and anti-persistence ($H < 1/2$). The multifractal ideology has ruined the dogma of its uniqueness and justified divergent estimates where not accidental, but it did result from a given physics, i.e. that of intermittency. Among the many ways of defining a generalized Hurst exponent (GHE, Gómez-Aguila and Sánchez-Granero, 2021), it is straightforward to consider the scaling exponent $H(q)$ of the $q$th root of the (absolute) $q$th-order statistical moment of the field $^2$:

$$H(q) = -K(q)/q. \qquad (A8)$$

This definition is very generic and is an effective measure of the evolution of the $q$th-order statistical moment, if any, with respect to its order $q$. In particular, the uniqueness of $H(H(q) \equiv H)$ is recovered for fractional integrations of homogeneous processes. By contrast, fractional integrations of multifractal processes yield a non-constant part, e.g. for a fractional integration of a conservative field with the scaling moment function $K_c(q)$:

$$H(q) = H - K_c(q)/q. \qquad (A9)$$

In the generic case of universal multifractals, $K_c(q)$ only depends on the universal parameters $C_1$ and $\alpha$ (Eq. A7). Because $H(2)$ is often considered the historical Hurst exponent

<footnote n="1">¹As discussed below, this is related to the historical Hurst exponent, although it is rarely identical to it; see Eq. (A9) and the associated comments.</footnote>

<footnote n="2">²The minus sign that appears in this relation is only due to the left-hand side being a scaling exponent with respect to scales, whereas the right-hand side is with respect to the resolution.</footnote>
(Kantelhardt, 2002), it is only equal to $H$ for homogeneous processes according to Eq. (A9).

The main drawback of the generalized Hurst exponent $H(q)$ is that it gives access to the statistics of the cascade generator less directly than with the scaling moment function $K(q)$. This may explain why many GHE studies have limited outputs due to a lack of theoretical guidance, e.g. by only providing raw statistics of $H(q)$ such as its minimum and maximum over a given range of $q$.

## A5 Detrending frameworks

Multi-fractal detrended fluctuation analysis (MFDFA; Kantelhardt, 2002) is a popular scaling analysis technique that explicitly uses the concept of the generalized Hurst exponent. However, this is not done directly on the field of interest but rather on the standard deviations of the residues of polynomial regressions on the running sum of the fluctuations of the original time series. It thus corresponds to a multifractal extension of the (fractal) DFA; Peng et al., 1994), as generally presented.

Let us provide some details about this. Let $Y(i)$ be the cumulative fluctuation of the original time series $x(k)$ of the mean value $<x>$:

$$Y(i) = \sum_{k=1}^{i} [x(k) - <x>]. \qquad (A10)$$

The series is split into $N_s$ non-overlapping sub-series of finite size $s$, and a detrending polynomial $y_v(i)$ with constant order $m$ is fitted into each sub-series by least squares. This yields the following root mean square variation over the $v^{th}$ sub-series:

$$F(v, s) = \left[ \frac{1}{s} \sum_{i=1}^{s} \{Y[(v-1)s + i] - y_v(i)\}^2 \right]^{1/2} . \qquad (A11)$$

Averaging over the $N_s$ sub-series of size $s$ yields the total variation for the DFA analysis:

$$F(s) = \left[ \frac{1}{N_s} \sum_{v=1}^{N_s} F^2(v, s) \right]^{1/2} \qquad (A12)$$

The generalization to MFDFA is straightforwardly obtained by introducing the statistical order $q$ instead of 2:

$$F_q(s) = \left[ \frac{1}{N_s} \sum_{n=1}^{N_s} [F^2(v, s)]^{q/2} \right]^{1/q} . \qquad (A13)$$

The estimate of $H(q)$ is obtained using the logarithmic slope of $F_q(s)$:

$$F_q(s) \approx s^{H(q)} \Leftrightarrow H(q) \approx \frac{\ln F_q(s)}{\ln s}. \qquad (A14)$$

It is similar but not identical to the scaling of the trace moment of the original field.

## A6 Deferring to future work

A priori, $H(q)$ is not unique, since it may depend like $F_q(s)$ on the order $m$ of the detrending polynomials $y_v(i)$, and there is no obvious theoretical guidance on how to choose this order. In addition, the cumulative fluctuation obviously increases the order of integration $H$ by a unit (Eq. A10). Conversely, the obtained estimates (Eq. A14) must be reduced by the same amount to indirectly estimate $H(q)$ for the original series. The most serious theoretical drawback is the linear decomposition into local polynomial trends and stochastic fluctuations. Moreover, as the former is maximized by least squares, the importance of the fluctuations is minimized even though they are initially at the centre of the analysis. Because of these many issues, we defer the MFDA and variants of our data to future work.

**Code availability.** Python scripts used in this paper to extract data from the database are available along with 3 months of meteorological data in the public repository at https://doi.org/10.5281/zenodo.5801900 (Gires et al., 2021). A full description can be found in Gires et al. (2022). Python scripts used for this paper to carry out multifractal and joint multifractal analysis are available in the public repository at https://doi.org/10.5281/zenodo.3707904 (Gires et al., 2020).

**Data availability.** Three months of the meteorological data used in this paper can be accessed in the public repository at https://doi.org/10.5281/zenodo.3707904 (Gires et al., 2020). A full description can be found in the data paper (Gires et al., 2022). Longer time series can be accessed by contacting the authors.

**Author contributions.** AG, JJ, IT, and DS designed the initial content of the study. JJ performed the data analysis and added the theoretical studies and numerical simulations following it. JJ wrote the main content of the paper under the supervision of AG, IT, DS, ES, and YR. All the authors contributed to the review of the paper.

**Competing interests.** At least one of the (co-)authors is a member of the editorial board of *Nonlinear Processes in Geophysics*. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
**Acknowledgements.** The authors acknowledge the partial financial support from the Chair of Hydrology for Resilient Cities (endowed by Veolia) of the École nationale des ponts et chaussées.

**Financial support.** This research has been supported by the ANR JCJC RW-Turb project (grant no. ANR-19-CE05-0022-01).

**Review statement.** This paper was edited by Luciano Telesca and reviewed by Felipe Segundo Abril Bermudez and one anonymous referee.

## References

Anselmet, F., Gagne, Y., Hopfinger, E., and Antonia, R.: High-order velocity structure functions in turbulent shear flows, J. Fluid Mech., 140, 63–89, 1984.

Beiter, P., Cooperman, A., Lantz, E., Stehly, T., Shields, M., Wiser, R., Telsnig, T., Kitzing, L., Berkhout, V., and Kikuchi, Y.: Wind power costs driven by innovation and experience with further reductions on the horizon, WIREs Energy Environ., 10, e398, https://doi.org/10.1002/wene.398, 2021.

Calif, R. and Schmitt, F. G.: Multiscaling and joint multiscaling description of the atmospheric wind speed and the aggregate power output from a wind farm, Nonlin. Processes Geophys., 21, 379–392, https://doi.org/10.5194/npg-21-379-2014, 2014.

Chambers, J. M., Mallows, C. L., and Stuck, B.: A method for simulating stable random variables, J. Am. Stat. Assoc., 71, 340–344, 1976.

Fitton, G., Tchiguirinskaia, I., Schertzer, D., and Lovejoy, S.: Scaling Of Turbulence In The Atmospheric Surface-Layer: Which Anisotropy?, J. Phys. Conf. Ser., 318, 072008, https://doi.org/10.1088/1742-6596/318/7/072008, 2011.

Fitton, G., Tchiguirinskaia, I., Schertzer, D., and Lovejoy, S.: Torque Fluctuations In The Framework Of A Multifractal 23/9-Dimensional Turbulence Model, J. Phys. Conf. Ser., 555, 012038, https://doi.org/10.1088/1742-6596/555/1/012038, 2014.

Gires, A., Tchiguirinskaia, I., Schertzer, D., and Lovejoy, S.: Influence of the zero-rainfall on the assessment of the multifractal parameters, Adv. Water Resour., 45, 13–25, https://doi.org/10.1016/j.advwatres.2012.03.026, 2012.

Gires, A., Tchiguirinskaia, I., Schertzer, D., Schellart, A., Berne, A., and Lovejoy, S.: Influence of small scale rainfall variability on standard comparison tools between radar and rain gauge data, Atmos. Res., 138, 125–138, https://doi.org/10.1016/j.atmosres.2013.11.008, 2014.

Gires, A., Tchiguirinskaia, I., and Schertzer, D.: Data for “Approximate multifractal correlation and products of universal multifractal fields, with application to rainfall data” by Auguste Gires, Ioulia Tchiguirinskaia, and Daniel Schertzer, NPG 2020, Zenodo [data set], https://doi.org/10.5281/zenodo.3707904, 2020.

Gires, A., Jose, J., Tchiguirinskaia, I., and Schertzer, D.: Data for “Three months of combined high resolution rainfall and wind data collected on a wind farm”, Zenodo [data set], https://doi.org/10.5281/zenodo.5801900, 2021.

Gires, A., Jose, J., Tchiguirinskaia, I., and Schertzer, D.: Combined high-resolution rainfall and wind data collected for 3 months

on a wind farm 110 km southeast of Paris (France), Earth Syst. Sci. Data, 14, 3807–3819, https://doi.org/10.5194/essd-14-3807-2022, 2022.

Gómez-Aguila, A. and Sánchez-Granero, M.: A theoretical framework for the TTA algorithm, Physica A, 582, 126288, https://doi.org/10.1016/j.physa.2021.126288, 2021.

Guezuraga, B., Zauner, R., and Pölz, W.: Life cycle assessment of two different 2 MW class wind turbines, Renew. Energ., 37, 37–44, https://doi.org/10.1016/j.renene.2011.05.008, 2012.

Halsey, T. C., Jensen, M. H., Kadanoff, L. P., Procaccia, I., and Shraiman, B. I.: Fractal measures and their singularities: The characterization of strange sets, Phys. Rev. A, 33, 1141, https://doi.org/10.1103/PhysRevA.33.1141, 1986.

Hentschel, H. G. E. and Procaccia, I.: The infinite number of generalized dimensions of fractals and strange attractors, Physica D, 8, 435–444, 1983.

Hubert, P. and Carbonnel, J.: Dimensions fractales de l'occurrence de pluie en climat soudano-sahélien, Hydrologie Continentale, 4, 3–10, https://www/documentation.ird.fr/hor/fdi:27277 (last access: 26 November 2024), 1989.

Hurst, H. E.: Long-term storage capacity of reservoirs, T. Am. Soc. Civ. Eng., 116, 770–799, 1951.

Jiang, Z.-Q., Xie, W.-J., Zhou, W.-X., and Sornette, D.: Multifractal analysis of financial markets: A review, Rep. Prog. Phys., 82, 125901, https://doi.org/10.1088/1361-6633/ab42fb, 2019.

Johnson, K. E.: Adaptive Torque Control of Variable Speed Wind Turbines, NREL/TP-500-36265, 107 pp., https://www.nrel.gov/docs/fy04osti/36265.pdf (last access: 26 November 2024), 2004.

Jose, J., Gires, A., Tchiguirinskaia, I., and Schertzer, D.: Influence of lower threshold on empirical data in estimation of multifractal parameters (using atmospheric extinction coefficient as the field of study), in: AGU Fall Meeting Abstracts, vol. 2021, NG45C-0591, https://enpc.hal.science/hal-04580214/ (last access: 26 November 2024), 2021.

Jose, J., Gires, A., Schnorenberger, E., Roustan, Y., Schertzer, D., and Tchiguirinskaia, I.: Multifractal analysis of wind turbine power and rainfall from an operational wind farm – Part 2: Joint analysis of available wind power and rain intensity, Nonlin. Processes Geophys., 31, 603–624, https://doi.org/10.5194/npg-31-603-2024, 2024.

Jung, C. and Schindler, D.: The role of air density in wind energy assessment – A case study from Germany, Energy, 171, 385–392, https://doi.org/10.1016/j.energy.2019.01.041, 2019.

Kantelhardt, J.: Multifractal deterended fluctuation analysis of non-stationary time series, Physica, 316, 81–91, 2002.

Lavallée, D., Lovejoy, S., Schertzer, D., and Ladyo, P.: Nonlinear variability and landscape topography: analysis and simulation, Fractals in Geography, edited by: De Cola, L. and Lam, N. 158–192, PTR, Prentice Hall, https://www.physics.mcgill.ca/~gang/eprints/eprintLovejoy/neweprint/topoall.pdf (last access: 26 November 2024), 1993.

Lévy, P.: Théorie de l'addition des variables aléatoires, in: Collection des monographies des probabilités, ac37003071, Gauthier-Villars, https://books.google.fr/books?id=1d6iuAEACAAJ (last access: 26 November 2024), 1937.

Li, J., Li, S., and Wu, F.: Research on carbon emission reduction benefit of wind power project based on life cycle assessment theory, Renew. Energ., 155, 456–468, https://doi.org/10.1016/j.renene.2020.03.133, 2020.
Lovejoy, S. and Schertzer, D.: Scale, Scaling and Multifractals in Geophysics: Twenty Years on, in: Nonlinear Dynamics in Geosciences, Springer, New York, NY, https://doi.org/10.1007/978-0-387-34918-3_18, 2007.

Lovejoy, S. and Schertzer, D.: The Weather and Climate: Emergent Laws and Multifractal Cascades, Cambridge University Press, https://doi.org/10.1017/CBO9781139093811, 2013.

Mandelbrot, B. B.: The fractal geometry of nature, vol. 1, WH freeman New York, ISBN 0-7167-1186-9, 1982.

Parisi, G. and Frisch, U.: On the singularity structure of fully developed turbulence, in: Turbulence and Predictability in Geophysical Fluid Dynamics and Climate Dynamics, North-Holland Publisher, Amsterdam, The Netherlands, 84–87, 1985.

Peng, C.-K., Buldyrev, S. V., Havlin, S., Simons, M., Stanley, H. E., and Goldberger, A. L.: Mosaic organization of DNA nucleotides, Phys. Rev. E, 49, 1685, https://doi.org/10.1103/PhysRevE.49.1685, 1994.

Picard, A., Davis, R. S., Gläser, M., and Fujii, K.: Revised formula for the density of moist air (CIPM-2007), Metrologia, 45, 149, https://doi.org/10.1088/0026-1394/45/2/004, 2008.

Schertzer, D. and Lovejoy, S.: Elliptical turbulence in the atmosphere, in: Symposium on Turbulent Shear Flows, Fourth International Symposium on Turbulent Shear Flows, University of Karlsruhe, Karlsruhe, FRG, https://link.springer.com/book/10.1007/978-3-642-69996-2 (last access: 26 November 2024), 12–14 September 1983.

Schertzer, D. and Lovejoy, S.: On the dimension of atmospheric motions, in: Turbulence and Chaotic Phenomena in Fluids, edited by: Tatsumi, T., Amsterdam, North Holland, 505–12, ISBN 0-444-87-594-8, 1984.

Schertzer, D. and Lovejoy, S.: The dimension and intermittency of atmospheric dynamics, in: Turbulent Shear Flows 4: Selected Papers from the Fourth International Symposium on Turbulent Shear Flows, University of Karlsruhe, Karlsruhe, FRG, 12–14 September 1983, Springer, 7–33, https://www.physics.mcgill.ca/~gang/eprints/eprintLovejoy/neweprint/TSFall.pdf (last access: 26 November 2024), 1985.

Schertzer, D. and Lovejoy, S.: Physical modeling and analysis of rain and clouds by anisotropic scaling multiplicative processes, J. Geophys. Res.-Atmos., 92, 9693–9714, https://doi.org/10.1029/JD092iD08p09693, 1987.

Schertzer, D. and Lovejoy, S.: Multifractal simulations and analysis of clouds by multiplicative processes, Atmos. Res., 21, 337–361, https://doi.org/10.1016/0169-8095(88)90035-X, 1988.

Schertzer, D. and Lovejoy, S.: Nonlinear Variability in Geophysics: Multifractal Simulations and Analysis, Springer US, Boston, MA, 49–79, ISBN 978-1-4899-3499-4, https://doi.org/10.1007/978-1-4899-3499-4_3, 1989.

Schertzer, D. and Lovejoy, S.: Hard and soft multifractal processes, Physica A, 185, 187–194, https://doi.org/10.1016/0378-4371(92)90455-Y, 1992.

Schertzer, D. and Lovejoy, S.: EGS Richardson AGU Chapman NVAG3 Conference: Nonlinear Variability in Geophysics: scaling and multifractal processes, Nonlin. Processes Geophys., 1, 77–79, https://doi.org/10.5194/npg-1-77-1994, 1994.

Schertzer, D. and Lovejoy, S.: Universal Multifractals Do Exist!: Comments on “A Statistical Analysis of Mesoscale Rainfall as a Random Cascade”, J. Appl. Meteorol., 36, 1296–1303, https://doi.org/10.1175/1520-0450(1997)036<1296:UMDECO>2.0.CO;2, 1997.

Schertzer, D. and Tchiguirinskaia, I.: A Century of Turbulent Cascades and the Emergence of Multifractal Operators, Earth and Space Science, 7, e2019EA000608, https://doi.org/10.1029/2019EA000608, 2020.

Tessier, Y., Lovejoy, S., and Schertzer, D.: Universal Multifractals: Theory and Observations for Rain and Clouds, J. Appl. Meteorol. Clim., 32, 223–250, https://doi.org/10.1175/1520-0450(1993)032<0223:UMTAOF>2.0.CO;2, 1993.

Ulazia, A., Gonzalez-Rojí, S. J., Ibarra-Berastegi, G., Carrenot-Madinabeitia, S., Sáenz, J., and Nafarrate, A.: Seasonal Air Density Variations over the East of Scotland and The Consequences for Offshore Wind Energy, in: Proceedings of the 2018 7th International Conference on Renewable Energy Research and Applications (ICRERA), Paris, France, 14–17 October 2018, 261–265, https://ieeexplore.ieee.org/document/8566716 (last access: 26 November 2026), 2018.

van Kuik, G. A. M., Peinke, J., Nijssen, R., Lekou, D., Mann, J., Sørensen, J. N., Ferreira, C., van Wingerden, J. W., Schlipf, D., Gebraad, P., Polinder, H., Abrahamsen, A., van Bussel, G. J. W., Sørensen, J. D., Tavner, P., Bottasso, C. L., Muskulus, M., Matha, D., Lindeboom, H. J., Degraer, S., Kramer, O., Lehnhoff, S., Sonnenschein, M., Sørensen, P. E., Künneke, R. W., Morthorst, P. E., and Skytte, K.: Long-term research challenges in wind energy – a research agenda by the European Academy of Wind Energy, Wind Energ. Sci., 1, 1–39, https://doi.org/10.5194/wes-1-1-2016, 2016.

Vestas Wind Systems A/S, V.: V90-2.0 MW™ IEC IIA/IEC S turbines, https://www.vestas.com/en/pages/backup-2-mw-platform/V90-2-0-MW (last access: 26 November 2024), 2023.

Wiser, R., Yang, Z., Hand, M., Hohmeyer, O., Infield, D., Jensen, P., Nikolaev, V., O'Malley, M., and Zervos, G.: Wind Energy, in: IPCC Special Report on Renewable Energy Sources and Climate Change Mitigation, edited by: Edenhofer, O., Pichs-Madruga, R., Sokona, Y., Seyboth, K., Matschoss, P., Kadner, S., Zwickel, T., Eickemeier, P., Hansen, G., Schlömer, S., and von Stechow, C., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA https://www.ipcc.ch/site/assets/uploads/2018/03/Chapter-7-Wind-Energy-1.pdf (last access: 16 November 2024), 2011.