# Multifractal analysis of oceanic chlorophyll maps remotely sensed from space

Louis de Montera, Manel Jouini, Sébastien Verrier, Sylvie Thiria, Michel Crépon

► To cite this version:

Louis de Montera, Manel Jouini, Sébastien Verrier, Sylvie Thiria, Michel Crépon. Multifractal analysis of oceanic chlorophyll maps remotely sensed from space. Ocean Science, 2011, 7 (2), pp.219-229. 10.5194/os-7-219-2011. hal-00567118

HAL Id: hal-00567118

https://hal.science/hal-00567118v1

Submitted on 10 Dec 2015

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 CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Multifractal analysis of oceanic chlorophyll maps remotely sensed from space

L. de Montera<fnref n="1" />, M. Jouini<fnref n="1" />, S. Verrier<fnref n="2" />, S. Thiria<fnref n="1" />, and M. Crepon<fnref n="1" />

Received: 14 October 2010 – Published in Ocean Sci. Discuss.: 21 January 2011

Revised: 24 March 2011 – Accepted: 25 March 2011 – Published: 28 March 2011

**Abstract.** Phytoplankton patchiness has been investigated with multifractal analysis techniques. We analyzed oceanic chlorophyll maps, measured by the SeaWiFS orbiting sensor, which are considered to be good proxies for phytoplankton. The study area is the Senegalo-Mauritanian upwelling region, because it has a low cloud cover and high chlorophyll concentrations. Multifractal properties are observed, from the sub-mesoscale up to the mesoscale, and are found to be consistent with the Corssin-Obukhov scale law of passive scalars. This result indicates that, in this specific region and within this scale range, turbulent mixing would be the dominant effect leading to the observed variability of phytoplankton fields. Finally, it is shown that multifractal patchiness can be responsible for significant biases in the nonlinear source and sink terms involved in biogeochemical numerical models.

## 1 Introduction

It is sometimes argued that turbulent mixing leads to homogeneous fields. However, Kolmogorov (1941), Obukhov (1949) and Corssin (1951) have shown that, on the contrary, turbulent mixing generates highly irregular structures that are heterogeneous at all scales. Their work was based on the hypothesis of scale invariance, which means, in simple words, that eddies can be expected to occur in a similar manner at all scales. In the case where the physical quantity is the concentration of a passive tracer, these authors demonstrated that its variability exhibits fractal properties which can be described

statistically using scale laws. This result is often referred to as the theory of passive scalars.

Since the phytoplankton's ability to swim is very limited, its displacements are mainly due to the velocity of the fluid in which it evolves. Phytoplankton patchiness is thus strongly related to turbulence. This consequence has led numerous authors to study the scale invariance properties of phytoplankton patches, and to confront experimental data with phenomenological models derived from, or inspired by, the theory of passive scalars. Early studies remained confined to second-order moments, such as the slope of the power spectrum (see, e.g., Platt, 1972), whereas more recent research takes into account the intermittent transfer of conservative quantities in scale space, such as energy and scalar variance, which give rise to multifractal statistics through cascade processes (Seuront et al., 1996a, b, 1999; Seuront and Schmitt, 2004, 2005a, b; Lovejoy et al., 2001a, b; Pottier et al., 2008).

At smaller scales, most of these studies found empirical proof for a passive scalar regime of phytoplankton patchiness, corresponding to the well-known "-5/3" power spectrum slope of homogenous and isotropic turbulence. This purely turbulent regime appears to be limited to spatial scales smaller than a particular scale of the order of 100 m, called the "planktoscale" by Lovejoy et al. (2001b). In fact, although phytoplankton can reasonably be described as passively advected, in the sense that its retroaction on the turbulent flow is negligible, it cannot be considered to be totally passive, since it is biologically active. One important biological process is zooplankton grazing, and the "planktoscale" is currently interpreted as corresponding to the scale at which changes take place in the grazing regime. This modification of the grazing regime appears to be related to the zooplankton's ability to swim. Contrary to phytoplankton, zooplankton is able to swim, although its speed remains limited. Therefore, there exists a scale above which its displacements

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

Correspondence to: S. Thiria
(sylvie.thiria@locean-ipsl.upmc.fr)

<footnote n="1">¹Laboratoire d'Océanographie et du Climat, Expérimentation et Approches Numériques, UMR 7159, Paris, France</footnote>

<footnote n="2">²Laboratoire Atmosphères, Milieux, Observations Spatiales, UMR 8190, Guyancourt, France</footnote>
are dominated by turbulent mixing. This hypothesis is supported by the fact that the zooplankton's concentration power spectrum whitens at scales smaller than the "planktoscale" (Currie and Roff, 2006).

At scales larger than the "planktoscale", the phytoplankton patchiness description is more confused. Some authors found a power spectrum slope steeper than -5/3 (around -2, Currie and Roff, 2006; Seuront et al., 1999), interpreting this result as a transition from Eulerian to Lagrangian statistics, due to the inertia of the boat carrying the instruments (cf. Seuront et al., 1996b). On the other hand, the analysis of remotely sensed phytoplankton fields from aircraft led to a smoother slope (around -1.2, Lovejoy et al., 2001b). From a theoretical point of view, the situation is even less clear: some studies reached the conclusion that growth and trophic interactions should decrease the slope of the power spectrum (Denman and Platt, 1976; Fasham, 1978), whereas some others predict that it should increase (Steele and Henderson, 1979) or that the power spectrum has no specific regime (Horwood, 1978).

In this context, to the best of our knowledge, the large volumes of data collected by remote sensing from space, over a period of more than two decades, have almost not been exploited in order to improve scientific understanding of the multi-scaling properties of phytoplankton fields (with the noticeable exception of Nieves et al., 2007). The aim of the present study is to analyze oceanic chlorophyll maps obtained through the use of this type of sensor. The first part of the paper briefly recalls the passive scalar theory and the notion of multifractal intermittency. The second and third parts describe both the cascade model and the analysis technique. The fourth part presents the dataset and the pre-treatments. The fifth and sixth parts are dedicated to the results and their interpretation. Finally, the last part of the paper provides an example of the importance of multifractal patchiness in oceanic tracers, by assessing the biases it produces in bio-geochemical numerical models.

## 2 Theoretical background: turbulence and multifractals

Richarsdon (1922) described turbulence as a cascade process that transfers kinetic energy from large scales to small scales by a hierarchy of imbricated eddies. The hypothesis of scale invariance relies on phenomenology and the invariance of the Navier-Stokes equations under dilatation or contraction of the reference system (see Appendix A of Schertzer and Lovejoy, 1987). On the basis of this hypothesis, and the conservation of energy in the inertial range, Kolmogorov (1941) used dimensionality and some general assumptions, such as homogeneity and isotropy, to derive his famous statistical scale law:

$$ \Delta v_l \sim \epsilon^{1/3} l^{1/3} \quad (1) $$

In this equation, $\Delta v_l \equiv \langle |v(x+l) - v(x)| \rangle$ represents the mean shear of (longitudinal) velocity between two points separated by a distance $l$, and $\epsilon$ represents the mean density of the energy flux, which is equal to the rate of energy dissipation per unit mass. A similar scale law has been derived for the concentration $C$ of a passive tracer (Obukhov, 1949; Corrsin, 1951):

$$ \Delta C_l \sim \varphi^{1/3} l^{1/3} \quad (2) $$

where $\varphi \equiv \chi^{3/2} \epsilon^{-1/2}$ represents the non-linear coupling between velocity and concentration, with $\chi \equiv - \frac{\partial(\Delta C_l)^2}{\partial t}$ the mean density of the concentration variance flux (for a review concerning these early turbulent models, see, e.g., Panchev, 1971).

Landau and Lifshitz (1944) pointed out that these fluxes have no reason to be homogeneous: although they are on average conserved during the cascade process, their transfer is a priori intermittent. This remark has led the transfer process to be described by stochastic multiplicative cascades (Novikov and Stewart, 1964; Yaglom, 1966). A multiplicative cascade can be constructed by iterating the following simple procedure: (i) distribute a quantity uniformly over an interval, (ii) divide this interval into sub-intervals, (iii) multiply these by a random variable in order to obtain the new quantity for each sub-interval, (iv) repeat steps (ii) and (iii) until the smallest scale of the cascade is reached. The important point here is that the distribution of the random variables, referred to as the multiplicative weights in the following, does not depend on the level of iteration of the construction algorithm. Thus, because the latter is not dependent on scale, the resulting mathematical object has fractal and even multifractal properties. It turns out that these properties can be described by the scaling of its statistical moments, of fractional order $q$ (for more details, see Schertzer et al., 2002):

$$ \langle \epsilon_l^q \rangle \sim \left( \frac{L}{l} \right)^{K_\epsilon(q)} \quad (3) $$

$$ \langle \chi_l^q \rangle \sim \left( \frac{L}{l} \right)^{K_\chi(q)} \quad (4) $$

where $\epsilon_l$ and $\chi_l$ are the fluxes averaged at scale $l$, $L$ is the largest scale of the cascade, and $K_\epsilon(q)$ and $K_\chi(q)$ are the so-called moment scaling functions.

In order to obtain realistic fields, the discrete cascade model described above has been generalized to continuous cascades, obtained by scale densification (Schertzer and Lovejoy, 1987). This generalization was necessary because two points separated by a given length in physical space do not always have the same distance to their closest common ancestor in the cascade (cf. Pecknold et al., 1993). An interesting property of continuous cascades is that their generators (i.e., the logarithm of the random multiplicative weights) converge towards infinitely divisible laws. However, there is still no consensus concerning the degree of convergence. Some authors proposed Poisson generators (She
and Lévêque, 1994; Dubrulle, 1994) whereas some others add an assumption of self-similar renormalization, so that the generator converges more accurately towards stable distributions (Schertzer and Lovejoy, 1987, 1997). Until this question finds a definitive answer, since the notion of scale invariance is at the root of the theory, the latter assumption seems plausible, and the decision was made to use stable laws to describe the generator of the process which correspond to Gaussian distributions if the variance is finite or Levy distributions if the variance is infinite. In this type of case, the moment scaling functions take the simple form (Schertzer and Lovejoy, 1987):

$$K_{\varepsilon}(q) = \frac{C_{1\varepsilon}}{\alpha_{\varepsilon}-1} (q^{\alpha_{\varepsilon}} - q) \qquad (5)$$

$$K_{\chi}(q) = \frac{C_{1\chi}}{\alpha_{\chi}-1} (q^{\alpha_{\chi}} - q) \qquad (6)$$

where $\alpha_{\varepsilon}$ and $\alpha_{\chi}$ are multifractality parameters varying between 0 and 2, and $C_{1\varepsilon}$ and $C_{1\chi}$ are intermittency parameters varying between 0 and the dimension of the embedding space, which here is equal to 2.

### 3 The FIF model

Concerning the chlorophyll concentration, these laws cannot be directly applied, the main reason being that biological activities may produce deviations from a purely passive scalar behaviour. Nevertheless, we expect that the variability of chlorophyll maps still presents some multifractal properties, and that it would be possible to use a cascade model similar to that presented above. We thus looked for a phenomenological model having the same form as Eq. (1), but in which the parameters are not known, i.e.:

$$\Delta\text{Chl}_l \approx \langle \zeta^a \rangle^l H \qquad (7)$$

where Chl is the chlorophyll concentration and $\Delta\text{Chl}_l \equiv \langle |\text{Chl}(x+l) - \text{Chl}(x)| \rangle$. $\zeta$ is a conserved flux and $a$ and $H$ are adjustable parameters. As described above, the conserved flux has to verify the basic multifractal relation:

$$\langle \zeta_l^q \rangle \approx \left( \frac{L}{l} \right)^{K_\zeta(q)} \qquad (8)$$

and it is assumed that it converges towards a log-stable law:

$$K_{\zeta}(q) = \frac{C_{1\zeta}}{\alpha_{\zeta}-1} (q^{\alpha_{\zeta}} - q). \qquad (9)$$

This model is described by four parameters: $a$, $H$, $\alpha_{\zeta}$ and $C_{1\zeta}$. However, it is possible to reduce the number of parameters to three, since taking the $a$-th power of $\zeta$ in Eq. (7) is equivalent to a simple shift of $H$ by $K_{\zeta}(a)$, and to the multiplication of $C_{1\zeta}$ by a factor $a^{\alpha_{\zeta}}$ (Lavallée et al., 1993). The proof uses the $q$-th-order structure functions of chlorophyll maps defined by the average of the $q$-th power of

chlorophyll concentration variations, denoted by $\Delta\text{Chl}_l^q \equiv \langle |\text{Chl}(x+l) - \text{Chl}(x)|^q \rangle$. According to Eq. (7), the $q$-th-order structure function is equal to:

$$\Delta\text{Chl}_l^q \approx \langle \zeta^{aq} \rangle l^{qh}. \qquad (10)$$

By introducing Eq. (8), this equation simplifies to:

$$\Delta\text{Chl}_l^q \approx l^{qH-K_\zeta(aq)}. \qquad (11)$$

Then, the term $K_{\zeta}(aq)$ can be decomposed into conservative and non-conservative parts using the following identity, which can be straightforwardly derived from Eq. (9):

$$K(ap) = qK(a) + a^{\alpha}K(q). \qquad (12)$$

This operation yields:

$$\Delta\text{Chl}_l^q \approx l^{q(H-K(a))-a^{\alpha_{\zeta}} K_{\zeta}(q)}. \qquad (13)$$

As expected, by defining:

$$\begin{cases} K'_{\zeta}(q) \equiv a^{\alpha_{\zeta}} K_{\zeta}(q) \\ H' \equiv H - K_{\zeta}(a), \end{cases} \qquad (14)$$

one obtains the usual form of the structure function corresponding to the case where $a=1$:

$$\Delta\text{Chl}_l^q \approx l^{qH'-K'_\zeta(q)}. \qquad (15)$$

In the following, we thus use a simplified form of Eq. (7) requiring only three parameters, namely, $H$, $\alpha_{\zeta}$ and $C_{1\zeta}$:

$$\Delta\text{Chl}_l \approx \zeta^l H \qquad (16)$$

This model is called the fractionally integrated flux (FIF) and was first presented by Schertzer and Lovejoy (1987) in the framework of their study of rain fields. In this regard, it is interesting to note the similarities between marine biogeochemistry and the cycle of water in the atmosphere. Firstly, both are strongly dependent on ascending currents: these currents bring nutrients to the surface layers of the ocean and water vapour to the upper layers of the atmosphere. Then, the first phase transition to heavier particles generally occurs in thin layers in which physical conditions are appropriate: phytoplankton is produced near to the ocean's surface because it needs light to grow, whereas clouds are formed in the atmospheric layer in which water vapour condenses. The next phase transition to heavier particles occurs when phytoplankton feeds zooplankton, and when cloud droplets are incorporated into raindrops. Finally, zooplankton may die and sink (or return to a nutrient form by mineralization), whereas raindrops may fall (or return to water vapour by evaporation). Although these atmospheric and oceanic processes have very different space and time scales, with one involving biology and the other physics, the comparison is striking. The interesting point here is that, in both cases, the evolution cycle is an alternating composition of turbulent mixing and phase transition processes. This analogy leads phytoplankton patchiness to be thought of as “clouds in the sea”. Besides, as will be shown below, the multifractal parameters obtained from chlorophyll maps are close to those obtained in the case of cloud or rain fields.
## 4 Analysis technique

The first step of the analysis consists in verifying the scale law given by Eq. (9), and in estimating its exponent $H$. This is generally performed by using the first-order structure function. Since the flux $\zeta$ is assumed to be conserved in scale space, whatever the scale $l$, $\langle \zeta_l \rangle$ is constant. Therefore, Eq. (9) reduces to:

$$ \Delta\text{Chl}_l \propto l^H. \qquad (17) $$

This equation allows $H$ to be estimated using the simple expression:

$$ H = \log_l (\Delta\text{Chl}_l). \qquad (18) $$

The second step consists in quantifying the multifractal properties of the flux, and in estimating $\alpha_\zeta$ and $C_{1\zeta}$. In order to do so, it is necessary to reconstruct the cascade and therefore to retrieve the flux at the finest available scale. According to Eq. (9), this requires a fractional derivative of order $H$. However, a simple derivation of integer order provides a good numerical approximation (Lavallée et al., 1993), such as taking the norm of the gradient of the field:

$$ \zeta_{l\max} \approx \sqrt{\left(\frac{d\text{Chl}}{dx}\right)^2 + \left(\frac{d\text{Chl}}{dy}\right)^2}. \qquad (19) $$

Note that, since the rest of the analysis is based on the gradient of the field, it is crucial to work with data affected by a low level of noise. Indeed, if the noise is strong, taking the gradient of the field will result in useless, noisy fields. Therefore the finest available scale does not necessarily correspond to the measurement scale: it is usually necessary to perform initial averaging of the data at a larger scale, before computing the gradient, in order to suppress the noisiest of the finest scales. The cut-off scale at which the fields have to be averaged is called the “effective measurement scale” in the following. This scale is determined by computing the power spectrum, and then estimating the wave number above which it flattens out.

Once the flux has been obtained at the “effective measurement scale”, the stochastic multiplicative cascade can be reconstructed by averaging (or “degrading”) the flux at larger scales. The statistical moments are then computed for various orders and scales in order to test Eq. (8). If the scaling of the statistical moments is verified, $K_\zeta(q)$ can be estimated. Finally, the parameters $\alpha_\zeta$ and $C_{1\zeta}$ are obtained by determining the least squares fit to this function.

## 5 Dataset

Particular attention was paid to the selection of chlorophyll maps, because multifractal analysis is very sensitive to the quality of the data. As explained above, the analysis technique is based on the gradient of the field. Therefore, if the

signal is too noisy, the fluctuations due to turbulence or other processes will be hidden. Moreover, the estimation of higher-order statistical moments can easily be biased by the presence of a few unrealistic values in the data, such as isolated pixels having abnormally high chlorophyll concentrations.

Another difficulty is that areas below clouds or high aerosol concentrations cannot be observed, because the sensor cannot see the sea surface. As a consequence, chlorophyll maps remotely sensed from space present many “holes” of different sizes (the set defined by the locations of these missing data may be fractal itself, because cloud and aerosol distributions are also fractal, see respectively Lovejoy and Schertzer, 2006 and Lilley et al., 2004). We also noticed that the values around the periphery of these “holes” were not reliable, presumably because of uncertainties in the correction of the atmospheric effect. Therefore, it was chosen to study only maps which had no missing values. This type of data is of course difficult to find, because of the abundance of clouds and aerosols in the atmosphere, such that a compromise needed to be found between the size of the maps and the sample size.

In order to optimize this compromise, the study area was carefully chosen. The most appropriate area was found to be the Senegalo-Mauritanian upwelling region, because it normally has a very low cloud cover (although this does not remain true during the summer months, when the InterTropical Convergence Zone (ITCZ) moves north). Another reason is the presence of upwelling, which provides high chlorophyll concentrations far from the coast, due to peculiar oceanic conditions (Aristegui et al., 2004; Lathuilière et al., 2008). Therefore, the choice of this area reduced the measurement noise with respect to the coherent signal. Finally, the chosen location lies between 10° N–26° N and 14° E–26° E, which corresponds to the area between the Cape Verde islands and the coast of West Africa, between Mauritania and Guinea-Bissau (see Fig. 1).

The choice of product level also has to be carefully considered. Classical 8-day composite maps could not be used, because they include a non-uniform time averaging in the data, depending on the amount of missing data for each pixel. Level L3 products (daily global maps mapped to a uniform scale grid) could not be used either, firstly because of projection effects, and secondly because this product is actually derived from sub-sampling of the original data: only one pixel is kept for each square of 4 × 4 pixels. This reduction in the amount of data is unavoidable, because SeaWiFS is positioned in Low Earth Orbit (LEO), meaning that the time it remains within the field of view of receiver stations it too short for full datasets to be transmitted to the ground. However, full resolution chlorophyll data was transmitted for some restricted areas, including the Senegalo-Mauritanian upwelling region. This dataset is called the local unmapped level L2 product, and is suitable for use in this study. In this product, the pixel resolution is around 1 km². However, the spot size over which the data are measured varies with
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Geographic map of the Senegalo-Mauritanian upwelling region.</figcaption></figure>

elevation angle. Therefore, only the inner part of the scans was considered, in order to limit this effect. Note also that some authors (Lovejoy et al., 2001b) recommend using direct analysis of marine reflectivities (level L1 product) because the fields' heterogeneity may bias chlorophyll concentration retrieval algorithms. Indeed, since these algorithms are generally non-linear, it is not correct to extrapolate them directly to the measurement scale which is much larger than the scale of homogeneity. However, this problem should not affect our analysis because the retrieval of chlorophyll concentration was performed without any extrapolation in scale space (the retrieval algorithm is based on an empirical relation derived from the comparison between remotely sensed marine reflectivities and in-situ chlorophyll concentrations). Moreover, working directly with marine reflectivities is difficult because of its lack of physical interpretation. Actually, the only physical quantity that can be related to a theoretical scale law is the chlorophyll concentration Chl. For example, consider a non-linear relationship of the form $Chl = f(R)$, where R denotes a marine reflectivity. If f is non-linear, then $f^{-1}$ is also non-linear. There is consequently no reason for the marine reflectivity $R = f^{-1}(Chl)$ to verify a scaling of the form of Eq. (9) because non-linear transformations do not generally conserve first-order structure functions.

Finally, 100 maps of 128 × 128 pixels of 1 km², with a minimum of 99.5% of available data, were extracted from SeaWiFS data over a period of one year running from July 2003 to June 2004 (a sample chlorophyll map is shown in Fig. 2). The few missing data were interpolated automatically by computing the mean of the surroundings pixels. All selected maps were checked manually. Some maps had to be rejected because of an offset affecting some parts of the field. The origin of this offset is not known. The gradient of the maps was also checked manually, in order to detect any

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Example of a 128 km² horizontal chlorophyll map (resolution 1 km²) extracted from the SeaWiFS local L2 product.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** First-order structure function of SeaWiFS chlorophyll maps compared with a linear curve of slope equal to 0.4. The departure from the theoretical fit observed at the finest scales is attributed to measurement noise, and corresponds to flattening of the power spectra at high wave numbers (see Fig. 4).</figcaption></figure>

isolated, unrealistically high values. Each of these unrealistic pixels was corrected using the mean value of the surrounding pixels.

# 6 Results

Figure 3 shows the first-order structure function for the 100 SeaWiFS chlorophyll maps. The smaller scales (1–4 km) were not taken into account when determining the fit, because they present a deviation from the scaling observed throughout the remainder of the scale range (4–128 km). The cause of this deviation does not appear to be physical,
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Angle-integrated power spectrum of SeaWiFS chlorophyll maps compared with a linear curve of slope equal to -1.67. The power spectrum flattens out at the fifth octave, corresponding to wavelengths smaller than 4 km.</figcaption></figure>

because such a break in the scaling has never been observed in other studies (cf. Lovejoy et al., 2001b). This break was therefore associated with the scale below which the measurement noise becomes dominant, when compared with the coherent signal (here, the definition of the noise is very large: it includes not only the sensor's sensitivity, but also atmospheric corrections and retrieval algorithms errors). This hypothesis is confirmed by the power spectrum (Fig. 4), which flattens out beyond a wave number corresponding to 4 km in the physical space. Finally, over the scale range 4–128 km, the empirical first-order structure function is consistent with Eq. (17), and *H* is estimated to be around 0.4 (the numerical fit yields 0.402 with a standard deviation of the estimator equal to 0.005).

Since the noise has to be removed before continuing the analysis, the data were averaged over 4 × 4 km² areas. Then, for each map, the norm of the gradient was computed and normalized in order to reconstruct the cascade. Figure 5 (left) shows the scaling of the statistical moments for various orders. This set of scale laws is found to be consistent with the basic multifractal relation given by Eq. (8). For each order *q*, the slope of the scale law provides an estimation of *K*<sub>ζ</sub>(*q*). The moment scaling function retrieved by this method is shown in Fig. 6. The fit of this function according to Eq. (9) yields *C*<sub>1ζ</sub> ≈ 0.12 and *α*<sub>ζ</sub> ≈ 1.92 (here, we announced to provide the standard deviations of the estimators because they would indicate an artificially high precision; the estimation error of the whole analysis technique has been tested with simulations and is found to be around 10% for

both parameters). Note that the values of these parameters, as well as that of *H*, are close to those obtained for rain and clouds, which are respectively *H* ≈ 0.4, *C*<sub>1</sub> ≈ 0.12, *α* ≈ 1.8 (Verrier et al., 2010) and *H* ≈ 0.4, *C*<sub>1</sub> ≈ 0.08, *α* ≈ 1.9 (Lovejoy and Schertzer, 2006). Another possibility is to normalize the norm of the gradient in the same manner for all maps by using the “climatological” mean computed over all maps. This technique has the advantage to provide an estimation of the outer scale of the cascade by extrapolating the scale laws of the moments (see, e.g., Lovejoy and Schertzer, 2006). Although the sample used in this study is limited and may not be representative, the results are presented in Fig. 5 (right) and yield an outer scale equal to 2000 km, which could be related to the size of oceanic gyres in terms of order of magnitude.

We also tried to perform the same type of analysis using SST (Sea Surface Temperature), which is another useful, remotely sensed oceanic tracer. However, this attempt failed because the spectrum of the SST maps was found to flatten out at larger scales (around 32 km) than that of chlorophyll maps, and the available range of scales was thus insufficient. This whitening effect, which hides the small scale fluctuations, may be due to air-sea exchanges, which tend to spatially homogenize the SST. However, Nieves et al. (2007) performed a multi-scale analysis of SST data with a larger scale range (level L3 product) and found that the observed multifractal spectra was very similar to the one obtained with chlorophyll concentration data. This result provides an additional argument in favour of a link between phytoplankton patchiness and turbulent mixing at large scales, which will be developed in the next section.

The use of statistical moments is a very convenient way of estimating multifractal parameters. However, as it is not very intuitive, we propose here to demonstrate the existence of a cascade process, through the use of the more classical concept of probability density. The algorithm used in this method is the following: (i) compute the flux at the finest available scale, (ii) perform averages over 2 × 2 squares, (iii) compute the multiplicative weights that relate the values of the coarse-grained flux to the previous ones, (iv) plot the Probability Density Function (PDF) of the logarithm of these multiplicative weights, (v) iterate steps (ii), (iii) and (iv) until the largest scale of the cascade is reached. This method is straightforward to implement and does not require any prior assumption concerning the data. The results obtained with our selection of chlorophyll maps are given in Fig. 7. The PDF of the logarithm of the multiplicative weights does not depend on the scale at which they are derived, thus confirming the use of a scale invariant cascade model. Figure 8 shows the left tail of the total PDF, compared with a Gaussian distribution having the same mean and variance. The empirical PDF decays as a power law (producing a straight line on a log-log graphic), which is much slower than the Gaussian behaviour. This result supports the fact that the generator follows a Lévy law with infinite variance, and allows *α* to
<figure><img src="image_8.png" /><figcaption>Fig. 5. Scaling of the statistical moments of the flux ζ for the orders q = 0, 0.1, 0.2, ..., 2, with corresponding theoretical fits. Here, L corresponds to the largest scale of the SeaWiFS chlorophyll maps, i.e. 128 km. Left: for each map, the flux was normalized to a mean value of 1. Right: the normalization was performed with the “climatological” mean value computed over all maps, which allows estimating the outer scale of the cascade by extrapolation to larger scales.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Moment scaling function K_ζ(q) of the flux ζ, with theoretical fit.</figcaption></figure>

be estimated using a different approach, since the theoretical slope of the asymptote this distribution is equal to -(1+α). The resulting value of α is found to be 1.95, which is consistent with the value previously obtained using statistical moments.

## 7 Interpretation

Since the parameter *H* was found to be close to 1/3, it is tempting to relate it to the theory of passive scalars. This theory is based on the hypothesis of a 3-D isotropic turbulence that does not hold for our selection of chlorophyll maps, because, in the considered scale range (1–128 km), the ocean is a stratified fluid with a horizontal dimension much larger than the vertical one. However, some recent studies (e.g., Lovejoy and Schertzer, 2010) suggest that the Corrsin-Obukhov scale law may still be valid in the horizontal. Therefore, if turbulent mixing is the dominant effect, we may expect that the horizontal variability of phytoplankton fields would verify the scale law given in Eq. (2). If this is correct, then, assuming the velocity and passive scalar fluctuations to be independent, Schmitt et al. (1996) have shown that the parameter *H* of the FIF model (Eq. 9) should be equal to:

$$H = 1/3 + K_{\epsilon}(1/6) - K_{\chi}(1/2). \quad (20)$$

The deviation of *H* with respect to the value 1/3 is due to the intermittency of the energy and scalar variance fluxes, since a conserved flux raised to a power exponent, not equal to 1, is no longer a conserved quantity. The term K_ε(1/6) depends only on the turbulence, and is well known; by assuming the
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** PDFs of the logarithm of the multiplicative weights for each level of the cascade (corresponding to contractions of the averaging area by a factor 2², from 128 km² until 4 km²). The PDFs are very similar, with the exception of the function corresponding to the last scale contraction (from 8 km² to 4 km²), which is flatter. This may be due to the presence of noise.</figcaption></figure>

parameters $\alpha_{\epsilon} = 1.5$ and $C_{1\epsilon} = 0.25$ proposed by Schmitt et al. (1996), its value is expected to be around -0.05. However, the estimation of the term $K_{\chi}(1/2)$ is more delicate, since the multifractal parameters of $\chi$ are not known a priori, and have to be estimated. One possible solution consists in using the empirical multifractal parameters obtained for $\zeta$ in the previous section, because they have a simple approximate relationship to those of $\chi$ (de Montera et al., 2010):

$$ \begin{cases} \alpha_{\chi} \approx \alpha_{\zeta} \\ C_{1\chi} \approx 2^{\alpha_{\zeta}} C_{1\zeta}. \end{cases} \qquad (21) $$

This yields $\alpha_{\chi} \approx 1.92$ and $C_{1\chi} \approx 0.45$, thus allowing $K_{\chi}(1/2)$ to be estimated at a value equal to -0.11. The (semi-)theoretical value of $H$ is therefore $1/3 - 0.05 + 0.11 \approx 0.39$, which is consistent with its experimental value of 0.4 obtained with the SeaWiFS chlorophyll maps.

This coherency led us to the conclusion that phytoplankton behaves like a passive scalar within the studied scale range, which includes the mesoscale and the sub-mesoscale. This does not mean that phytoplankton is a purely passive scalar, however it implies that biological activity does not affect the scale law generated by turbulent mixing. This is consistent with the previous finding of Seuront et al. (1999) and Currie and Roff (2006), who showed that biological activity affected the scaling over a limited range only, between 30 m and 500 m, which is smaller than the resolution of remotely sensed satellite data.

However, as explained in the introduction, other studies (e.g., Lovejoy et al., 2001b) found a parameter $H$ equal

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Log-log graph of the left tail of the total PDF of the logarithm of multiplicative weights (blue), compared with a Gaussian having the same mean and variance (red). The PDF decays as a power law, with a slope -2.95 (green fit), corresponding to a Lévy law of index $\alpha = 1.95$. The Gaussian function decays much faster, and would therefore be inappropriate for cascade generation.</figcaption></figure>

to 0.12 and concluded to a combined turbulent/growth-dominated process. Therefore, the question is still open and future studies should try to understand precisely in which particular seasons or locations this departure from the turbulent scaling is likely to occur. According to the model proposed in Lovejoy et al. (2001a), this departure should be observed in area which have a weak turbulent activity combined with a high growth rate.

## 8 Bias in biogeochemical numerical models

The forecasting of coupled turbulent/biogeochemical systems is currently performed by means of 3-D numerical simulations. The main shortcoming of this technique is that it necessarily implies the use of high-pass filtering in scale space (or “scale truncation”), which strongly affects the estimation of non-linear advection terms in the fluid mechanics equations. This truncation of scale space is unavoidable because of the limited power of computers. It means, for example, that a small length interval, considered as a differential element $dx$ in the equations, has a value much larger than the scale of homogeneity in the numerical simulation (generally 10–100 km for global models, whereas dissipation occurs at scales of the order of a millimeter). The impact of this drastic simplification remains unknown. Although it is generally believed that it can be compensated for, for example by increasing the viscosity (Boussinesq hypothesis), this remains
<figure><img src="image_12.png" /><figcaption>Fig. 9. Assessment of the distribution of the relative error percentage resulting from the hypothesis of homogeneity over 128 km² areas, for a quadratic source term in a biogeochemical numerical model.</figcaption></figure>

to be demonstrated (for a test of the Boussinesq hypothesis, see Schmitt, 2007).

If biogeochemical processes are involved, the situation is even worse, because the estimation of these interactions is also affected by the truncation error. Moreover, the parameters of biogeochemical models are often obtained by means of laboratory experiments performed at a typical scale of one meter. Therefore, since the relations in which these parameters are involved are generally non-linear, it is not correct to use them at larger scales if the real fields are heterogeneous. It can thus be useful to assess the bias generated by the assumption of homogeneity over larger scales. For this, we consider a global numerical model operating with a 1° grid scale (roughly corresponding to the 128 km² maps analyzed in the present paper), which includes a quadratic source term of the form $βC²$, where C is the concentration of a tracer and β is a parameter assumed to be derived under stable conditions, at the scale of one meter in a laboratory. If it is assumed that the 100 SeaWiFS chlorophyll maps are realizations of the sub-grid heterogeneity of the tracer, then for each map we compute the source term at the finest available scale (which is 1 km in this case, whereas a 1 m scale would be needed!), and average these values over the whole 128 km² map. Finally, we estimate the value of this source term that would result from the hypothesis of homogeneity, by averaging the concentration over the whole 128 km² map and then computing the source term. The source term is then estimated with a relative error E equal to:

$$E = \frac{\langle \text{Chl}^2 \rangle - \langle \text{Chl} \rangle^2}{\langle \text{Chl}^2 \rangle} \qquad (22)$$

The PDF of the percentage of this relative error is shown in Fig. 9. Its mean value is approximately 22%, which is far from being negligible. One possible approach for reducing this error would be to derive an analytic expression for the scale dependency of the biological parameters (such as β in the example above), using the multifractal parameters of the tracer patchiness, if available.

## 9 Conclusions

Multifractal properties of oceanic chlorophyll maps have been investigated with remotely sensed data recorded from space. The FIF model has been validated, showing that chlorophyll maps can be modelled statistically, through the use of a fractionally integrated multiplicative cascade. In the study area, the Senegalo-Mauritanian upwelling region, the parameters of this model were found to be $H \approx 0.4$, $C_1 \approx 0.12$ and $\alpha \approx 1.92$. The estimates of the scale law exponent H is consistent with passive scalar behaviour, indicating that phytoplankton variability is dominated by turbulent mixing over the studied scale range (4–128 km), and that biological activity do not modify this scaling. This result confirms previous studies that reached this conclusion based on in-situ data. However, it cannot be generalized to other locations because it may not be correct in areas having a high growth rate combined with a weak turbulent activity.

Finally, it has been shown that, as a consequence of this multifractal patchiness, the non-linear source and sink of biogeochemical numerical models could be strongly biased. Future studies should therefore be dedicated to the use multifractal techniques to improve the accuracy of numerical simulations. This could be performed, for example, by predicting the scale dependence of the model parameters or by refining the assimilation of data measured at different scales. Although the effect of grazing was not observed in this study because of the low resolution of satellite data, the development of such techniques implies to take it into account since the scaling is modified at lower scales, in particular at scales of the order of the so-called “planktoscale”.

*Acknowledgements.* This study was funded by the French Centre National d'Études Spatiales (CNES).

Edited by: E. J. M. Delhez

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

The publication of this article is financed by CNRS-INSU.
## References

Aristegui, J., Alvarez-Salgado, X. A., Barton, E. D., Figueiras, F. G., Hernandez-Leon, S., Roy, C., and Santos, A. M. P.: Oceanography and fisheries of the Canary current/Iberian Region of the eastern North Atlantic, in: The Sea, Vol. 14, Chap. 23, edited by: Robinson, A. R. and Brink, K. H., John Wiley and Sons, New York, 877–931, 2004.

Corrsin, S.: On the spectrum of isotropic temperature fluctuations in an isotropic turbulence, J. Appl. Phys., 22, 469–475, 1951.

Currie, W. J. S. and Roff, J. C.: Plankton are not Passive Tracers: Plankton in a Turbulent Environment, J. Geophys. Res., 111, C05S07, doi:10.1029/2005JC002967, 2006.

de Montera, L., Verrier, S., Mallet, C., and Barthès, L.: A passive scalar-like model for rain applicable up to storm scale, Atmos. Res., 98(1), 140–147, 2010.

Denman, K. L. and Platt, T.: The variance spectrum of phytoplankton in a turbulent ocean, J. Mar. Res., 34, 593–601, 1976.

Dubrulle, B.: Intermittency in fully developed turbulence: Log-Poisson statistics and generalized scale-covariance, Phys. Rev. Lett., 73, 959–963, 1994.

Fasham, M. J. R.: The application of some stochastic processes to the study of plankton patchiness, in: Spatial Pattern in Plankton Communities, edited by: Steele, J. H., Springer, New York, 131–156, 1978.

Horwood, J. W.: Observations on spatial heterogeneity of surface chlorophyll in one and two dimensions, J. Mar. Biol. Assoc. UK, 58, 487–502, 1978.

Kolmogorov, A. N.: Local structure of turbulence in an incompressible liquid for very large Reynolds numbers, Proc. Acad. Sci. URSS, Geochem. Section, 30, 299–303, 1941.

Landau, L. D. and Lifshitz, E. M.: Fluid mechanics, 1st edn., Edt. MIR, Moscow, 1944.

Lathuilière, C., Echevin, V., and Levy, M.: Seasonal and intraseasonal surface chlorophyll-a variability along the northwest African coast, J. Geophys. Res., 113, C05007, doi:10.1029/2007JC004433, 2008.

Lavallée, D., Lovejoy, S., Schertzer, D., and Ladoy, P.: Nonlinear variability and landscape topography: analysis and simulation, in: Fractals in geography, edited by: de Cola, L. and Lam, N., Prentice-Hall, 171–205, 1993.

Lilley, M., Lovejoy, S., Strawbridge, K., and Schertzer, D.: 23/9 dimensional anisotropic scaling of passive admixtures using lidar aerosol data, Phys. Rev. E, 70, 036301–036307, 2004.

Lovejoy, S. and Schertzer, D.: Towards a new synthesis for atmospheric dynamics: space-time cascades, Atmos. Res., 96(1), 1–52, 2010.

Lovejoy, S. and Schertzer, D.: Multifractals, cloud radiances and rain, J. Hydrol., 322, 59–88, 2006.

Lovejoy, S., Currie, W. J. S., Tessier, Y., Claereboudt, M. R., Roff, J. C., Bourget, E., and Schertzer, D.: Universal Multifractals and ocean patchiness: phytoplankton, physical fields and coastal heterogeneity, J. Plankton Res., 23, 117–141, 2001a.

Lovejoy, S., Schertzer, D., Tessier, Y., and Gaonac'h, H.: Multifractals and Resolution independent remote sensing algorithms: the example of ocean colour, Int. J. Remote Sens., 22, 1191–1234, 2001b.

Nieves, V., Llebot, C., Turiel, A., Solé, J., García-Ladona, E., Estrada, M., and Blasco, D.: Common turbulent signature in sea surface temperature and chlorophyll maps, Geophys. Res. Lett.,

34, L23602, doi:10.1029/2007GL030823, 2007.

Novikov, E. A. and Stewart, R.: Intermittency of turbulence and spectrum of fluctuations in energy-dissipation, Izv. Akad. Nauk. SSSR, Ser. Geofiz., 3, 408–412, 1964.

Obukhov, A.: Structure of the temperature field in a turbulent flow, Izv. Akad. Nauk SSSR, 13, 55–69, 1949.

Panchev, S.: Random Functions and Turbulence, Pergamon Press, London, 1971.

Pecknold, S., Lovejoy, S., Schertzer, D., Hooge, C., and Malouin, J. F.: The simulation of universal multifractals, in: Cellular Automata: Prospects in astrophysical applications, edited by: Perdang, J. M. and Lejeune, A., World Scientific, 228–267, 1993.

Platt, T.: Local phytoplankton abundance and turbulence, Deep-Sea Res., 19, 183–187, 1972.

Pottier, C., Turiel, A., and Garçon, V.: Inferring missing data in satellite chlorophyll maps using turbulent cascading, Remote Sens. Environ., 112, 4242–4260, 2008.

Richardson, L. F.: Weather prediction by numerical processes, Cambridge University Press, Cambridge, 1922.

Schertzer, D. and Lovejoy, S.: Physically based rain and cloud modeling by anisotropic, multiplicative turbulent cascades, J. Geophys. Res., 92, 9692–9714, 1987.

Schertzer, D. and Lovejoy, S.: Universal Multifractals do Exist!, J. Appl. Meteorol., 36, 1296–1303, 1997.

Schertzer, D., Lovejoy, S., and Hubert, P.: An Introduction to Stochastic Multifractal Fields, in: ISFMA Symposium on Environmental Science and Engineering with related Mathematical Problems, edited by: Ern, A. and Weiping, L., Higher Education Press, Beijing, 106–179, 2002.

Schmitt, F., Schertzer, D., Lovejoy, S., and Brunet, G.: Multifractal temperature and flux of temperature variance in fully developed turbulence, Europhys. Lett., 34, 195–200, 1996.

Schmitt, F. G.: About Boussinesq's turbulent viscosity hypothesis: historical remarks and a direct evaluation of its validity, C.R. Mecanique, 335, 617–627, 2007.

Seuront, L., Schmitt, F. G., Lagadeuc, Y., Schertzer, D., Lovejoy, S., and Frontier, S.: Universal Multifractal structure of phytoplankton biomass and temperature in the ocean, Geophys. Res. Lett., 23, 3591–3594, 1996a.

Seuront, L., Schmitt, F. G., Schertzer, D., Lagadeuc, Y., and Lovejoy, S.: Multifractal Analysis of Eulerian and Lagrangian Variability of Physical and Biological Fields in the Ocean, Nonlinear Proc. Geoph., 3, 236–246, 1996b.

Seuront, L., Schmitt, F. G., Lagadeuc, Y., Schertzer, D., and Lovejoy, S.: Universal Multifractal analysis as a tool to characterize multiscale intermittent patterns: example of phytoplankton distribution in turbulent coastal waters, J. Plankton Res., 21, 877–922, 1999.

Seuront, L. and Schmitt, F. G.: Eulerian and Lagrangian properties of biophysical intermittency in the ocean, Geophys. Res. Lett., 31, L03306, doi:10.1029/2003GL018185, 2004.

Seuront, L. and Schmitt, F. G.: Multiscaling statistical procedures for the exploration of biophysical couplings in intermittent turbulence; Part I. Theory, Deep-Sea Res. Pt. II, 52, 1308–1324, 2005a.

Seuront, L. and Schmitt, F. G.: Multiscaling statistical procedures for the exploration of biophysical couplings in intermittent turbulence; Part II. Application, Deep-Sea Res. Pt. II, 52, 1325–1343, 2005b.
Steele, J. H. and Henderson, E. W.: Spatial patterns in North Sea plankton, Deep-Sea Res., 26, 955–963, 1979.

She, Z.-S. and Lévêque, E.: Universal scaling laws in fully-developed turbulence, Phys. Rev. Lett., 72, 336–339, 1994.

Verrier, S., de Montera, L., Barthès, L., and Mallet, C.: Multifractal analysis of African monsoon rain fields, taking into account the zero rain-rate problem, J. Hydrol., 389(1–2), 111–120, 2010.

Yaglom, A. M.: The influence on the fluctuation in energy dissipation on the shape of turbulent characteristics in the inertial interval, Sov. Phys. Dokl., 2, 26–30, 1966.