Multifractal analysis of oceanic chlorophyll maps remotely sensed from space

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. 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, 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.


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 Correspondence to: S. Thiria (sylvie.thiria@locean-ipsl.upmc.fr) 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(Seuront et al., , b, 1999Schmitt, 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 Published by Copernicus Publications on behalf of the European Geosciences Union. 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 biogeochemical numerical models.

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: v l ε 1/3 l 1/3 In this equation, v l ≡ |v(x + l) − v(x)| represents the mean shear of (longitudinal) velocity between two points separated by a distance l, and ε 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): where ϕ ≡ χ 3/2 ε −1/2 represents the non-linear coupling between velocity and concentration, with χ ≡ − ∂( C l ) 2 ∂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): where ε l and χ l are the fluxes averaged at scale l, L is the largest scale of the cascade, and K ε (q) and K χ (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 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): where α ε and α χ are multifractality parameters varying between 0 and 2, and C 1ε and C 1χ are intermittency parameters varying between 0 and the dimension of the embedding space, which here is equal to 2.

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.: Chl l ζ a l H where Chl is the chlorophyll concentration and Chl l ≡ |Chl(x + l) − Chl(x)| . ζ is a conserved flux and a and H are adjustable parameters. As described above, the conserved flux has to verify the basic multifractal relation: and it is assumed that it converges towards a log-stable law: This model is described by four parameters: a, H , α ζ and C 1ζ . However, it is possible to reduce the number of parameters to three, since taking the a-th power of ζ in Eq. (7) is equivalent to a simple shift of H by K ζ (a), and to the multiplication of C 1ζ by a factor a α ζ (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 Chl q l ≡ |Chl(x +l)−Chl(x)| q . According to Eq. (7), the q-th-order structure function is equal to: Chl q l ζ aq l qh .
(10) By introducing Eq. (8), this equation simplifies to: Then, the term K ζ (aq) can be decomposed into conservative and non-conservative parts using the following identity, which can be straightforwardly derived from Eq. (9): This operation yields: As expected, by defining: one obtains the usual form of the structure function corresponding to the case where a = 1: In the following, we thus use a simplified form of Eq. (7) requiring only three parameters, namely, H , α ζ and C 1ζ : 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.

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 ζ is assumed to be conserved in scale space, whatever the scale l, ζ l is constant. Therefore, Eq. (9) reduces to: This equation allows H to be estimated using the simple expression: The second step consists in quantifying the multifractal properties of the flux, and in estimating α ζ and C 1ζ . 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: 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 ζ (q) can be estimated. Finally, the parameters α ζ and C 1ζ are obtained by determining the least squares fit to this function.

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 higherorder 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 Schertzer, 2006 andLilley 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 2 . However, the spot size over which the data are measured varies with Ocean Sci., 7, 219-229, 2011 www.ocean-sci.net/7/219/2011/ 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 2 , 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 (1 km) (2 km) (4 km) (8 km) (64 km) (128 km) (16 km) (32 km) 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 theorical 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). isolated, unrealistically high values. Each of these unrealistic pixels was corrected using the mean value of the surrounding pixels. 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, 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).

Results
Since the noise has to be removed before continuing the analysis, the data were averaged over 4 × 4 km 2 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 ζ (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 1ζ ≈ 0.12 and α ζ ≈ 1.92 (here, we renounced 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 1 ≈ 0.12, α ≈ 1.8  and H ≈ 0.4, C 1 ≈ 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 Ocean Sci., 7, 219-229, 2011 www.ocean-sci.net/7/219/2011/ 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.

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: 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  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 2 , from 128 km 2 until 4 km 2 ). The PDFs are very similar, with the exception of the function corresponding to the last scale contraction (from 8 km 2 to 4 km 2 ), which is flatter. This may be due to the presence of noise. parameters α ε = 1.5 and C 1ε = 0.25 proposed by Schmitt et al. (1996), its value is expected to be around −0.05. However, the estimation of the term K χ (1/2) is more delicate, since the multifractal parameters of χ are not known a priori, and have to be estimated. One possible solution consists in using the empirical multifractal parameters obtained for ζ in the previous section, because they have a simple approximate relationship to those of χ (de : This yields α χ ≈ 1.92 and C 1χ ≈ 0.45, thus allowing K χ (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 ≈ 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 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 α = 1.95. The Gaussian function decays much faster, and would therefore be inappropriate for cascade generation. to 0.12 and concluded to a combined turbulent/growthdominated 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.

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  Fig. 9. Assessment of the distribution of the relative error percentage resulting from the hypothesis of homogeneity over 128 km 2 areas, for a quadratic source term in a biogeochemical numerical model. 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 2 maps analyzed in the present paper), which includes a quadratic source term of the form βC 2 , 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 2 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 2 map and then computing the source term. The source term is then estimated with a relative error E equal to: 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.

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 ≈ 0.4, C 1 ≈ 0.12 and α ≈ 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".