References

Responses to comments from Dr. Zhang Question (Q) Your paper is very interesting and well written. However, I have several concerns to discuss with you Page 11323, section 2.3 is confusing. It is described that an algorithm for estimating landscape and regional C fluxes including following four steps. Steps one and two are carried out in your paper. I am not clear which parameters are optimized in this study. Which method used to conduct parameter optimiza-

Upscaling of gross ecosystem production to the landscape scale using multi-temporal Landsat images, eddy covariance measurements and a footprint model B. Chen Growing interest in climate change has stimulated recent research that aims to quantify components of the natural carbon (C) cycle.The eddy-covariance technique (EC) is commonly used to directly measure the CO 2 , water vapor and energy exchange between the atmosphere and terrestrial ecosystems (Baldocchi, 2008).Today, there exist more than 400 EC-flux towers across continents.EC measurements are a rich source of information on temporal variability and environmental controls of CO 2 exchange between the atmosphere and terrestrial ecosystems (Law et al., 2000).These global EC datasets provide investigators opportunities and information to 1) explore emergentscale properties by quantifying how the metabolism of complex ecosystems respond to perturbations in climate variables on diurnal, seasonal, interannual and decadal time scales and elucidate physical and biological controlling factors (Law et al., 2000;Baldocchi, 2008); 2) examine carry-over effects that may be introduced by either favorable or deleterious conditions during antecedent years (Barford et al., 2001); 3) observe a disturbance and the recovery from it or to span a natural sequence of ecological development coupled with fluctuations in climate (Amiro et al., 2006;Stoy et al., 2006); and 4) test and validate ecosystem process models (Chen et al., 2007;Urbanski et al., 2007), since most of these models span timescales from hours to decades.Although the available EC data have been rapidly accumulating, much of this information is of limited use because of the difficulties/uncertainties in i) assessing/interpreting the associated measuring biases of EC data and ii) upscaling of the EC fluxes at the ecosystem (typically less than 1-3 km 2 for each site) to larger scales, e.g.landscape and regional scales.
The EC method is based on measurements of turbulent fluctuations of the vertical velocity and the concentration of a passive tracer.Knowledge of the area's soil and vegetation that impacts the EC flux is clearly important both in planning the site tower location and in the interpretation of measured fluxes (Finnigan, 2004).The adoption of the EC technique to estimate surface exchange is based on the assumption that certain meteorological conditions (e.g.horizontal homogeneity, steady-state, and nonadvection) are satisfied (G öckede et al., 2004).Since such conditions are often violated in complex terrain, e.g. at flux monitoring sites in forests, correct interpretation of EC-data is still a matter of some difficulty (Sogachev et al., 2004).In particular, the spatial variability of vegetation density influences the lower atmospheric circulation and surface exchange of energy, water and C over a wide range of scales (e.g., Shen and Leclerc, 1995;Buermann et al., 2001;Cosh and Brutsaert, 2003).As a result, evaluation of the spatial representativeness of long-term accumulated EC-flux measurements is still challenging (Chen et al., 2009a).
Footprint analysis is a recognized part of the establishment and siting of flux towers and the analysis of their output (Finnigan, 2004) measurements over a heterogeneous surface depends largely on the footprint over which the fluxes are sampled.The temporal and spatial variability of footprints and the associated influence of varying site heterogeneities on tower flux measurements has yet to be fully been investigated, although this information is critical needed for interpretation and for making a more wide use of the globally available EC datasets which have been rapidly accumulating.
In order to use EC flux measurements to provide estimates of components of the natural C cycles at landscape ( 101 -10 2 km 2 ), regional (10 3 -10 6 km 2 ), or hemispheric to global land surface (10 7 -10 8 km 2 ) scales, they must be reasonably "upscaled" using either models and/or remote sensing measurements (e.g.Earth observation (EO) data).However, it has been proved that, it is an extremely challenging task to scale up those EC measurements from stand-level to a region or global scales because of the large spatial heterogeneity and temporal dynamics of ecosystems across complex landscapes and regions and the nonlinearity inherent in ecophysiological processes (Levy et al., 1999;Chen et al., 2007;Hilker et al., 2008).
Satellite remote sensing can provide consistent and systematic observations of vegetation and ecosystems over large spatial extents on variable spatial and temporal resolutions.For example, the Moderate Resolution Imaging Spectroradiometer (MODIS) is a 36 band spectrometer providing a global data set every 1-2 d.The spatial resolution of MODIS (pixel size at nadir) is 250 m for channel 1 and 2 (0.6-0.9 µm), 500 m for channel 3 to 7 (0.4-2.1 µm) and 1000 m for channel 8 to 36 (0.4-14.4 µm), respectively.The Landsat Enhanced Thematic Mapper Plus (ETM+) is a sensor carried onboard the Landsat 7 satellite and has acquired images of the Earth at 30 m spatial resolution with a 16-d repeat cycle.Data from the satellite-borne MODIS are currently used in the calculation of global weekly gross primary productivity (GPP) at 1-km spatial resolution (Running et al., 2004;Coops et al., 2007).These data with variable spatial and temporal resolutions, however, require appropriate methods for upscaling and interfacing EC-measurements to satellite observations (Drolet et al., 2008;Hall et al., 2008).It is challenging to compare the estimated GPP using MODIS satellite Figures

Back Close Full Screen / Esc
Printer-friendly Version Interactive Discussion data at 1-km resolution with the EC-derived GPP because of mismatch between them in spatial scales (Xiao et al., 2004).Whereas, the Landsat imagines have fine resolutions (i.e. 30 m), which are ideal for spatial scaling and comparison of C fluxes derived from satellite-borne and EC data with assistance of flux footprint analysis (Chen et al., 2009a).
The objective of this study is twofold: 1) to examine biophysical performance of vegetation indices in relation to seasonal dynamics of CO 2 fluxes; and 2) to develop a practical approach for upscaling GPP to the landscape scale and assessing the EC sensor location biases (i.e. the influence of patch-scale heterogeneities and the spatial representativeness of the EC flux footprint on the uncertainty in EC CO 2 flux data) using the Landsat images and a recently developed footprint model (Chen et al., 2009a).Qian Yanzhou (QYZ) EC flux tower site, locating in South China subtropical monsoon climatic zone at 26 • 44 48 N, 115 • 04 13 E, which belongs to the Chinaflux network, was selected as an experiment site.The EC data measured at this tower and the Landsat images in 2004 were used in this study.Firstly, footprint climatology across multi-temporal scales (i.e.daily, biweekly, monthly and annual) was calculated using the Simple Analytical Footprint model on Eulerian coordinates (SAFE, Chen et al., 2008Chen et al., , 2009a)); secondly, biweekly high-spatial-resolution GPP maps were produced based on the Landsat data using a light-use efficiency modeling approach (e.g., Xiao et al., 2004); and finally these year-round estimates of GPP were compared to that directly derived from EC measurements using two alternative weighting approaches, i.e. footprint weighting or "equal" weighting.These comparisons provided opportunities to optimize the parameters used in the light-use efficiency model and to estimate the tower sensor location biases caused by the variations in size and orientation of footprint climatology and patch-scale heterogeneities. 2 Materials and methods

Site description
The study site, QYZ Experimental Station, which belongs to Chinese Ecosystem Research Network (CERN) and ChinaFLUX network, is located in southeastern China (elevation 102 m).The mean annual air temperature was 17.9 • C, and mean annual precipitation was 1485 mm (1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004).Most areas in the same latitude zone as QYZ around the world are arid steppes or deserts.The warm and humid environment in QYZ is the result of unique southeast monsoon.QYZ is located on gently undulating terrain with slopes between 2.88 and 13.58 (Wen et al., 2006).The QYZ ecological experimental station was established in 1983, at that time the area was mainly covered by wild grasslands, shrub lands and some sparse Pinus massoniana.After the experimental station setup, the land covers of QYZ were changed much by scientists (Huang et al., 2007) 1.According to the field measurements in August of 2003, the leaf area index (LAI) of the plantation was 4.5 (Huang et al., 2007).The understory shrub mainly includes Loropetalum Chinense and Lyonia compta.The soil is red soil, which weathered from red sand rock.to equal NEE in well-mixed conditions, i.e., friction velocity (u * ) larger than a threshold value (u th * =0.2 m s −1 ).Daytime R eco and values to fill nighttime gaps were calculated using the relationship between nighttime NEE in well-mixed conditions and soil temperature at the 5-cm depth (Wen et al., 2006).Net ecosystem productivity (NEP) was calculated as NEP=−NEE.GPP was calculated as the sum of daytime NEP and calculated daytime R eco .

Data-model fusion and upscaling framework based on Landsat images
An algorithm for estimating landscape and regional C fluxes (e.g.GPP) includes the following four steps (Fig. 1): i) A satellite-based vegetation photosynthesis model was adopted to produce GPP maps at 30 m resolutions using Landsat and climate data; ii) EC flux footprints for the corresponding Landsat repeating periods (normally biweekly) were calculated using a recently developed footprint model (SAFE, Chen et al., 2009a); iii) By assuming the footprint integration of remotely-sensed GPP to be comparable with the EC-derived GPP values, several key parameters in the satellite-based vegetation photosynthesis model were optimized using the data-model fusion technique; and iv) The updated satellite-based vegetation photosynthesis model was used for data fusion with other satellite data (e.g.MODIS) or directly used for estimating landscape/regional GPP.

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion (Running et al., 2004) and the vegetation photosynthesis model (Xiao et al., 2004), the algorithm (Fig. 2) used in this study relies on the light-use efficiency (ε) approach relating GPP to the amount of absorbed photosynthetically active radiation (APAR) (Monteith, 1966(Monteith, , 1972) ) such that, where PAR is the photosynetically active radiation (in µmol photosynthetic photon flux density, PPFD), f PAR chl is the fraction of PAR absorbed by leaf chlorophyll in the canopy, and ε is the light use efficiency (µmol CO 2 /µmol PPFD).Light use efficiency (ε) is affected by leaf phenology, temperature, and water: where ε 0 is the apparent quantum yield or maximum light use efficiency (µmol CO 2 / µmol PPFD) for a given land cover type or vegetation function type, and P m , W m and T m are the modifiers for the effects of leaf phenology, water and temperature on light use efficiency of vegetation, respectively.

Model parameter estimates
Different parameters and inputs for the satellite-based algorithm are estimated in different ways: i) the fraction of PAR absorbed by leaf chlorophyll in the canopy (f PAR chl ) and the modifiers (P m , W m ) are estimated using Landsat imagery data; ii) PAR and temperature modifier (T m ) are calculated using climate data (either from tower measurements or climate models); and iii) the maximum light use efficiency (ε 0 ) is referred to the land-cover-related look-up table and then modified/optimized using EC tower C measurements and footprint climatology.
To accurately estimate f PAR chl in forests is a challenge to both radiative transfer modeling and field measurements.Significant effort and progress have been made in developing advanced vegetation indices that are optimized for retrieval of f PAR from individual optical sensors (Gobron et al., 1999;Govaerts et al., 1999).In this study, the Introduction

Conclusions References
Tables Figures

Back Close
Full f PAR within the photosynthetically active period of vegetation is estimated as a linear function of the Enhanced Vegetation Index (EVI), where EVI directly normalizes the reflectance in the red band as a function of the reflectance in the blue band (Huete et al., 1997): where G=2.5, C 1 =6, C 2 =7.5, and L=1 (Huete et al., 1997); and ρ nir , ρ red and ρ blue are the reflectance of near infrared bands, red bands and blue bands, respectively.The leaf phenology modifier (P m ) is estimated using the Normalized Difference Vegetation Index (NDVI) and the Land Surface Water Index (LSWI).NDVI (Tucker, 1979;Field et al., 1995) and LSWI (Xiao et al., 2002) are calculated, respectively using Eqs.( 5) and ( 6): NDVI = (ρ nir − ρ red )/(ρ nir + ρ red ), (5) LSWI = (ρ nir − ρ swir )/(ρ nir + ρ swir ), where ρ nir , ρ red and ρ swir are the reflectance of near infrared bands, red bands and short infrared bands, respectively.P m is calculated at two different phases, depending upon life expectancy of leaves (deciduous versus evergreen): During bud burst to leaf full expansion 1 After leaf full expansion (7) The timings of bud burst and leaf full expansion can be identified using NDVI.The effect of water on plant photosynthesis (W m ) has been estimated as a function of available soil content in plant root zone and water vapor pressure deficit (VPD) in a number of process-based ecosystem models (e.g., Chen et al., 2007) and remote-sensing based models (e.g., Running et al., 2000).and canopy, and VPD represents evaporative demand in the atmosphere.Leaf and canopy water content is largely determined by the dynamics of both soil moisture and VPD.As the first order of approximation, here following the alternative and simple approach that uses a satellite-derived water index (Xiao et al., 2004), we use Eq. ( 8) to estimate the seasonal dynamics of W m : where LSWI max is the maximum LSWI within the plant growing season for individual pixels.The temperature modifier T m was estimated at each time step, using the equation developed for the terrestrial ecosystem model (Raich et al., 1991): where T min , T max and T opt are the minimum, maximum and optimal temperature for photosynthetic activities, respectively.Their values are respectively set to be 0, 35 and 20 • C in this study.If air temperature falls below T min , T m is set to be zero.
The ε 0 values vary with vegetation types, and the information about ε 0 for individual vegetation types can be obtained from a survey of the literature (Ruimy et al., 1995) and optimized using EC tower measurements.According to the work (Zhang et al., 2006), the ε 0 value was estimated to be 0.032 µmol CO 2 /µmol PPFD in this study stand in 2004.

Scaling of remotely sensed GPP to EC-flux derived GPP
The ecosystem-scale overall GPP (F GPP,rs ) was up-scaled from the spatial distributed GPP field (F GPP,rs ) using Eq. ( 10), F GPP,rs (x,y)φ pure (x,y)dxd y . (10) In Eq. ( 10), φ pure is the pure footprint and Ω Π is the upwind footprint source area.Both F GPP,rs and φ pure were estimated at 30-m resolution.The EC-derived GPP value is expected to be close to the value of F GPP,rs weighted with φ pure over Ω Π ."Equal" integration of F GPP,rs was also calculated using Eq. ( 10) by setting φ pure to be 1 over Ω Π .Both footprint integration and "Equal" integration of F GPP,rs were calculated at biweekly time steps first and then summed up to gain annual values for comparisons.
Comparisons of ecosystem-scale GPP estimates between integrated from F GPP,rs and directly derived from EC measurements were made for the year of 2004 at the QYZ site.

Sensor location bias
The spatial representativeness of the footprint is given by the sensor location bias (∆) following Schmid (1997), where F GPP,φ and F GPP,φ=1 is the footprint-weighted and "equally"-weighted GPP, respectively.Alternatively, the sensor location bias was also estimated by replacing the "equally" integrated F GPP,φ=1 with the GPP values of the tower location pixel.The values of ∆ and of root bias δ (δ = √ ∆) were calculated at biweekly time steps and then averaged to gain monthly and annual values for the year of 2004 at QYZ.

Dataset used
The Landsat image data in 2004 for the area (domain) of 6 km×6 km centered at the QYZ tower acquired from the website of US Geological Survey at http://glovis.usgs.gov/.The NDVI (Fig. 3), EVI and LSWI were calculated based on the Landsat original data after atmospheric and other corrections.The QYZ tower data in 2004 were used for footprint calculations and several key parameters (e.g.ε 0 ) of the satellite-based GPP algorithm estimates and for GPP comparisons.The seasonal dynamics of C flux components measured at QYZ in 2004 were shown in Fig. 4.

Variations in seasonal footprint climatology
Figure 5 shows the mean monthly daytime pure footprints and the corresponding cumulative footprint contours for every other months in 2004 for QYZ.The seasonal variations in size and orientation of footprint for QYZ are significant.The areas of 90% cumulative footprints in winter months were less than 1 km 2 whist were as large as 2.5 km 2 in summer months.

Comparing estimated values of remotely sensed and EC derived GPP
The remotely-sensed biweekly mean GPP maps at 30 m resolution were produced based on Landsat image data, which overlaid by the corresponding period daytime 11328 Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion footprint.Figure 6 shows an example at annual time step.The spatial variations of remotely-sensed GPP were not obvious for this plantation stand.The elliptical shape of the annual footprint distributed along the NNW-SSE prevailing wind directions.The size of the annual mean 90% cumulative footprint was about 1.5 km 2 .The ecosystem-scale overall biweekly GPP (namely footprint-integrated GPP) was up-scaled from these spatially distributed GPP field using Eq. ( 10).The "equally" integrated GPP was also calculated using Eq. ( 10) by setting φ pure =1 for comparisons.The remotely-sensed GPP using Landsat data reasonably followed the seasonal dynamics of observed GPP (Figs. 7 and 8).During winter time (when GPP<=∼3 g m −2 d −1 ), the remote-sensing algorithm underestimated the measured GPP; while during summer time (when GPP>=∼7 g m −2 d −1 ), the remote-sensing algorithm overestimated the observed GPP (Figs. 7 and 8).The footprint integrated GPP values were closer to ECderived GPP values than the "equally" integrated GPP and the tower pixel's GPP values though their differences were small.The annual mean GPP values were 5.15, 5.06, 4.51 and 5.09 g C m −2 s −1 for EC-derived, footprint-integrated, "equally"-integrated and the tower pixel's GPP, respectively.As shown in Figs. 3 and 6, the land surface heterogeneity is quite small for this plantation stand, therefore the annual mean tower location biases was less than 5%.Even for this "ideal" EC tower location, the estimated biases in annual sum of GPP at the landscape scale between with or without footprint consideration could be as large as 200 g C m −2 .

Discussion
The seasonal phenologically dynamics of canopy development (leaf flush, expansion, senescence, fall), in relation to their biophysical, biochemical (e.g., chlorophyll and other pigments, nitrogen) and optical properties, which in turn influence both biophysical parameters (e.g., albedo, latent and sensible heat fluxes) and biogeochemical parameters (e.g., photosynthesis) of the land surface (Xiao et al., 2004;Li et al., 2007).Therefore, the time series of the vegetation indices have potential to provide valuable insight into the processes (e.g., growing season length and water condition) that regulate ecosystem carbon exchange.NDVI and EVI are widely applied to detect the information on leaf area index (Chen et al., 2006) and f PAR (Running et al., 2004;Xiao et al., 2004).A number of studies demonstrated that the seasonal dynamics of MODIS EVI agreed better with GPP than MODIS NDVI for different ecosystems (e.g.evergreen needleleaf forest, Xiao et al., 2004; alpine ecosystems in Qinghai-Tibetan Plateau, Li et al., 2007).The availability of time-series data of SWIR and NIR bands from the new generation of optical sensors (e.g., VGT, MODIS, Landsat 7) offer new opportunity for quantifying canopy water content at large spatial scales (Xiao et al., 2004).These sensor-specific advanced vegetation indices (e.g.EVI and LSWI) have been optimized for the Moderate Resolution Imaging Spectroradiometer (MODIS), the Global Imager (GLI) and the VEGETATION sensors (Li et al., 2007).Clearly, there is a need to examine those advanced vegetation indices for all available satellite sensors in relation to leaf phenology and the seasonal dynamics of GPP across the flux tower sites in various biomes.Limited number of studies had evaluated radiometric and biophysical performance of vegetation indices (e.g.EVI, NDVI) from Landsat data, probably because of mismatch of spatial resolutions between Landsat data and EC flux towers.
The performance of a new algorithm based on the advanced vegetation indices using Landsat data was evaluated with assistance of a footprint model.A good agreement between predicted and EC measured biweekly GPP in 2004 in an evergreen needleleaf forest at QYZ (R 2 =0.9, p<0.001,Fig. 8) indicates that there exist a good quantitative relationship between the Landsat vegetation indices and CO 2 flux data, in terms of the seasonal phase and magnitude of photosynthesis.However, the discrepancies between remotely-sensed and EC-derived GPP are still large, especially in winter and in the middle of summer (Fig. 7).Those large discrepancies may be attributed to three sources of errors.The first source is the sensitivity of the remote sensing algorithm to PAR and T m .The parameters for estimating T m in Eq. ( 9 phases.In this study, we simply assume those parameters have no seasonal variations.This may lead to ill parameterization of T m , such as over-corrected (smaller values of T m ) during the low-temperature periods while under-corrected (larger values of T m ) during the high-temperature periods.The second source is the time-series data of vegetation indices from Landsat satellite images.We used the biweekly Landsat data that have no BRDF correction or normalization, and thus, the effect of angular geometry on surface reflectance and vegetation indices remained.EVI is a semi-empirical mathematic transformation of observed reflectance from individual spectral bands (blue, red and NIR) of optical sensors (Huete et al., 2002).The parameters in Eq. ( 4) may vary seasonally.For simplification, however, we only optimized one value for each parameter and applied it to the whole year period.The third source is the error in EC-derived GPP.The EC measurements themselves are not free from error.The partitioning of NEE into its component fluxes (i.e.GPP and R eco ) has large uncertainties (Falge et al., 2001;Chen et al., 2009b).The two major steps to derive GPP are the gap filling of NEE and estimation of daytime ecosystem respiration, and both of them carried on a lot of uncertainties.The maximum light use efficiency (ε 0 ) is the basis for the remote-sensing based algorithms or models and the accurate estimating of ε 0 is one of the key steps for using the satellite data to estimate either GPP or NPP (Running et al., 1999).In nature, ε 0 is determined by many biological and biophysical factors and soil nutrient conditions.Much attention should be given to the variability of ε 0 among vegetation types across a heterogeneous landscape (Li et al., 2007).Since ε 0 differs significantly among vegetation types, these differences should be accounted for when estimating GPP using remotely-sensed data.These global EC datasets provide investigators opportunities to estimate the ecosystem-scale photosynthetic (including ε 0 ) and respiratory parameters.A widely-used method for those parameters estimation is the Michaelis-Menten approach (Falge et al., 2001).Caution should be paid when we apply the EC-derived ε 0 to remote-sensing-based algorithms for GPP estimates because: i) there are large uncertainties in the EC-derived ε 0 ; ii) the values of EC-derived ε 0 vary seasonally and  14, 3536-3550, 2001. Chen, B., Chen, J. M., Mo, G., Yuen, C-W., Margolis, H., Higuchi, K., and     Remotely-sensed (footprint integration) Remotely-sensed ("equal" integration) Remotely-sensed (tower pixel) Linear (Remotelysensed (footprint integration)) 1:1

CO 2 ,
water vapour and energy exchange, and meteorological variables were measured continuously at the QYZ site from 2003.Measurements, instruments, calculation procedures and gap filling methodologies were described inWen et al. (2006).Net ecosystem exchange (NEE) was calculated as the sum of the EC CO 2 flux above the canopy and the change in CO 2 storage in the air column between the EC-sensor height and the ground(Wen et al., 2006).Ecosystem respiration (R eco ) at night was assumed Introduction Soil moisture represents water supply to the leaves 11325 We use the Simple Analytical Footprint model on Eulerian coordinates(SAFE, Chen     et al., 2008, 2009a) to calculate QYZ EC tower's footprints.The flux footprints were calculated at a grid size of 30 m×30 m (consistent with the Landsat spatial resolutions) covering the area (domain) centred on the towers of 6 km×6 km.The model was run at half-hourly time steps.The half-hourly footprint f (x,y) was rotated along the wind Introduction to biweekly, monthly and annual time steps to yield the information on footprint climatology, φ(x,y).The total footprint Φ = Ω Π φdxd y within the model domain (Ω Π ) equals 1.The calculated footprint provides a map of the contribution for the area around the tower to the integral EC-measured flux component.The detailed description of footprint model and the footprint climatology calculations can been found inChen et al. (2009a).
), T min , T max and T opt are vegetation-type and climate-zone dependent, and may vary with different seasonal

Fig. 8 .
Fig. 8.A comparison of biweekly gross primary production (GPP) values between the observed (EC-derived) and predicted (remotely sensed) in 2004 in an evergreen needleleaf forest at QYZ, China.

. Ge 1 , D. Fu 1 , G. Liu 1 , G. Yu 1 , X. Sun 1 , S. Wang 1 , and H. Wang 1 Introduction
1,2 , Q . The interpretation of EC flux Figures . The EC flux tower was established in late August of 2002.The dominated trees in the flux footprint area are Pinus elliottii, Pinus massoniana and Cunninghamia lanceolata.The stand characteristics at QYZ based on a survey made in 2005 (Wen et al., 2006) are list in Table iii) the EC-derived ε 0 represents the EC flux footprint area, whose sizes and orientations vary with time, in other words, the spatial representativeness of the EC-derived ε 0 is normally different from the satellite image pixels.Making use of the satellite data with fine resolution (e.g.Landsat data) to estimate GPP, flux footprint modeling should be involved to optimize the remote-sensing-based algorithms' parameters, such as ε 0 , T m etc.These optimized parameters then can transfer to the applications of other satellite images with coarse resolution (e.g.MODIS) by applying data-fusion techniques.Combining vegetation indices (e.g.EVI, LSWI) from different multi-temporal/spatial satellite sensors' data, climate data (PAR, temperature), optimized ε 0 parameter for individual vegetation types, and flux footprint modeling and data-fusion, the remote-sensing-based algorithms is a powerful tool for estimating landscape/regional or global GPP.Introduction