Reconstruction of super-resolution ocean pCO 2 and air-sea ﬂuxes of CO 2 from satellite imagery in the Southeastern Atlantic

. An accurate quantiﬁcation of the role of the ocean as source/sink of Green House Gases (GHGs) requires to access the high-resolution of the GHG air-sea ﬂux at the interface. In this paper we present a novel method to reconstruct maps of surface ocean partial pressure of CO 2 (pCO 2 ) and 5 air-sea CO 2 ﬂuxes at super resolution (4 km) using Sea Surface Temperature (SST) and Ocean Colour (OC) data at this resolution, and CarbonTracker CO 2 ﬂuxes data at low resolution (110 km). Inference of super-resolution of pCO 2 , and air-sea CO 2 ﬂuxes is performed using novel nonlinear signal 10 processing methodologies that prove efﬁcient in the context of oceanography. The theoretical background comes from the Microcanonical Multifractal Formalism which unlocks the geometrical determination of cascading properties of physical intensive variables. As a consequence, a multiresolution 15 analysis performed on the signal of the so-called singularity exponents allows the correct and near optimal cross-scale inference of GHGs ﬂuxes, as the inference suits the geometric realization of the cascade. We apply such a methodology to the study offshore of the Benguela area. The inferred rep-20 resentation of oceanic partial pressure of CO 2 improves and enhances the description provided by CarbonTracker, capturing the small scale variability. We examine different combinations of Ocean Colour and Sea Surface Temperature products in order to increase the number of valid points and the 25 quality of the inferred pCO 2 ﬁeld. The methodology is validated using in-situ measurements by means of statistical errors. We obtain that mean absolute and relative errors in the inferred values of pCO 2 with respect to in-situ measurements are smaller than for CarbonTracker.


Introduction
The ocean can be thought of as a complex system in which a large number of different processes (e.g.physical, chemical, biological, atmosphere-ocean interactions) interact with each 35 other at different spatial and temporal scales (Rind, 1999).These scales extend from millimeters to thousands of kilometers and from seconds to centuries (Dickey, 2003).Recently there is a growing body of evidence that the upper few hundred meters of the oceans are dominated by submesoscale 40 activity, covering the range 1-10 km, and that this activity is important to understand global ocean properties (Klein and Lapeyre, 2009).Accurately estimating the sources and sinks of GHGs at the air-sea interface requires to resolve these small scales (Mahadevan et al., 2004).However, the scarcity of oceanographic cruises and the lack of available satellite products for GHG concentrations at high resolution prevent us from obtaining a global assessment of their spatial variability at small scales.For example, from the in-situ ocean measurements the uncertainty of the net global ocean-50 atmosphere CO 2 fluxes is between 20 and 30% (IOCCP, 2007), and could be higher in the Oxygen Minimum Zones (OMZ) of the Eastern Boundary Upwelling Systems (EBUS) due to the extreme regional variability in these areas (Paulmier et al., 2008;Franco et al., 2014).This indeed suggests 55 the design of proper methodologies to infer the fluxes at high resolution from presently available satellite image data, in order to improve current estimates of gas exchanges between the ocean and the atmosphere.
The most commonly used methods to estimate air-sea CO 2 fluxes are based either on statistical methods, inverse modeling with atmospheric transport models or global coupled physical-biogeochemical models.Among others, Takahashi et al. (2002Takahashi et al. ( , 2009) ) that interpolate sea surface pCO 2 measurements with advanced statistical methods to provide cli-65 matological monthly maps of air-sea fluxes of CO 2 in the global surface waters at a spatial resolution of 4 • × 5 • .Global maps at the same spatial resolution but at higher temporal resolution (daily) have been estimated by Rödenbeck et al. (2014) by fitting the mixed-layer carbon budget equation to ocean pCO 2 observations.An international effort to compile global surface CO 2 fugacity (fCO 2 ) measurements has recently been performed and reported in Pfeil et al. (2013); Bakker et al. (2014), and interpolated by Sabine et al. (2013) generating a monthly gridded product with fCO 2 values in a 1 • × 1 • grid cell.Other statistical approaches based on the neural-network statistical method have been shown to be useful to estimate climatological and monthly 1 • × 1 • maps of pCO 2 by Landschützer et al. (2014) and Telszewski et al. (2009) respectively.Gruber et al. (2009) used inverse modeling of sources and sinks from the network of atmospheric CO 2 concentrations jointly with transport models.The third type of methods is based on direct computations of the air-sea CO 2 fluxes in coupled physical-biogeochemical models incorporating the biogeochemical processes of the carbon dioxide system.In the latter, simulated surface ocean pCO 2 can be constrained with available ship observations as shown by Valsala and Maksyutov (2010).
Another new avenue to infer air-sea GHG fluxes is through inverse modeling applied to vertical column densi-90 ties (VCD) extracted from satellite spectrometers, i.e.Greenhouse gases Observing SATellite (GOSAT) and SCanning Imaging Absorption SpectroMeter for Atmospheric CHar-tographY (SCIAMACHY), at low spatial resolution (Garbe and Vihharev, 2012).A global estimation of CO 2 fluxes in 95 the ocean has been derived at 1 • × 1 • spatial resolution from global atmosphere observations used into a data assimilation system for CO 2 called CarbonTracker (Peters et al., 2007).In all these datasets the rather coarse spatial resolution leads to uncertainties in the actual estimate of the sources and sinks 100 of CO 2 , calling for an improvement of the resolution of CO 2 flux estimates.
In this regard, the last few years have seen the appearance of interesting new developments on multiscale processing techniques for complex signals coming from Earth Observations (Yahia et al., 2010).These methods make use of phenomenological descriptions of Fully Developed Turbulence (FDT) in nonlinear physics, motivated by the values taken on by Reynolds number in ocean dynamics.As predicted from the theory and also observed in the ocean, in 110 a turbulent flow the coherent vortices (eddies) interact with each other stretching and folding the flow, generating smaller eddies or small scale filaments and transition fronts characterized by strong tracer gradients (Frisch, 1995).This results in a cascade of energy from large to smaller scales.There-115 fore the inherent cascade of tracer variance under the turbulent flow dominates the variability of the geometrical distribution of tracers such as temperature or dissolved inorganic carbon, as shown by Abraham et al. (2000), Abraham and Bowen (2002) and Turiel et al. (2005).Geometrical organi-120 zation of the flow linked to the energy cascade allows study of its properties from the geometrical properties of any tracer for which the advection is the dominant process.The relationships between the cascade and the multifractal organization of FDT has been set up either in a canonical (Arneodo 125 et al., 1995;Frisch, 1995) or microcanonical (Turiel et al., 2005;Bouchet and Venaille, 2012) descriptions.Within the microcanonical framework (MMF) the singularity exponents unlock the geometrical realization of the multifractal hierarchy.Setting up a multiresolution analysis on the singular-130 ity exponents computed in the microcanonical framework allows near optimal cross scale inference of physical variables (Sudre et al., 2015).
These advances open a wide field of theoretical and experimental research and their use in the analysis of complex data 135 coming from satellite imagery has been proven innovative and efficient, showing a particular ability to perform fusion of satellite data acquired at different spatial resolutions (Pottier et al., 2008) or to reconstruct from satellite data current maps at submesoscale resolution (Sudre et al., 2015).In this 140 paper we apply these novel techniques emerging from nonlinear physics and nonlinear signal processing for inferring submesoscale resolution maps of the air-sea CO 2 fluxes and associated sinks and sources from available remotely sensed data.We use this methodology to derive cross scale inference 145 according to the effective cascade description of an intensive variable, through a fusion process between appropriate physical variables which account for the flux exchanges between the ocean and the atmosphere.This approach is not only very novel in signal processing, but also connects the 150 statistical descriptions of acquired data with their physical content.This makes the approach useful to reconstruct all GHGs.
Unlike the Lagrangian approach to reconstruct tracer maps at high resolution (Berti and Lapeyre, 2014), our methodol-155 ogy works in the Eulerian framework and we do not need to know the trajectories of oceanic tracer particles but only high resolution instantaneous maps of tracers which can be directly obtained from remote sensing.
The Eastern Boundary Upwelling Systems (EBUS) and 160 Oxygen Minimum Zones (OMZs) are likely to contribute significantly to the gas exchange between the ocean and the atmosphere (Hales et al., 2005;Waldron et al., 2009;Paulmier et al., 2011).The Benguela upwelling system, the region of interest in this study, is one of the highest produc-  icantly to the global air-sea CO 2 flux.Some studies using data from in-situ samples have found the region of Benguela to be an annual sink of CO 2 with -1.70 (in 1995 and 1996) and -2.02Mt C/year in 2005 (Santana-Casiano et al., 2009;Monteiro, 2010), with a strong variability between 2005 and 2006 from -1.17 to -3.24 mol C/m 2 per year, respectively (González-Dávila et al., 2009).This paper is organized as follows: Sect. 2 describes the datasets used as input in our algorithm.Sect. 3 is devoted to describe the methodology used through the study.Statistical description of the input datasets is presented in Sect. 4. Results of the inference method are given in Sect. 5 by providing outputs of our algorithm, then evaluating the various satellite products and assessing the performance of the method using in situ measurements.

Data
The input data combines air-sea CO 2 fluxes at low resolution and satellite ocean data at high resolution.To validate the method we use in-situ measurements of oceanic pCO 2 .
2.1 Input data: Air-sea CO 2 fluxes at low resolution It is known that the evolution of a concentration, c, in the atmosphere is given by the advection-reaction-diffusion equation: 190 with the wind field u, the density of the air ρ, the turbulent diffusivity tensor T d , the chemical reaction rate g and the net flux at the air-sea interface F (Garbe et al., 2007(Garbe et al., , 2014)).Using optimal control and inverse problem modeling, a map of F can be derived using Earth Obser-195 vation data (Garbe and Vihharev, 2012) (Peters et al., 2007).CarbonTracker system assimilates and integrates a diversity of atmospheric CO 2 data into a computation of sur-210 face CO 2 fluxes, using a state-of-the-art atmospheric transport model and an ensemble Kalman filter.We obtain the partial pressure of ocean CO 2 by using the equation of the net flux in the air-sea interface: where α is the gas solubility, which depends on SST and Sea Surface Salinity SSS, and K, the gas transfer velocity, is a function of wind, salinity, temperature, sea state, which can be obtained from satellite data.To estimate the gas transfer velocity we use the well accepted relationships for the rived from LEGOS (Laboratoire d'Etudes en Géophysique et Océanographie Spatiales) product compiled by Delcroix et al. (2011) and winds from Cross Calibrated Multi-Platform Ocean surface winds from JPL (Jet Propulsion Laboratory) PO.DAAC (Physical Oceanography Distributed Active Archive Center, http://podaac.jpl.nasa.gov/).We assume a p air CO2 to be constant in the domain of study and it is derived from the Globalview-CO2 product of the Cooperative Atmospheric Data Integration Project coordinated by Carbon Cycle Greenhouse Gases Group (GLOBALVIEW, 2013) (www.esrl.noaa.gov/gmd/ccgg/globalview/).We use values taken at the closest sea-level station to the Benguela, located at Ascension Island (7.97 • S and 14.40 • W) as our reference atmospheric CO 2 .The partial pressure of CO 2 is determined from its mole fraction (xCO 2 ), using the following equation (Dickson et al., 2007): pCO 2 = xCO 2 •p, where p is the total pressure of the mixture.We assume this pressure to be close to 1 atm for the conversion following ORNL/CDIAC 105 report (Program developed for CO 2 Systems calculation).
The raw data of CarbonTracker fluxes of CO 2 in the area of interest are strongly binned and exhibit strong gradients across those bins.This turns out to be suboptimal for our super-resolution approach.Garbe and Vihharev (2012) have developed an optimal control approach to invert interfacial fluxes using a simplified inverse problem of at-250 mospheric transport.The inverse problem is solved using the Galerkin finite element method and the Dual Weighted Residual (DWR) method for goal-oriented mesh optimization.An adaptation of this approach has been applied to the CarbonTracker data set.However, the estimations are expensive and computing results for all the time frames of interest was infeasible.Therefore, an anisotropic diffusion-based approach has been applied to the raw fluxes of the Carbon-Tracker data set.The diffusion is steered by the direction of the low-altitude wind field.The results thus retain the structure of the CarbonTracker fluxes very well while suppressing artifacts.Examples of this process are shown in Fig. 1.
Results are comparable to the physically more accurate approach of Garbe and Vihharev (2012).
2.2 Input data: Satellite Ocean data at high resolution 265 Oceanic pCO 2 is a complex signal depending, at any spatial resolution, on sea surface temperature, salinity, chlorophyll concentration, dissolved inorganic carbon, alkalinity and nutrients concentrations.Both the biological pump, with chlorophyll a as a proxy, and the physical pump, driven by 270 the temperature and salinity (e.g.solubility, water mass), govern the evolution of pCO 2 in the surface ocean.
We use here the high resolution satellite ocean data for chlorophyll a, as a proxy for the biological carbon pump and for Sea Surface Temperature (SST), as a proxy for the ther-275 modynamical pump, (see Section 3.2 for more details on the connection of these oceanic variables).

Chlorophyll-a (Chl-a) from Ocean Colour (OC)
In this study we use Chl-a concentrations from two different Ocean Colour products: MERIS and GLOBCOLOUR.-Averaging from single-instrument chl-a concentration.In this case CHL1 daily level 3 (L3) products 295 are generated for each instrument using the corresponding L2 data.At the beginning of the averaging process, an inter-calibration correction is applied to the MODIS and SeaWiFS (Sea-Viewing Wide Field-of-View Sensor) CHL1 daily L3 products in order to get compatible 300 concentrations with respect to the MERIS sensor.The merged CHL1 concentration is then computed as the average of the MERIS, MODIS and SeaWiFS quantities, both as: an arithmetic mean or a weighted average value (AVW).In the AVW method, values of CHL1 are 305 weighted by the relative error for each sensor on the results of the simple averaging.
-Garver-Siegel-Maritorena model (GSM).In this method single-instrument daily L3 fully normalized water leaving radiances (individually computed for each 310 band) and their associated error bars are used by the GSM model.These radiances are not inter-calibrated before incorporation in the model (see Maritorena and Siegel (2005) for more details).
Snapshots of both Chl-a fields derived from MERIS and GSM GLOBCOLOUR corresponding to September 21, 2006 are displayed in Fig. 2 a) and b), respectively.This example shows the clear difference in the remote sensing coverage between the two products.The merged GLOBCOLOUR product yields a more covered Chl field than the one obtained from MERIS.The merging algorithm in GLOBCOLOUR product tends to decrease the missing points induced by clouds for each individual instrument.

Sea Surface Temperature (SST)
We use SST derived from OSTIA and MODIS products.
OSTIA (Operational SST and Sea Ice Analysis system) is a new analysis of SST that uses satellite data provided by the GHRSST (Group for High Resolution SST) project, together with in situ observations to determine the SST with a global coverage and without missing data.The datasets are produced daily and at spatial resolution of 1/20 • (∼6km) performing a multi-scale optimal interpolation using correlation length scales from 10 km to 100 km (more details in Donlon et al. (2012)).The other SST product used in this study is derived from MODIS (MODerate Resolution Imaging Spectroradiometer) sensors carried on board the Aqua satellite since December 2002.This SST product is derived from the MODIS mid-infrared (IR) and thermal IR channels and is available in various spatial and temporal resolutions.We use Level-3 daily maps of SST at the spatial resolution of 1/24 • (∼4 km) (Savtchenko et al., 2004).In Fig. 3 a) and b), we show one snapshot of SST from OSTIA and MODIS respectively corresponding to the same day on September 21, 2006.In the case of OSTIA products, the SST field is fully covered of points while for MODIS products there are gaps due to cloudiness.On other hand, MODIS product offers a more detailed visualization of the small structures.All SST data have been regridded at 1/32 • by bilinear interpolation.

Validation data: in-situ measurements
Among the available data in SOCAT version 2 (Bakker et al., 2014)  QUIMA2008-5, QUIMA2008-6, QUIMA2008-7 -2010, one cruise: ANT27-1 The small number of cruises found in one decade (24 cruises) shows that the scarcity of cruise in the Benguela 370 region is a fact.This indeed demonstrates the crucial need of developing a robust method to infer high resolution pCO 2 from space.Moreover for some of these cruises, for instance, the track of GALATHEA cruise is too close to the coast and is out of the original CarbonTracker domain.Due to this re-

Method
The idea behind the methodology hinges on the fundamental discovery of a simple functional dependency between the transitions -those being measured by the dimensionless values of the singularity exponents computed within the framework of the Microcanonical Multifractal Formalismof the respective physical variables under study : SST, Ocean Colour and oceanic partial pressure (pCO 2 ).With that functional dependency being adequately fitted into a linear re-400 gression model, it becomes possible to compute, at any given time, a precise evaluation of pCO 2 singularity exponents using SST, Ocean Colour and low resolution acquired pCO 2 .Once these singularity exponents are computed, they generate a multiresolution analysis from which low resolution 405 pCO 2 can be cross-scale inferred to generate a high resolution pCO 2 product.In this study we choose SST and Chl, and not other variables such as Sea Surface Height, because we focus on the use of physical variables which are correlated spatially and temporally to pCO 2 and that can be obtained

Singularity exponents and the multifractal hierarchy of turbulence
In the ocean, the turbulence causes the formation of unsteady eddies on many scales which interact with each other (Frisch, 1995).Most of the kinetic energy of the turbulent motion is contained in the large scale structures.The energy cascades from the large scale structures to smaller scale structures by an inertial and essentially inviscid mechanism.This process continues, creating smaller and smaller structures which produces a hierarchy of eddies.Moreover, the ocean is a system displaying scale invariant behavior, that is, the correlations of variables do not change when we zoom in or we zoom out the system, and can be represented by power-laws in particular, with the scaling exponents h.
It can be shown that the scaling exponents are the values taken on by localized singularity exponents, which can be computed at high precision in the acquired data using the Microcanonical Multifractal Formalism.Hence, within that framework, the multifractal hierarchy of turbulence, defined 430 by a continuum of sets F h indexed by scaling exponents h, is obtained as the level sets of the geometrically localized singularity exponents.
We will not review here the details in the computation of the singularity exponents h(x), leaving the reader to consult 435 references (Turiel et al., 2005(Turiel et al., , 2008;;Pont et al., 2011b;Maji and Yahia, 2014;Sudre et al., 2015) for an effective description of an algorithm able to compute the h(x) at every point x in a signal's domain.

Functional dependencies between the singularity exponents of intensive physical variables
Another important idea implemented in the methodology is 450 the coupling of the physical information contained in SST and OC images with the ocean pCO 2 .For instance, it is known, that marine primary production is a key process in the oceanic carbon cycling, and variations in the concentration of phytoplankton biomass can be related to variations 455 in the carbon concentrations.Surface temperature is also related to the gas solubility in the ocean, and areas with high temperatures are more suitable for releasing CO 2 to the at-mosphere.We have studied the relationship of SST and Chla variables with pCO 2 using the outputs of a coupled Regional Ocean Modeling System (ROMS) with the BIOgeochemical model of the Eastern Boundary Upwelling System (BIOEBUS) (Gutknecht et al., 2013).The ROMS includes several levels of nesting and composed grids, which makes it an ideal model for the basis of our methodology in 465 working in two spatial resolutions.BIOEBUS has been developed for the Benguela to simulate the first trophic levels of the Benguela ecosystem functioning and also to include a more detailed description of the complete nitrogen cycle, including denitrification and anammox processes as well as the oxygen cycle and the carbonate system.This model coupled to ROMS has been also shown to be skillful in simulating many aspects of the biogeochemical environment in the Peru upwelling system (Montes et al., 2014).When one compares SST and Chl with pCO 2 one finds undetermined 475 functional dependency.However, when comparing their corresponding singularity exponents one obtains a clear simpler dependency.This is due to the fact that SST, Chl and pCO 2 are variables of different dimensions while singularity exponents are dimensionless quantities.These results show that there is a good correlation between the turbulent transitions given by the singularity exponents and that singularity exponents are good candidates for a multiresolution analysis performed on the three signals SST, Chl and pCO 2 .Furthermore, the log-histograms and singularity spectrum show that singularity exponents of pCO 2 images possess a multifractal character.Therefore, such signals are expected to feature cascading, multiscale and other characteristic properties found in turbulent signals as described in Turiel et al. (2008) and Arneodo et al. (1995).Consequently the use of non-linear and multiscale signal processing techniques is justified to assess the properties of the pCO 2 signal along the scales.Therefore, in our methodology, the local connection between different tracer concentrations (SST and Chl-a) with pCO 2 , is performed in order to obtain a proxy for pCO 2 at high resolution by using the following linear combination of multiple linear regressions: where S(pCO 2 )(x) refers to the singularity exponent of pCO 2 at x, S(SST )(x) to singularity exponent of SST at x, S(Chl − a)(x) to singularity exponent of Chl-a signal at x.In order to propagate the pCO 2 signal itself along the scales in the multiresolution analysis we introduce S(pCO 2 LR ) to refer to the singularity exponent from pCO 2 at low resolution interpolated on the high resolution grid.a(x), b(x) and c(x) are the regression coefficients associated to singularity exponents, and d(x) is the error associated to the multiple-linear regression.These regression coefficients are estimated using simulated data from the ROMS-BIOEBUS model developed for the Benguela upwelling system and described above.

510
Once we have introduced these coefficients in the linear combination on satellite data, we obtain a proxy for singularity exponents of pCO 2 at high resolution and we can perform the multiresolution analysis to infer the information across the scales.interesting are those which are optimal with respect to inferring information along scales, in particular, in a context where information is to be propagated along the scales from low resolution to high resolution.
The effective determination of an optimal wavelet for a 530 given category of turbulent signals is, in general, a very difficult open problem.This difficulty can be contoured by considering multiresolution analysis performed on the signal of the singularity exponents h(x) themselves.Indeed, since the most singular manifold (the set F h associated to the low-535 est singularity exponents) is associated with the highest frequencies in a turbulent signal, and since the multifractal hierarchy F h converges to this set, it is physically evident that the multifractal hierarchy corresponds to a description of the detail spaces of a multiresolution analysis performed on a 540 turbulent signal.Consequently, designing by V j and W j respectively the approximation and detail spaces computed on S(pCO 2 )(x) signal, and by A j and P j their corresponding orthogonal projections from space L 2 (R 2 ), the following reconstruction formula: consists in reconstructing a signal across the scales using the detail spaces of the singularity exponents, hence regenerating a physical variable according to its cascade decomposition.From these ideas, which are described more 550 fully in the paper by Sudre et al. (2015), we can deduce the following algorithm for reconstructing a super-resolution pCO 2 signal from available high-resolution SST, Chl-a, and low-resolution pCO 2 : i) After selecting a given area of study, compute the singu-555 larity exponents of SST, Chl and pCO 2 at low and high resolution from ROMS-BIOEBUS output.This is done once and then they can be used for every computation performed over the same area.
-K: gas transfer velocity obtained by the parametrization developed by Sweeney et al, 2007, as a function of the wind.
iii) Obtain the regression coefficients a, b, c and d of Eq. 3 for the singularity exponents obtained in step ii) iv) Calculate the singularity exponents of available satellite SST, Chl at high resolution and ocean pCO 2 at low resolution (step i).
v) Use coefficients obtained in step iii) and apply Eq. 3 to the singularity exponents from satellite data (step iv) to estimate a proxy of singularity exponents of high resolution ocean pCO 2 , S(pCO 2 ).
vi) Using Eq. 4 reconstruct pCO 2 at high resolution from the multiresolution analysis computed on signal S(pCO 2 ) and cross-scale inference on pCO 2 at low resolution.
vii) Use Eq. 2 to calculate air-sea CO 2 fluxes from the inferred pCO 2 obtained in step vi) The methodology has been successfully applied to dual ROMS simulation data at two resolutions, obtaining a mean absolute error of pCO 2 reconstructed values with respect to ROMS simulated high-resolution pCO 2 equal to 3.2µatm (0.89% of relative error) (V.Garc ¸on 2014, pers.comm.).

Preliminary analysis of Sea Surface Temperature (SST) and Chlorophyll images
Since the key element for the application of our inferring algorithm relies on the ability to obtain the singularity exponents and their quality, the success of our methodology applied to satellite data depends on the quality and the properties of the input data.In order to assess such properties we perform a statistical analysis of the different datasets.First, we analyze the Chl and SST Probability Distribution Functions (PDFs).In Fig. 4a) we present the PDFs for Chl from MERIS, GLOBCOLOUR-GSM and GLOBCOLOUR-AVW; the required histograms are built using daily Chl values over 2006 and 2008 at each point of the spatial grid in the area of Benguela.Each one of these PDFs is broad and asymmetric, with a small mode (i.e. the value of Chl at which the probability reaches its maximum) between 0.1 and 0.2 mg/m 3 and a heavy tail.The heavy tail (i.e.non-gaussianity) means that the extreme values can not be neglected.In this case Chl values are mostly low (small mode) but there is a significant number of isolated and dispersed patches with very high Chl values producing intermittency (long tails in the PDF).Intermittency in the context of turbulence is the 610 tendency of the probability distributions of some quantities to develop long tails, i.e. the occurrence of very extreme events.Further information can be obtained by computing statistical quantities such as standard deviation, skewness and kurtosis.Table 1 shows that standard deviation is rather the same 615 for the three OC products while skewness and kurtosis values hugely differ.The degree of intermittency is measured by the kurtosis, the higher the kurtosis, the higher the intermittency.We found that kurtosis is almost ten times higher in GLOBCOLOUR products than in MERIS.

620
We have repeated the same analysis for SST datasets.The PDFs of the SST values for OSTIA and MODIS products are shown in Fig. 4b).In this case both PDFs possess similar shape, broad with the mode around 18 • C with a much less deviation from gaussianity as compared to Chl values.This 625 is confirmed with the computation of the statistical moments showed in Table 1.We obtain small values of the standard deviation and kurtosis in both cases, although slightly higher in the case of MODIS.The kurtosis is less than 3, meaning that there is not an important number of atypical values of 630 SST and therefore weak and short tails in the PDFs.If turbulence is dominated by coherent structures localized in space and time, then PDFs are not Gaussian, and the kurtosis will be higher than 3. To analyze this feature we turn to the statistical analysis of the singularity exponents, which, 635 as explained before, have the ability to unveil the cascade structures given by the tracer gradients.In Fig. 4c), it can be seen that the PDFs of the singularity exponents of the Chl for the three products are rather similar with almost the same standard deviation and with a slightly higher value of 640 the kurtosis in the GLOBCOLOUR-GSM product, 4.3, than for MERIS, 3.1, and GLOBCOLOUR-AVW, 3.1, (see Table 2).This shows that Chl from GLOBCOLOUR-GSM product contains more extreme values which produce intermittency likely given by the strongest structures.The PDFs of 645 the singularity exponents of the SST for OSTIA is narrower and with a highest peak than for MODIS SST.However, surprisingly the kurtosis is larger for singularity exponents of OSTIA SST, 5.1, than for MODIS SST, 3.2.Finally, we obtain the singularity spectra from the em-650 pirical distributions of singularity exponents shown in Fig. 4c) and d).One can see in Fig. 4e) that for the two GLOB-COLOUR products the shape of the spectrum is closer to binomial cascade of multiplicative processes than for MERIS.This is discussed in more depth in the next sections.We now apply the methodology to infer ocean pCO 2 maps at For the inference we use the following three combinations of Chl and SST products described in Section 2.1: MERIS-OSTIA, GLOBCOLOUR-OSTIA, GLOBCOLOUR-MODIS.We do not include the MERIS-MODIS combination in the analysis due to the fact that the use of such satellite data results in a too drastic reduction of the coverage of the resulting pCO infer 2 field, but using merged products offers wider coverage instead.The inferred pCO 2 obtained from two merged products for Chl-a, GLOBCOLOUR GSM and GLOBCOLOUR AVW is rather the same, with a slight improvement when GSM is used.Thus for the sake of clarity, we only show Figures for GLOBCOLOUR-GSM and some statistical results making comparisons with AVW.Therefore from now on we refer GLOBCOLOUR to the Chl-a obtained by the GSM merged method.Figure 5d shows one example of pCO infer 2 field corresponding to March 22, 2006 when we use SST data from OSTIA (Fig. 5a), Ocean Colour from GLOBCOLOUR (Fig. 5b) at high resolution and pCO 2 at low resolution (Fig. 5c) derived from CarbonTracker air-sea flux of CO 2 (Fig. 5e) and using the Eq. 2. The air-sea flux of CO 2 at super-resolution (Fig. 5f) is obtained from the pCO infer 2 field and a constant value of atmospheric pCO 2 equal to 385.6µatm.On this day the images of the pCO infer 2 and fluxes of CO 2 combine a good coverage and a clear identification of small scale structures and gradients, as described below.Note that the air-sea CO 2 flux from CarbonTracker presents a large land mask close to the coast and consequently, we rather study the offshore area of the Benguela upwelling.Comparing the figures one can see that values of pCO 2 and CO 2 flux over the domain (from 4.5 • E to coast (taking out the mask of the CarbonTracker domain and from 20.5 • S to 35 • S) vary between 360 and 380µatm and between -4x10 −8 and 0.5x10 −8 mol C m −2 s −1 , respectively.The resultant flux of CO 2 is positive (to-wards the atmosphere) in the region 25 • -28 • S and from 7 • E eastward to the coast and is negative (into the ocean) south of 30 • S and east of 6 • E. Thus, we see that in the southern part offshore of the Benguela area there is a strong CO 2 sink and the northern part behaves as a weak CO 2 source.

710
What is new in the reconstructed pCO 2 is, for instance, that the cascade of information across the scales enhances gradients in the field of pCO 2 .It is striking that the highresolution map provides the position of the North-South dipole "front" located at 30 • S (i.e.-1.5x10 −8 isoline in 715 green) which could not be inferred accurately from the low resolution map.The low resolution map would provide an estimate of the location of the "front" that is ∼1.5 • northern of the location inferred from the high-resolution map.Moreover one can see small structures in the pCO infer 2 field 720 between 33-35 • S and 9-12 • E in the pCO infer 2 field (Fig. 5d) .The small spatial scale variability is captured in the superresolution pCO 2 field and not in pCO Ctrack 2 as shown in the longitudinal profile of the images plotted in Fig. 5 at latitude 33.5 • S (see Fig. 6).The same high spatial variability 725 given by the small scale structures of the SST and OC images can be appreciated in their corresponding longitudinal profiles displayed in the panel a) and b) of Fig. 6.It is worthy to note the change in the shape of the profiles between the pCO infer 2 and pCO Ctrack 2 and fluxes of CO 2 at large scale, 730 from 5.5 • E to 10.5 • E, showing that the method not only introduces small scale features but also modifies the large scale spatial variability.

Evaluation of using different satellite products
Since the underlying aim of this work is to develop a method-735 ology to infer super-resolution pCO 2 from space using remote observations, we perform a validation study of the different data used in the inferring computations.This provides us an evaluation of which satellite products are more suitable for our methodology and thus a gain in confidence in 740 our method as well as a better understanding of its limitations.The evaluation analysis is addressed taking into account two main concerns: one related to the number of valid points yielded in the pCO infer 2 field, and another with regard to the degradation of the information contained in the tran-745 sition fronts.A valid point is a pixel where we have simultaneously Chl, SST and pCO 2 values from CarbonTracker, from which we can obtain a value of pCO infer 2 , in other words without missing information.One example comparing the reconstructed pCO 2 field obtained from the above-mentioned 750 product combinations is plotted in Fig. 7.The general pattern is quite similar in all of them with some differences in the details of the small scales and in the missing points due to cloudiness (white patches).This example clearly shows how different coverage of the pCO 2 can be in the field depending The same pattern with an area of higher pCO 2 between 24 • S and 30 • S and lower pCO 2 values outside this region is produced with the three combinations.The most noticeable differences are located in the most northern region and in the south-eastern region off Benguela.This can be quantified by computing the standard deviation of the reconstructed pCO 2 values among the different combination of datasets.Fig. 8 d) shows the spatial distribution of the time average over 2006 and 2008 of the standard deviation computed in each pixel among the pCO infer 2 values obtained from the three products combinations.The larger values of the dispersion (not greater than 5µatm) are found in the northern region from 23 • S to the north and and in the southern region, in particular, in the area from 31.5 • S to the south and from 11 • E to the east.The low value of the dispersion indicates that the method is robust when different datasets are used in the inference.
First, we compute the number of valid points in the pCO infer 2 field for each product combination.Table 3 summarizes the total number of valid points for each products combination for both years 2006 and 2008.As expected, the number of valid points is found to be the highest for the combination of merged products OSTIA SST and GLOBCOLOUR-GSM with N GO =27313043 points, N OM =9800776 points.Looking at the different proportions, we find that the number of valid points is 2.78 times larger when using the merged products OSTIA and GLOBCOLOUR-GSM than using OSTIA and MERIS, 1.33 times larger than using MODIS and GLOBCOLOUR-GSM 790 and 1.08 times larger using OSTIA SST and GSM Chla than using MODIS SST and GSM Chl a. Further, if we know that the total number of pixels in the domain taking out the points of the CarbonTracker mask and for the two years is N p =55711378, one can estimate the loss 795 of valid points for each combination, LP x .LP x is computed by dividing the relative difference between the number of total available pixels in the domain N p and the number of points in the inferred pCO 2 field obtained for each product combination,N x , with respect to the total 800 number of pixels N p , LP x = Np−Nx Np 100%.Here the subscript x refers to the product combination (e.g.LP x = LP OM , LP OG and LP M G for the loss of valid points with the OSTIA-MERIS, OSTIA-GLOBCOLOUR and MODIS-GLOBCOLOUR products combination, respectively).The 805 loss of valid points due to cloudiness in the ocean colour and SST images is less severe for the OSTIA-GLOBCOLOUR combination with a loss of 51%, and is more affected by the cloudiness the OSTIA-MERIS combination with a loss of 82%.

810
Next we explore the quality of the information contained in the transition fronts, in particular, in the non-merged products such as MERIS OC and MODIS SST as compared to the merged products: GLOBCOLOUR OC and OSTIA SST.The PDFs of pCO 2 values from CarbonTracker and .Indeed the histograms show also a better agreement between merged products and CarbonTracker: the peak of the PDF for pCO infer 2 is closer to CarbonTracker peak in the case of OSTIA and GLOBCOLOUR than when using MERIS and MODIS products.Furthermore, to examine the transition fronts for the different products we compute the singularity spectra for the three product combinations (see Figure 10).At low values of h (singularity exponent), related to the most singular manifolds, the shape of singularity spectrum for inferred data from merged products better matches a binomial cascade, with an improved description of the dimension of the sharpest transition fronts.We know from the theory, that tracers advected by the flow in the turbulent regime, as happens in the ocean, shows a multifractal behavior with a 835 characteristic singularity spectrum D(h) similar, for some types of turbulence, to D(h) for the binomial multiplicative process.

Validation with in-situ measurements 840
Next, we perform a validation analysis of the results of our algorithm to infer pCO 2 at super-resolution with field observations of oceanic pCO 2 .We perform the validation using pCO 2 ocean data from in-situ measurements (pCO insitu 2 ) taken in the Benguela region (see Section 2.3).We decided 845 to carry out directly the validation on pCO 2 rather than on the air-sea CO 2 flux since the field measurements do provide oceanic pCO 2 data.The number of valid intersections is the largest with the OSTIA-GLOBCOLOUR combination (Table 4).To quantify the loss of valid intersections between the in-situ measurements and points in the pCO inf er 2 field, likely due to cloudiness, we compute the relative difference between the number of measurements into the CarbonTracker domain and the valid points in the inferred pCO 2 field with respect to the number of intersections measurements of each cruise and the pCO Ctrack where N is the number of intersections.at the same intersection, -Mean Relative Error (RE): average over all the inter-910 sections of the errors of the estimated values of pCO 2 (CarbonTracker or inferred) with respect to the refer- We started the statistical validation by analyzing each QUIMA cruise separately (not shown) and we found that in most of the cruises, the absolute error for inferred pCO 2 is relatively small (less than 15 µatm) except on August 21, 2006 andMay 17, 2008 with an error of 44 µatm and 30 µatm, respectively.Then we address the global validation using all available cruises during these years.
We summarize in Table 4 the results of the computations of the errors given by Eq. 7 to Eq.10 by making averages over all valid intersections found during 2005,2006  In this work we have presented a method to infer high resolution CO 2 fluxes by propagating the small scale information given in satellite images across the scales of a multi-resolution analysis determined on the critical transitions given by singularity exponents.More specifically, we have also been obtained in the inferred pCO 2 , showing that the inferring algorithm is capturing the small scales features of the pCO 2 field.The examination of different combina-975 tions of Ocean Colour and Sea Surface Temperature (SST) products reveals that using merged products, i.e.GLOB-COLOUR, the quality and the number of valid points in the pCO 2 field are increased.We show that mean absolute errors of the inferred values of pCO 2 with respect to in-situ mea-980 surements are smaller than for CarbonTracker.The statistical comparison of inferred and CarbonTracker pCO 2 values with in-situ data shows the potential of our method as well as the shortcomings of using CarbonTracker data for the estimation of air-sea CO 2 fluxes.These results indicate that the outputs 985 of our algorithm will be only as good as the inputs.
We are aware that further investigations can be performed in order to improve the algorithm.On one hand the multiple linear regression coefficients could be derived differentiating the seasons (i.e.coefficients would vary as a function of cal-990 endar month) considering the marked seasonal cycle in the Benguela upwelling system.Additionally, future work will be focused on the extension of the computations to larger areas in order to infer global high resolution CO 2 fluxes.This will allow more comprehensive and robust validation from 995 more in-situ measurements.

Fig. 1 .
Fig. 1.Estimated fluxes from CarbonTracker data.Shown are the results on the Benguela upwelling system on March 23, 2006.Left are the CarbonTracker fluxes, right are our results.

375
striction we only document the offshore conditions of this upwelling system.Owing to the relatively large number of cruises during 2005, 2006 and 2008 (a total of 20 cruises, representing 83% of all available cruise data from 2000 through 2010), in this validation, we focus the analysis on 380 the set of QUIMA-cruises during 2005 (QUIMA2005), 2006 (QUIMA2006) and 2008 (QUIMA2008) and we present the global analysis using all available cruises during these three years.Santana-Casiano et al. (2009) analyzed this data to study the sea surface pCO 2 , fCO 2 and CO 2 air-sea fluxes 385 offshore of the Benguela upwelling system between 2005 and 2006 (for each month from July 2005 up to November 2006) and González-Dávila et al. (2009) extended the study including cruises data from 2007 to 2008.The QUIMA line crosses the region between 5 • S and 35 • S, with all the cruises 390 following the same track.

Fig. 2 .
Fig. 2. Snapshot of Chl-a fields corresponding to September 21, 2006, regridded at 1/32 • of spatial resolution from MERIS (a) and GSM GLOBCOLOUR (b).c) and d) are the spatial distribution of Singularity Exponents of the Chl-a plotted in a) and b) respectively

Fig. 3 .
Fig. 3. Snapshot of SST fields corresponding to September 21, 2006 regridded at 1/32 • of spatial resolution from OSTIA (a) and MODIS (b).c) and d) are the spatial distribution of Singularity Exponents of the SST plotted in a) and b) respectively.
-scale inference of pCO 2 data Among the functionals that are most commonly used to analyze the scaling properties of multifractal systems, wavelets occupy a prominent position.Wavelets projections are integral transforms that separate the relevant details of a signal at 520 different scale levels, and since they are scale-tunable, they are appropriate to analyze the multiscale behavior of cascade processes and to represent them.However, as shown inPottier et al. (2008),Yahia et al. (2010) andPont et al. (2011a) not all multiresolution analyses are equivalent but the most 525 Fig. 4. a) Probability distribution functions (PDF) of Chl-a values derived from the three products: MERIS, GLOBCOLOUR-AVW and GLOBCOLOUR-GSM.(b) PDF of SST values for OSTIA and MODIS products.c) PDFs for the singularity exponents of Chl for the different Ocean Colour products.d) PDFs for the singularity exponents of Chl for the different SST products.e) Singularity spectra corresponding to c). f) Singularity spectra corresponding to d).
of super-resolution pCO 2 and air-sea fluxes of CO 2 offshore of the Benguela upwelling system

Fig. 5 .Fig. 6 .
Fig. 5. Maps of a) SST from OSTIA at 1/32 • of spatial resolution, b) Chl at 1/32 • of spatial resolution from GSM GLOBCOLOUR products, c) ocean pCO2 from CarbonTracker at the spatial resolution of 1 • , d) inferred pCO2 at super-resolution (1/32 • ) derived from OSTIA SST and GLOBCOLOUR-GSM Chl-a shown in a) and b) respectively, e) Air-sea CO2 flux as derived from CarbonTracker and f) Air-sea CO2 flux computed from super-resolution pCO2 shown in d) at 1/32 • .All images correspond to March 22, 2006.White color corresponds to invalid pixels due to cloudiness and points inside of the CarbonTracker land mask Fig. 7. a) Map of pCO2 field at low resolution from CarbonTracker.Reconstructed pCO2 field at super-resolution using: b) OSTIA SST and MERIS Chl-a, c) OSTIA SST and GSM-GLOBCOLOUR Chl-a and d) MODIS SST and GSM-GLOBCOLOUR Chl-a.All maps correspond to September 21 2006.

Fig. 8 ..Fig. 9 .
Fig. 8. Spatial distribution of the pCO infer 2 values to the time averages over both 2006 and 2008 years using: a) OSTIA SST and MERIS Chla, b) OSTIA SST and GSM-GLOBCOLOUR Chl-a and c) MODIS SST and GSM-GLOBCOLOUR Chl-a.d) Map with spatial distribution of the standard deviation for the pCO infer 2 among the different combination of the datasets.

2-
field, L infer = N Ctrack − N infer N Ctrack 100%.We repeat such a computation for the three product combinations.The percentage of losses of intersections in inferred field L infer becomes twice as large than in the case of the OSTIA-SST and MERIS-Chl combination, and even higher than with the CarbonTracker domain mask.In order to quantitatively study the difference between values of pCO Ctrack 2 Mean Difference (MD): average over all the intersections of the difference between pCO Ctrack 2 2 Fig. 10.a) Empirical PDFs for the singularity exponents of pCO2 fields from CarbonTracker and from the cascade of the three product combinations.b) Associated singularity spectra.In these computations we use all the pCO2 values obtained in 2006 and 2008.

965
have reconstructed maps of CO 2 fluxes at high resolution (4 km) offshore of the Benguela region using SST and ocean colour data at this resolution, and CarbonTracker CO 2 flux data at low resolution (110 km).The inferred representation of ocean surface pCO 2 improves the description provided by 970 CarbonTracker, enhancing the small scale variability.Spatial fluctuations observed in latitudinal profiles of in-situ pCO 2 (Surface Ocean CO 2 Atlas, http://www.socat.info)over the 2000-2010 period in our region of interest we find the following cruises with pCO 2 measurements:

Table 1 .
Values of the standard deviation, skewness and kurtosis for the different products.

Table 2 .
Values of the standard deviation, skewness and kurtosis of the singularity exponents for the different products.
660super-resolution from pCO 2 at low resolution derived from

Table 3 .
Number of valid points in the pCO2 fields and their difference between the three combinations of MERIS or GLOBCOLOUR CHL with OSTIA or MODIS SST in the area of Benguela.followed by the combination MODIS SST and GLOB-COLOUR Chl with N M G =20397047 points and finally by the OSTIA SST and MERIS Chl combination with

Table 4 .
and 2008.The absolute error, AE is smaller in the three cases of pCO infer 2 (17.77, 16.47 and 16.62 µatm for OSTIA-MERIS, OSTIA-GLOBCOLOUR and MODIS-GLOBCOLOUR combinations, respectively) than for pCO Ctrack Mean error, absolute error and relative error of pCO2 values obtained from CarbonTracker and pCO2 values inferred at super-resolution with respect to values of pCO2 measurements during the QUIMA2005/QUIMA2006/QUIMA2008 cruises in the Benguela region.22.08 and 22.07 µatm, respectively), showing that the 930 estimated pCO 2 field at super-resolution using our algorithm is improving the pCO 2 field obtained from CarbonTracker.The smallest AE is for the combination of SST and Chl provided by merged products.The values of pCO Ctrack values compensate each other (M D infer = 0.15, 3.42 and 8.42 µatm).In all cases the M D Ctrack and M D infer are positive, meaning that the pCO 2 values are overestimated.Finally, comparing the relative error of and with the in-situ measurements (see Table 5), we obtain 458 mutual intersections.We obtain similar results when taking into account all the intersections.The absolute error is smaller 950 in the case of pCO infer 2 , 17.65 µatm, than with pCO Ctrack 2 , 20.24 µatm, indicating that our algorithm is improving the estimation of ocean pCO 2 .The smallest AE is again for the combination with merged products.M D is positive showing that the most of the time pCO infer It can be appreciated in Figure11).Again the relative error is small, less than 0.06, for all the product combinations.

Table 5 .
Mean error, absolute error and relative error of pCO2 values obtained from CarbonTracker and pCO2 values inferred at super-resolution with respect to values of pCO2 measurements during the QUIMA2005/QUIMA2006/QUIMA2008 cruises in the Benguela region at the same intersections.