High-resolution land surface fluxes from satellite and reanalysis data ( HOLAPS v 1 . 0 ) : evaluation and uncertainty assessment

Surface water and energy fluxes are essential components of the Earth system. Surface latent heat fluxes provide major energy input to the atmosphere. Despite the importance of these fluxes, state-of-the-art data sets of surface energy and water fluxes largely differ. The present paper introduces a new framework for the estimation of surface energy and water fluxes at the land surface, which allows for temporally and spatially high-resolved flux estimates at the quasi-global scale (50 S, 50 N) (High resOlution Land Atmosphere Parameters from Space – HOLAPS v1.0). The framework makes use of existing long-term satellite and reanalysis data records and ensures internally consistent estimates of the surface radiation and water fluxes. The manuscript introduces the technical details of the developed framework and provides results of a comprehensive sensitivity and evaluation study. Overall the root mean square difference (RMSD) was found to be 51.2 (30.7) W m−2 for hourly (daily) latent heat flux, and 84 (38) W m−2 for sensible heat flux when compared against 48 FLUXNET stations worldwide. The largest uncertainties of latent heat flux and net radiation were found to result from uncertainties in the solar radiation flux obtained from satellite data products.


Introduction
Water and energy fluxes between the land surface and atmosphere are essential components of the Earth system.At the ecosystem scale, the land-atmosphere fluxes have been mainly measured by a network of flux-tower sites within the frame of FLUXNET (Baldocchi et al., 2001;Baldocchi, 2008).However, to generate global data sets of water and energy fluxes, the use of satellite data as well as models has become indispensable.
Different approaches exist to infer land turbulent surface fluxes by either one of the following methods (Kalma et al., 2008;Wang and Dickinson, 2012): (1) simulations by an offline land surface model (Roads and Betts, 2000); (2) empirical statistical models, e.g., obtained by machine learning techniques or neural networks (Jung et al., 2011); (3) surface energy balance models forced either by satellite remote sensing or reanalysis data (Bastiaanssen et al., 1998;Su, 2002); (4) methods based on Penman-Monteith or Priestley-Taylor equations (Fisher et al., 2008;Miralles et al., 2011;Mu et al., 2007;Zhang et al., 2015); and (5) spatial variability methods (Roerink et al., 2000;Peng et al., 2013b;Peng and Loew, 2014).Novel long-term satellite data records as well as increasing computing capacities allow one to generate spatially (< 10 km) and temporally (< 3 h) high-resolved estimates of surface fluxes at the global scale.The currently existing global data sets have spatial resolutions between 0.01 and 2.5 • and are focused on hourly to monthly timescales (Fisher et al., 2008;Miralles et al., 2011Miralles et al., , 2016;;Mu et al., 2007;Vinukollu et al., 2011;Zhang et al., 2010).The multidecadal trends in global terrestrial latent heat flux have also been investigated and analyzed based on these newly generated products (Jung et al., 2011;Mao et al., 2015;Miralles et al., 2014;Zhang et al., 2015;Zhang et al., 2016).For field and continental scale agricultural applications, ALEXI/DisALEXI (Anderson et al., 2007;Norman et al., 2003) already have the ability to provide very high spatial resolution surface fluxes (up to 10 m resolution) with the Published by Copernicus Publications on behalf of the European Geosciences Union.use of thermal observations from a combination of polar and geostationary orbiting satellites (Anderson et al., 2011).
The Global Energy and Water Cycle Experiment (GEWEX) LandFlux initiative aims for the analysis of existing global land surface flux products and the generation of new data sets of land surface fluxes (McCabe et al., 2016).A comparison of existing global latent heat flux data sets from either land surface models, reanalysis or satellite estimates was conducted within the GEWEX LandFlux-EVAL initiative (Mueller et al., 2011;Jiménez et al., 2011) and a synergy data set has been compiled, which provides latent heat fluxes at monthly timescale and a spatial resolution of 1 • (Mueller et al., 2013).
However, large discrepancies remain in the existing data products.The global mean latent heat flux over land was diagnosed as 45 ± 5 W m −2 with a spread as large as 20 W m −2 and substantial regional and seasonal differences (Jiménez et al., 2011).
These discrepancies might be either related to the different methods applied to estimate the surface fluxes as well as due to different ancillary data sets used (Ershadi et al., 2014;Vinukollu et al., 2011).Recently, McCabe et al. (2016) examined the performance of four commonly used methods for the estimation of surface evaporation with FLUXNET tower-based and globally-gridded forcing data.They found that the root mean square difference (RMSD) ranges from 61 to 101 W m −2 for 3-hourly data from 45 FLUXNET towers.As a parallel and complementary effort to the GEWEX LandFlux initiative, the ESA WACMOS-ET project aimed to identify the appropriate methods for the estimation of latent heat flux and maximizing the use of European Earth Observation data sets.The accuracy of the WACMOS-ET results have been validated against a set of FLUXNET sites.Compared to McCabe et al. (2016), a set of different FLUXNET sites and forcing data sets are investigated by Michel et al. (2016).They found accuracies between 40.8 and 88.5 W m −2 for 3-hourly values comparing against data from 24 eddy-covariance towers (Miralles et al., 2016;Michel et al., 2016).Another important finding from both recent projects is that no single algorithm always outperforms any other method.In addition, the existing models do not capture well the early-morning and late-afternoon transitions in the atmospheric boundary layer (Ershadi et al., 2014).In order to develop a more accurate global latent heat flux product, improvement of the parameterization and sensitivity analysis of the model to forcing data set are still needed (McCabe et al., 2016;Michel et al., 2016).
Only a limited numbers of studies provide evaluation of latent as well as sensible heat fluxes.Previous studies estimated sensible heat flux at the regional scale and validated against limited in situ measurements with accuracies ranging from ∼ 10 to ∼ 100 W m −2 (Jia et al., 2003;Marx et al., 2008;Tang et al., 2011b;Wang et al., 2013;Zhuang et al., 2016).
The present paper introduces a novel framework for the generation of global high-resolution land surface fluxes from satellite and reanalysis data sets.The High resOlution Land Atmosphere surface Parameters from Space (HO-LAPS) framework makes use of meteorological drivers coming from globally available satellite and reanalysis data sets and integrates many of the different components developed in previous studies within a single framework.A state-ofthe-art land surface scheme is used for the estimation of the surface energy and water fluxes.HOLAPS allows for internally consistent estimates of the surface radiation and water fluxes at high temporal (< 1 h) and spatial (< 5 km) resolutions.In particular, the shortwave and longwave surface radiation fluxes are consistently estimated, which is often not the case when satellite-based forcing data from different sources are used, as these can differ, e.g., in their cloud coverage or characterization of the atmospheric humidity profile.The different components of the HOLAPS framework are easily exchangeable as they are coupled through well-defined interfaces.This allows, for instance, for the integration of different approaches for the estimation of surface turbulent fluxes while building on the general HOLAPS infrastructure for providing all required forcing data.The required drivers for HOLAPS comprise satellite data at different processing levels and reanalysis data for a limited number of variables.The modular framework allows for integrating different land surface schemes.
The objectives of the present study are mainly 2-fold.First, we introduce and validate the surface fluxes from the novel HOLAPS framework at quasi-global scales (50 • S, 50 • N).Second, we perform a thorough sensitivity analysis of the impact of different forcing data sets on the accuracy of surface heat flux estimates.The latter is motivated by the question: how much uncertainty is introduced when using globally available satellite and reanalysis data as a driver for land surface models compared to local measurements.The HOLAPS results are validated using tower-based eddycovariance measurements for a wide range of ecosystems and climates.
We first briefly introduce the overall HOLAPS concept and framework developed in Sect. 2. The data sets and methods are introduced in Sects.3 and 4, respectively, followed by the summary and conclusions.

Model
The HOLAPS v1.0 framework is used for the estimation of quasi-global surface water and energy fluxes.It is based on a state of the art land surface model and was in particular designed to make use of satellite and reanalysis data as drivers as well as to maximize internal consistency of the different energy and water fluxes.HOLAPS is used for the estimation of quasi-global surface fluxes at high spatial and temporal resolutions.It is based on a radiation module, a planetary boundary layer model, a soil module and a general module for the exchange of energy and moisture at the surface layer.
All framework components are modular and are easily exchangeable.
Figure 1 shows the general surface state and fluxes simulated by HOLAPS and Fig. 2 shows the general interdependency between the different variables as described briefly in the following.A very detailed technical documentation of the entire model formulation is provided in Appendix B.
The all sky surface solar irradiance R g (W m −2 ) is either obtained from remote sensing products or is directly calculated internally by the HOLAPS radiation module using the MAGIC radiative transfer model (Mueller et al., 2009).The algorithm requires information on aerosol properties and surface albedo (α) as well as total column water vapor content (TCW) (kg m −2 ).Aerosol properties are taken from an aerosol climatology (Kinne et al., 2013).Total column water vapor content can be either used from climatologies or reanalysis data.Details on the accuracy of the MAGIC radiative transfer model is provided by Posselt et al. (2012).When radiation data are used as input, the radiation module calculates in addition the cloud coverage, which is further required for the calculation of consistent longwave radiation fluxes.
The land surface scheme is explicitly coupled to a onedimensional (1-D) mixed layer model for the planetary boundary layer (PBL), which is used to calculate the surface downwelling radiation consistently with the surface heat fluxes.As the PBL temperature and height are directly linked to the surface turbulent fluxes, a combination of the surface heat fluxes with a PBL model helps to better constrain the surface heat flux estimates as has been shown in previous studies, e.g., the ALEXI model (Margulis and Entekhabi, 2001;Anderson et al., 2007).However, while it has been shown that such an approach helps to better constrain the surface heat fluxes, it is rarely used in common methods for the estimation of surface heat fluxes.A mixed boundary layer model is used within HOLAPS (Kim and Entekhabi, 1998;Margulis and Entekhabi, 2001;Smeda, 1979), which calculates the boundary layer height and temperature using prog-nostic equations (see Sect.B2.6 in Appendix B), whereas the boundary layer temperature can be nudged towards available air temperature observations.The soil temperature is calculated using a force-restore approach (Ren and Xue, 2004), which gives the surface temperature (T S ) that is required for the calculation of the longwave surface net radiation budget.
The surface water fluxes comprise vegetation interception and soil moisture dynamics as well as evaporation and transpiration processes.The currently implemented land surface scheme calculates the latent heat flux following the Priestley-Taylor formulation.The surface aerodynamic and canopy resistances are estimated as a function of wind speed, air temperature, soil moisture and surface solar radiation flux.Calculated sensible heat flux feeds directly back into the PBL model, which constrains the diurnal evolution of the surface fluxes as discussed earlier.
The present paper will focus exclusively on the validation of HOLAPS v1.0 results using in situ flux-tower measurements as well as the assessment of the sensitivity of HO-LAPS to forcing perturbations, namely, different forcing data sets.An assessment of spatiotemporal dynamics estimated from HOLAPS and cross-comparison against other existing global data sets, e.g., the LandFlux-EVAL data set (Mueller et al., 2013), will be performed in a separate study.All symbols used throughout the manuscript are summarized in Appendix A.

Data
The HOLAPS framework was in particular designed to (a) make use of globally available satellite and reanalysis data and (b) ensure internally consistent flux estimates.The drivers required to force HOLAPS are summarized in Table 1.These consist of satellite remote sensing and reanalysis data sets, which have been thoroughly validated and which are briefly introduced in the following.The data sets have in common: they (a) provide long-term observations of the required driver variables and (b) provide this information at comparably high temporal and spatial resolutions, which is a major prerequisite.Data sets that are based on geostationary satellite measurements are therefore given preference.Static information on land cover and soil properties is required as well.All data need to be regridded to the computational grid, and temporal interpolation to the HOLAPS timescale is required.Details about the employed interpolation techniques are provided in Appendix B8.

FLUXNET data
Measurements of surface turbulent fluxes are obtained from eddy-covariance towers of the FLUXNET network.These measure the exchange of carbon dioxide, water vapor and energy between terrestrial ecosystems and the atmosphere (Baldocchi, 2003).Standard meteorological measurements are collected as at most stations.The most comprehensive compilation of these flux-tower measurements is available from the "La Thuile 2007" database (Papale et al., 2012).A subset of FLUXNET stations was used for the analysis in the present study.Stations were selected where (a) all variables required to run the HOLAPS model (Table 1) were available, (b) the station provided data with limited data gaps (> 80 % coverage).All data are carefully quality checked and available quality flags are applied to ensure the highest quality of the reference data.
The stations used in the present study are depicted in Fig. 3.A major number of stations are located in Europe and North America, and only a few stations are located in other regions.Table C1 lists all stations (N = 48) that fulfilled the above-described criteria and provides detailed information about data availability and relevant references for each station.The total number of measurement years, which is used for the present analysis, is M = 101 years.FLUXNET data are currently distributed under different data policies.For the present study we only use data from stations that provide their data under a "free fair use" license (http://www.fluxdata.org).
Eddy-covariance measurements are subject to uncertainties from various sources.A common problem is that the eddy-covariance measurements typically do not allow one to close the surface energy balance (R N − G − H − LE = 0) (see Table A1 for definition of acronyms throughout the paper).The energy imbalance for eddy-covariance measurements can be as high as 20 to 30 % on average (e.g., Wilson et al., 2002).The reason for this energy balance closure problem is still not fully understood and subject of ongoing research (e.g., Ingwersen et al., 2015).Several approaches have been developed to empirically correct for the energy closure (Foken et al., 2011;Twine et al., 2000;Wilson et al., 2002;Ingwersen et al., 2015).A simple energy balance correction (Bowen ratio method) is applied in this study following the approach as described in Twine et al. (2000).Further uncertainties in the FLUXNET data occur under stable conditions, as the eddy-covariance method requires turbulent conditions (Berbigier et al., 2001).It should be noted that the eddycovariance measurements are less accurate under rainfall conditions.Previous studies have therefore removed measurements during rain events (Ershadi et al., 2014;Michel et al., 2016).As we applied the quality flags available from the FLUXNET data, many rainfall events were masked already.
A sensitivity study was performed to evaluate if additional masking of rainfall events affects the results of the present study, but no deterioration of the HOLAPS performance during rainfall events could be identified.Therefore, we did not explicitly exclude any rainfall data from the analysis.

Large-scale forcing data
In the following we will briefly summarize the different forcing data sets used within the HOLAPS framework.

Radiation data
The surface solar radiation flux (R g ) is either prescribed from existing satellite data products or can be calculated internally within the HOLAPS framework (cf.Appendix B2.2).In both cases a maximum consistency between the shortwave and longwave radiation fluxes is ensured as the same ancillary data (TCW, cloud fractional coverage) are used.This explicit internal consistency of the radiation flux estimates is unique to the HOLAPS framework.
As the surface solar radiation is a major input to the surface energy balance, it is expected that uncertainties in radiation data will also affect the accuracy of the derived water and energy fluxes.Different approaches to estimate R g are therefore analyzed in the present study.The following radiation data sets are used: -FLUXNET: The radiation data measured at each FLUXNET station are used as a reference as these local measurements are expected to provide the most accurate surface solar radiation estimates for the FLUXNET locations.They also capture local changes in R g at high temporal frequencies (e.g.cloud shadowing) and might also be affected by local effects like topographic conditions.
A. Loew et al.: High-resolution land surface fluxes from satellite and reanalysis data -CM SAF-SIS: The EUMETSAT Climate Monitoring Satellite Application Facility (CM-SAF) has specialized in the generation of long-term climate data records from satellites.As part of their suite of radiation data products (www.cmsaf.eu), the CM SAF provides solar incoming surface (SIS) radiation data at hourly timescales and with a spatial resolution of 0.03 • (Posselt et al., 2012;Müller et al., 2015) for all sky conditions.The CM SAF-SIS is based on data from the series of ME-TEOSAT satellites.It therefore provides only a limited area coverage (see Fig. 3).
-GRIDSAT: The Gridded Satellite data set (GRIDSAT) (Knapp et al., 2011) provides a long-term (January 1980 to present) record of top-of-atmosphere (TOA) radiances in the visible and thermal spectral domains.It is based on the International Satellite Cloud Climate Project (ISCCP) (Rossow and Schiffer, 1991;Knapp, 2008) and provides data every 3 h on an equal angular grid with a resolution of 0.07 • .
The TOA radiances in the visible channels are used to estimate a cloud effective albedo (CAL) (Posselt et al., 2012), which is then used subsequently for the calculation of R g and cloud cover fraction (cf.Sect.B2.2).

Precipitation data
Satellite precipitation data sets are produced from satellite only or combined satellite and ground-based measurements at a variety of spatial (0.25 to 2 • ) and temporal (3-hourly to monthly) resolutions at the global scale.Ground-based precipitation estimates, e.g., from ground-based rain radars provide even higher temporal and spatial resolution, but are available only for limited areas.A comprehensive review and inter-comparison of existing satellite-based precipitation products and their application is provided by Kidd et al. (2012) and Kucera et al. (2013) The TRMM Multisatellite Precipitation Analysis (TMPA) product (3B42 v7) is used for the present study (Huffman et al., 2007).It combines microwave sounding and infrared observations and compensates product biases using rain gauge information on monthly timescales.TMPA provides 3-hourly precipitation information at a spatial resolution of 0.25 • .It has been available since 1998 and covers the geographical extent of 50 • N, 50 • S.
The high temporal frequency of the measurements is a major advantage for flux estimates and the main reason why TMPA is currently used within HOLAPS.The spatial extent of TMPA, however, currently limits the application of HO-LAPS to that same extent (±50 • latitude).

Vegetation data
Leaf area index (LAI) data products from the Moderate Resolution Spectroradiometer (MODIS) instruments (Justice et al., 2002) are used in the present study.We use an enhanced product from Beijing Normal University1 (Yuan et al., 2011), which provides enhanced temporal and spatial consistency of the MODIS LAI fields by post-processing the original MOD15A2 products (Myneni et al., 2002).This results in much more consistent LAI fields than in the original product, which contains abrupt changes in the time series.Surface albedo information is obtained from the ESA GlobAlbedo project (Muller et al., 2012;Potts et al., 2013).Both, LAI data and surface albedo are available every 8 days.As both variables are varying slowly in time, they are linearly interpolated to the model time step.However, it needs to be emphasized that the used LAI and albedo products are not necessarily consistent between each other, as they are derived from different instruments and using different inversion techniques.Such a consistency of land surface parameters could only be achieved through joint surface parameter retrieval approaches, e.g., that provided by Pinty et al. (2011), and is part of ongoing research activities, e.g., within the QA4ECV project (http://www.qa4ecv.eu/).

Reanalysis data
A number of additional fields (temperature, wind speed, total column water vapor path, pressure) are required from global reanalysis as these variables are not available from remote sensing data at the required temporal and spatial scales.Therefore, ERA-Interim reanalysis (Dee et al., 2011) fields are used for that purpose, which provide 6-hourly data on a regular global grid with 512 × 256 grid points, which corresponds to a spatial sampling of ∼ 0.7 • .The reanalysis fields are remapped to the flux-tower locations using bilinear interpolation.The scale mismatch between the used reanalysis field data and the local scale HOLAPS simulations might result in additional uncertainty in the simulations and is investigated in the present study.

Land cover data
Global land cover information is available with a spatial resolution of 300 m from the ESA Climate Change Initiative land cover project (Bontemps et al., 2012;Defourny et al., 2014).The land cover information is used for the spatial discretization of land-cover-dependent parameters in HOLAPS, e.g., roughness length or surface resistance parameters.These are summarized in Table B1.
However, for the present study, no global land cover data set is used as the experiments conducted are only performed on the point scale.The land cover type is known for each FLUXNET station and is therefore used in the present study.

Soil data
Information on soil properties is obtained from the Harmonized World Soil Database (HWSD) (FAO, 2012).The HWSD is based on soil mapping units with varying sizes.Thus, no fixed resolution can be given, but the map is gridded with a spatial spacing of 30 arcsec.The information on soil texture (sand, clay content) is used to derive soil hydrological properties using pedo-transfer functions (Cosby et al., 1984;Rawls and Brakensiek, 1985;Lee, 2005).
As the HWSD is a global data set, the local soil properties might differ from those of the used mapping units.Further uncertainties are introduced by the applied pedo-transfer functions to derived soil hydraulic parameters from soil texture information (e.g., Wösten et al., 2001).

Experimental setup
To quantify the accuracy of HOLAPS and the uncertainties related to the usage of different satellite and reanalysis data sets as drivers, we conduct a series of sensitivity experiments.Using the different data sets introduced in Sect.3.2, we aim to investigate the uncertainty introduced by replacing a locally measured forcing with satellite-based drivers.First a control simulation (CTRL) is conducted, which is based exclusively on local measurements from FLUXNET only.This allows one to quantify HOLAPS accuracies without additional uncertainties from the satellite and reanalysis data sets.Thus, the CTRL simulation is considered as the baseline accuracy of the current HOLAPS land surface scheme.For each site multiple years are used for the simulations (see Table 2).Results are then compared against reference measurements from FLUXNET and the accuracy of the simulations is quantified using various skill scores (cf.Sect.4.2).
Further experiments are conducted by replacing individual drivers (e.g., radiation, precipitation) with data from either satellite observations or reanalysis.This allows one to quantify the additional uncertainty introduced by the usage of these particular data products.The different experiment names allow one to identify the variable that was replaced by satellite/reanalysis data (e.g., experiment T a = air temperature was replaced).
However, as the different data sets cover different spatial domains (cf.Fig. 3) we generated subsets of stations representing the following different spatial domains: -Global (G): global coverage uses the maximum number of FLUXNET stations available.
-±50 • (50): as the precipitation data currently used are available only between 50 • S and 50 • N, we use this spatial domain to analyze the sensitivity to changes in the precipitation forcing.
-Meteosat disc (M): the analysis of the impact of satellite surface radiation data sets on HOLAPS results is investigated for the Meteosat spatial domain, as long-term radiation data sets are only available from the CM SAF for Meteosat so far.
-A few FLUXNET stations are located within the Meteosat disc, but within latitudes of 50 While this experimental setup allows one to quantify the impact of different drivers on the HOLAPS results, it does not allow one to explicitly disentangle different components of the overall mismatch between reference data and model results, which are affected by, e.g., model parameterization uncertainties, uncertainties in ancillary data (e.g., soil information) and spatial representativeness of the used reference and forcing data as well as uncertainties in the reference data itself.This could be achieved, e.g., by perturbing the model input parameters and usage of different ancillary data sets.Nevertheless, for the present study we keep the HOLAPS model setup fixed as described in Appendix B.

Analysis
We compare the net radiation and HOLAPS turbulent heat fluxes with the corresponding reference data from FLUXNET at hourly, daily and monthly timescales using standard statistical skill scores.The variance of the difference between the model simulations and FLUXNET data is a function of (a) the uncertainties of the HOLAPS model itself, (b) the sensitivity of the HOLAPS model to uncertainties in the forcing data (including representativeness error) and (c) uncertainties in the FLUXNET reference data.Uncertainties in the FLUXNET measurements might also result from varying temporal and spatial footprints of the fluxtower measurements (Chen et al., 2011).

Statistical metrics
The mean squared difference between in situ observations (x) and model results (y) is given as The RMSD is defined as the square root of Eq. (1).For the calculation of the centered root mean square difference (cRMSD), the bias is removed in advance.It is then defined as whereas the overbar indicates temporal averaging.This is also related to the Pearson correlation coefficient (r) (Taylor, 2001).The above-defined metrics (r, cRMSD, RMSD) are calculated for each FLUXNET station over the entire analysis period.We then normalize each metric by the corresponding metric obtained from the control experiment to obtain relative deviations of the error skill scores of an experiment and the same score from the CTRL simulation for the same station.

Temporal aggregation and data gaps
The comparison between FLUXNET and HOLAPS is performed on hourly, daily and monthly timescales and the above metrics are calculated for these different aggregation periods.
As the FLUXNET measurements also contain data gaps these might introduce sampling biases.A traceable approach is therefore required to derive the temporally aggregated reference.A daily mean is therefore only calculated if at least 16 h (i.e., two-thirds) of valid data were available from the FLUXNET measurements on that particular day.Given halfhourly data, this requires that at least 32 valid data samples are available from the eddy-covariance data set.Once daily mean fluxes have been calculated these are used to estimate monthly mean statistics.A monthly mean is calculated if at least two-thirds of the days of a month contained valid values.This approach was chosen as the data gaps might introduce biases for daily and monthly values and it was found that the calculated error statistics could be largely influenced by a few dates with insufficient reference data.Therefore, the chosen approach provides a traceable procedure to provide reference data for different temporal resolutions.

Results
The HOLAPS validation results are summarized in the following.We hereby focus on the accuracy of the surface energy and water fluxes estimated by HOLAPS and evaluate the surface net radiation (R N ), solar radiation (R g ) and the surface latent (LE) and sensible heat (H ) fluxes for all experiments.

Evaluation of surface net radiation (R N )
The estimated surface net radiation from all 48 stations is compared against the corresponding measurements from FLUXNET in Fig. 4 for the CTRL experiment and all FLUXNET stations.Overall, HOLAPS provides very accurate estimates of R N at hourly and daily timescales.The correlation between reference data and HOLAPS is r = 0.96 (0.91) for hourly (daily) data.All correlations are significant (p < 0.05).The corresponding RMSD is 54.5 (27.2) W m −2 for hourly (daily) data with almost no bias.While the correlation coefficients for the different CTRL simulations are very high (r > 0.95), the correlation coefficients for the experiments using METEOSAT or GRIDSAT radiation are lower, still amounting to r > 0.8 for most cases.Only minor differences can be observed between the RMSD and cRMSD, which indicates that the hourly estimates of R N have only a small bias.
The accuracy of the daily and monthly net surface radiation shows a picture similar to the hourly values (see Figs. D1  and D2).The RMSD for the daily fluxes ranges between 18 and 52 W m −2 for the majority of the results and correlations are typically larger than r = 0.95.In the cases where satellite data are used as a radiation driver, the RMSD also increases and the correlation coefficient reduces.However, for monthly mean fluxes (Fig. D2) the discrepancy between CTRL simulations and the METEOSAT and GRIDSAT experiments reduces.

Evaluation of surface solar radiation flux (R g )
As shown before, major uncertainties in the surface net radiation flux are introduced by using satellite radiation products within HOLAPS.The accuracy of the radiation data it-self are therefore investigated at the FLUXNET stations in the following.Figure 6 shows the RMSD and cRMSD for hourly surface global radiation fluxes.For the CTRL simulations, the deviations are close to zero as these experiments are based on the same radiation data as that used as reference.Minor deviations still occur in these cases as the FLUXNET measurements are not available at exactly the same time steps as HOLAPS simulations.As HOLAPS interpolates the driver data to equal time steps, small interpolation differences might occur, which result in non-zero RMSD values.
The RMSD of the satellite radiation data (METEOSAT, GRIDSAT) ranges between 75 and 143 W m −2 at hourly timescales.This is partly related to a negative bias between the FLUXNET radiation data and the satellite radiation data.Thus, the deviations in the radiation data have by far the strongest effect on the surface net radiation flux and are also likely to affect the surface turbulent heat flux estimates, which will be analyzed subsequently.

Evaluation of latent (LE) and sensible (H ) heat fluxes
The overall relationship between HOLAPS latent heat flux estimates and FLUXNET measurements is illustrated in In principle, the accuracy of the results obtained might depend on additional factors, e.g., the land cover type, the cloudiness of the sky or the local time.Additional analysis of the HOLAPS results were therefore performed to analyze in more detail the impact of these additional factors.
In order to explore if the model performance is influenced by the biome types, the overall HOLAPS error statistics across biomes are shown in Figs.E1 and E2.It can be seen that the performance of HOLAPS is generally stable across biomes.Relatively high RMSD (∼ 60 W m −2 ) was found over croplands, deciduous broadleaf forests and savannas.Michel et al. (2016) investigated the accuracy of surface latent heat flux at specific times of a day.We therefore also investigated if the HOLAPS error statistics vary between daytime and nighttime compared to the entire day.The day and night separation was based on a global radiation threshold of 20 W m −2 as suggested by Reichstein et al. (2005).Figures E3 and E4 show the HOLAPS latent heat flux error statistics over daytime and nighttime.Compared to full-day statistics (r = 0.87, RMSD = 51.2W m −2 ),

Summary of HOLAPS accuracies
So far we have summarized the overall accuracies of HO-LAPS for the different experiments.As the HOLAPS framework is designed to be used at the global scale with a maximum of satellite and reanalysis data as drivers, we summarize in the following the accuracy of the HOLAPS results for the GRIDSAT_G experiment, which corresponds to the case where only satellite and reanalysis drivers are used for HOLAPS flux estimates.Results are compared against the accuracy of the CTRL_G experiment that exclusively uses FLUXNET station data and the same stations.The overall accuracies at hourly, daily and monthly timescales for these two experiments are summarized in Table 3.On monthly timescales, the results for the latent heat flux of the CTRL simulations and GRIDSAT-based estimates are rather comparable.The correlation is r = 0.80 and r = 0.81 and RMSDs are 25.5 and 26.3 W m −2 for the GRIDSAT_G and CTRL_G experiments, respectively.However, at the hourly and daily timescales the RMSD can be 10-20 % larger for the GRIDSAT_G experiment than for the CTRL_G experiment, which is likely to be a result of the uncertainties of the surface shortwave radiation fluxes.
The accuracy of the two surface solar radiation data sets was estimated for the stations that were located within the Meteosat footprint.The RMSD and correlations for R g are summarized in Table 3 as well.For the METEOSAT experiment, the hourly (daily, monthly) RMSD for the surface solar radiation flux is 83.9 (24.7, 15.3) W m −2 , whereas it is 109.6 (52.9, 31.8)W m −2 for GRIDSAT.

Discussion
The HOLAPS framework provides estimates of surface net radiation and latent heat flux at accuracies that are comparable to those obtained in other studies (Ershadi et al., 2014;McCabe et al., 2016;Miralles et al., 2016).It was found that the major source of uncertainty is the surface solar radiation data used as a forcing.When using tower only measurements (CTRL), the RMSD of HOLAPS latent heat flux is 51.2 (30.7)W m −2 for hourly (daily) fluxes.Michel et al. (2016) and Miralles et al. (2016) evaluated the performance of four different algorithms to estimate the surface latent heat flux, within the WACMOS-ET project, using either tower-based forcings or satellite data.As this is probably one of the most comprehensive studies existing, we compare our results against results from that study.The RMSD for the algorithms investigated in the study of Michel et al. (2016) ranges between 40.8 and 88.5 W m −2 when comparing their results at 3-hourly time step and using tower data as a driver.At daily timescales, the RMSD obtained for the same four algorithms ranged between 22.7 W m −2 and 52.2 W m −2 .Correlations were found to range between 0.76 and 0.88 (0.66 and 0.78) for 3-hourly (daily) values.Under the support of the GEWEX LandFlux project, McCabe et al. (2016) evaluated the same methods but with a different number of tower stations.They found that the correlations range from 0.71 to 0.85, and RMSD range from 61 to 101 W m −2 for towerbased 3-hourly data.Similar statistic scores (RMSD between 64 and 105 W m −2 ) have also been reported by Ershadi et al. (2014), who also evaluated similar methods (SEBS, PT-JPL, PM, advection-aridity) with tower-based half-hourly or hourly data.For HOLAPS we have provided the accuracy measures when using all data samples (all stations + all years) at once.These were provided in Table 3.The HO-LAPS hourly (daily) RMSD is 51.2 (30.7)W m −2 with correlations of r = 0.87 (r = 0.79).However, these values are not exactly comparable with the study of Miralles et al. (2016) as (a) the HOLAPS statistic is based on hourly values instead of 3-hourly values for the WACMOS-ET project.Further, the information provided by Michel et al. (2016) is given as the mean value from results of all investigated stations.Thus, instead of calculating the RMSD for all data samples, these authors calculated first the error statistics and then provided the mean skill score.When following an approach similar to the 48 stations investigated in the present study, the mean RMSD of HOLAPS corresponds to 46.6 (26.5)W m −2 with mean correlations of r = 0.89 (0.85) for hourly (daily) timescales.Thus, following an approach similar to the one by Michel et al. (2016) the results of the present study are very similar to those of WACMOS-ET.
Similar differences are also obtained when using satellite data as a driver for the latent heat flux estimates.The RMSD obtained for 3-hourly (daily) estimates by Michel et al. (2016) ranges from 47.6 to 88.5 (24.5 to 59.0) W m −2 while HOLAPS hourly (daily) RMSD is 62.3 (29.1)W m −2 with correlations of r = 0.79 (r = 0.72), whereas Michel et al. (2016) found correlations of 0.69 < r < 0.82 (0.59 < r < 0.79) for 3-hourly (daily) comparisons.Overall, HOLAPS seems to provide improved correlations, which might be due to the enhanced temporal resolution of HOLAPS.It needs to be emphasized, however, that results of the present study are not fully comparable with Michel et al. (2016), due to the different temporal sampling, and the different number of stations investigated (N = 48 in this study instead of N = 24).
Overall, a small bias was observed, for both the simulations with flux-tower and satellite forcings (see Table 3).While the CTRL and GRIDSAT experiments differ on hourly and daily timescales, the RMSD for the monthly results is very similar.This indicates that the uncertainties due to the large-scale forcing are minimized at longer timescales.
Replacing station precipitation data with the TMPA largescale satellite forcing as well as using ERA-Interim for temperature and wind speed has a minor effect on the accuracy of the results obtained.By far the largest uncertainties are introduced when using satellite-based surface solar radiation data, whereas similar accuracies are obtained using either the METEOSAT or GRIDSAT data.The accuracy for the surface solar radiation flux from METEOSAT was found to have an RMSD of 83.9 (24.7)W m −2 for hourly (daily) timescales using the FLUXNET stations located within the Meteosat footprint (N = 19), which is slightly larger than the daily RMSD of 17.9 W m −2 reported by Müller et al. (2015) based on BSRN observations.As a further improvement of the surface solar radiation flux is expected to improve the latent heat flux estimates, a thorough investigation of the impact of different surface solar radiation data set will be performed in a future study.This could then also include the analysis of reanalysis-based radiation data, which was excluded from the present study, as Posselt et al. (2012) had already shown that the METEOSAT radiation data used in the present study has an overall better agreement with ground measurements than the ERA-Interim reanalysis radiation data.Overall, best results were obtained for clear-sky conditions.Decreasing performance of HOLAPS estimates was observed for increased cloudiness, which is likely to be caused by the increased uncertainties in the satellite-based radiation data under cloudysky conditions.No systematic differences between different biome types could be identified in this study.A more comprehensive sensitivity analysis of HOLAPS to different biomespecific model parameters might be subject of a further study, where the vegetation parameter of each biome will be perturbed and the relevant HOLAPS performance will be assessed.

Conclusions
This study has introduced a new framework for the estimation of high-resolution land surface water and energy fluxes, HOLAPS v1.0.The framework was developed to make use of existing satellite data records and to allow for the generation of temporal and spatial high-resolved and consistent quasi-global water and energy fluxes.Key features of the HOLAPS framework comprise internally consistent estimation of shortwave and longwave radiation fluxes; capability to directly use top-of-atmosphere radiances for surface solar flux estimations; constrained surface fluxes using a mixed boundary layer model in combination with the surface flux estimates; flexible framework for the generation of high-resolution land surface energy and water fluxes that allows one to use a multitude of different land surface schemes within the same framework.
This study analyzed the accuracy of HOLAPS v1.0 using data from 48 eddy-covariance towers.A sensitivity analysis was performed to investigate the tradeoff in using satellite data as drivers instead of locally measured tower-based data.
The results of this study can be summarized as follows: -The accuracy of the HOLAPS surface fluxes was found to be comparable or even better than results obtained in other studies for the surface net radiation as well as turbulent fluxes.
-The hourly (daily) RMSD for the surface net radiation flux was 54.5 (27.2) W m −2 with correlations of r = 0.96 (r = 0.91) when using tower data as drivers for HOLAPS.
-Using satellite and reanalysis data as only drivers, the RMSD and correlations were found to be 61.8W m −2 and r = 0.79 (33.1, r = 0.71) for the latent heat flux -Accuracy of turbulent flux estimates decreases with increasing cloudiness due to higher uncertainties in the surface solar radiation flux, which is consistent with previous studies.
-The largest uncertainties resulted from the uncertainties of the surface solar radiation flux.However, on monthly timescales, these uncertainties were minimized, which indicates that comparable accuracies can be obtained when using satellite-based drivers instead of local in situ data.
A first quasi-global data set generated using HOLAPS v1.0 is planned to be released to the scientific community after a thorough validation and cross-comparison against other data sets, e.g., the LandFlux-Eval (Mueller et al., 2013) data.Further improvements of the HOLAPS framework will comprise the capability to assimilate land surface temperature data from geostationary satellite observations to better constrain the surface latent heat flux estimates as well as the usage of new satellite observations, e.g., provided by the new SEN-TINEL series of satellites.Recent advances in available computational resources allow for the first time to exploit these high spatial resolution sensors at a global scale and might lead to operational services provided, e.g., in the frame of Copernicus services.
A major constraint is nevertheless the lack of consistent and harmonized geostationary satellite data records.The mosaic of geostationary satellites, known as GEORING, is currently operated by individual space agencies and so far no long-term climate or operational data set of harmonized and well-intercalibrated geostationary radiance and brightness temperature data is available at the original sensor resolution.The GRIDSAT data set, used in the present study is currently the only long-term GEORING data set available, but is limited in its spatial resolution.Further developments towards Fundamental Climate Data records from geostationary satellite data are therefore required.
Further studies using HOLAPS will investigate the potential to use the novel SENTINEL data streams and to further reduce the dependency on reanalysis data by using, e.g., the total column water vapor information from satellite data and exploit the potential of internally consistent land surface parameters currently developed, e.g., by different European projects (QA4ECV, MULTIPLY).
While the present study provides a sensitivity analysis of using the HOLAPS framework with different forcing data, it would be important to conduct further in-depth studies to disentangle the different components of the overall error budget (model uncertainties, forcing uncertainties, scale mismatches, reference data uncertainty), which still remains a major challenge to be addressed by the research community.

Code availability
The HOLAPS code used for this manuscript can be accessed via https://github.com/pygeo/holaps.PBL module ) is estimated using the MAGIC radiative transfer model (Mueller et al., 2009).The shortwave surface downwelling solar flux (R g ) for all sky conditions is then obtained from the clear-sky downwelling solar flux and the clear-sky index k as (Posselt et al., 2012) (B3) The clear-sky index is related to CAL through the following relationship (Hammer et al., 2003) where a = 2.0667, b = −3.667,c = 1.6667.

Longwave surface radiation fluxes
The longwave surface downwelling radiation flux (L ↓ ) depends on the near-surface moisture and temperature profile as well as the cloud coverage.The clear-sky longwave downwelling radiation flux L ↓ slab is calculated using the PBL model (Margulis and Entekhabi, 2001).L ↓ slab is then corrected for cloud coverage as (Brubaker and Entekhabi, 1995) (B5)

B3 Soil module
The surface temperature T S [K] is obtained by a revised force restore approach (Ren and Xue, 2004) as where A (K) is the diurnal temperature amplitude of T S , ) is the thermal inertia coefficient and is the thermal inertia, which is estimated as function of soil moisture conditions (Murray and Verhoef, 2007) and ω = 2π 86 400 is the diurnal angular frequency.The parameters B and a in Eq. ( B6) are set to a = 0.45π and B = 0.158 (Ren and Xue, 2004).The prognostic equation for the deep soil layer temperature T d is where d is the soil temperature damping scale depth with typical values on the order of d = 0.15 (m).The lapse rate between the mean surface and deep-layer temperature γ S (K m −1 ) is estimated from the differences between T S and T d and τ = 86 400 (s) is the time period, 1 day in our case.

B4 Water balance module
The surface water balance is defined as The soil moisture dynamics is calculated using a multilayer soil scheme, discretized into five layers.The soil layers have a thickness of dz = [0.05,0.1, 0.25, 0.6, 1.0] (m).Soil moisture fluxes between the different soil layers are simulated by solving numerically the Richards equation (Richards, 1931) whereas only vertical moisture fluxes are considered: The water fluxes between the different soil layers is solved using a numerical approach.The net soil water flux in a soil layers is hereby determined by the fluxes into and from the layers above any below, whereas the model allows for both downward (percolation) and upward (capillary rise) fluxes.Surface runoff Q is obtained as the excess of water that can not infiltrate the soil when maximum infiltration capacity is reached.The relationship between volumetric soil moisture content and soil suction head ψ is calculated using the model of van Genuchten (1980).
The water interception by the canopy is estimated by (Valente et al., 1997) where ET I = λ −1 LE I is the transpiration from the canopy interception storage and D is the through fall and drainage of water from the canopy layer to the soil.

B5 Turbulent flux module
For a vegetated patch with fractional vegetation coverage f c the surface latent heat flux is calculated as the weighted sum of the evaporation from soil (LE S ) and the transpiration from the canopy (LE C ) as well as evaporation from water intercepted by the canopy layer (LE I ) as where w I = (I /I max ) b is a weighting factor dependent on the current canopy interception storage I , the potential maximum interception storage I max ( ) (von Hoyningen-Huene, 1981) and an empirical parameter b = 0.5 (Chen and Dudhia, 2001).The vegetation cover fraction f c is obtained from leaf area index ( ) as (Norman et al., 1995) which assumes a random leaf distribution with spherical leaf angle distribution.The different latent heat flux components in Eq. ( B11) are then estimated using the Priestley-Taylor approach as where α pt = 1.26 is the Priestley-Taylor parameter for equilibrium evapotranspiration and , γ are the slope of the water vapor saturation curve and psychrometer constant (Pa K −1 ), respectively.The inhibition function 0 ≤ [φ] ≤ 1 describes the reduction of LE due to limiting factors like radiation, temperature and soil moisture.The soil net radiation is estimated as (Norman et al., 1995) R N,S = R N e 0.9 ln(1−f c ) (B14) and the canopy net radiation is then calculated as The sensible heat flux is estimated as where the aerodynamic surface resistance r a (s m −1 ) is calculated as where k ≈ 0.41 is the von Karman constant and u z corresponds to the win speed at canopy height and is obtained from wind speed data assuming a logarithmic wind profile and a displacement height d corresponding to two-thirds of the vegetation height (Maidment, 1993).The stability correction functions m,h are calculated after (Paulson, 1970) using the Richardson number Ri as an indicator for atmospheric stability.The roughness lengths for momentum and heat (z 0,m , z 0,h ) are parameterized for each land cover type (Table B1).

Surface inhibition functions
The canopy inhibition function 0 ≤ ϕ c ≤ 1 is defined as (Chen and Dudhia, 2001) where R r is a function of surface air temperature and pressure, C h is the surface exchange coefficient for heat and moisture and R C is the canopy resistance, given as r rad , where r rad is a radiation-specific parameter W m −2 and r min and r max are the minimum and maximum canopy resistance (s m −1 ), which are all land cover specific parameters (Table B1).The relative degree of soil saturation is given by and w 0 =1, w f = 800 and µ = 12 are empirical parameters (Anderson et al., 2007).f T a and f ↓ S are based on Chen and Dudhia (2001).

B6 Planetary boundary layer module
The prognostic equations of the PBL model are given by (Kim and Entekhabi, 1998;Smeda, 1979) with the proportionality constant a = 10 −5 (s −1 ) (Smeda, 1979).Alternative approaches to simulate the radiative cooling have been proposed (Kim and Entekhabi, 1998;Margulis and Entekhabi, 2001).The relationship between PBL air temperature (T ) and θ is given by with R/c p ≈ 0.286 for air.The details of the model formulations are based on Smeda (1979) and are given as follows: with δ = 0 in stable conditions and δ = 1 in unstable conditions.We set ζ = 0.01 to ensure a realistic collapse of the PBL (Kim and Entekhabi, 1998).
During daytime, the growth of the PBL is determined by the right side in Eq. (B30).During the transition between unstable and stable conditions, the PBL collapses because of turbulence dissipation.The PBL height during this transition phase is given as (Smeda, 1979) when assuming that H top = 0. Equation (B29) is applied in this transition phase until The mixed layer is capped by an inversion with inversion strength δ θ m (K), which determines the entrainment of overlying dry air from the free atmosphere as (McNaughton and Spriggs, 1986) Dry air entrainment causes the inversion strength itself to change according to where γ θm K m −1 is the potential temperature lapse rate above the PBL and is assumed to be constant.

B7 Model parameterization
The land cover specific model parameters are summarized in Table B1.They are based on the publications of Chen andDudhia (2001) andHagemann (2002).

B8 Interpolation methods
Different interpolation approaches are used to interpolate the input data onto the HOLAPS computational grid and time step.The used techniques are summarized in Table B2 for each of the HOLAPS drivers.The nearest neighbor remapping as well as bilinear interpolation are currently used for spatial remapping.The temporal interpolation is based on a linear interpolation of measurements (y 1 , y 2 ) between two observation times (t 1 t 2 ) whereas the weight w depends on the sampling times and the actual model time step.
To handle data gaps, the HOLAPS framework currently provides the following options: ignore: the data gap is ignored and filled by interpolation.
-last_valid: last valid value of the variable is used and the data gap is filled with this value.
-last_valid_same_time: use the last valid data at the same time of the day.This option is in particular useful for data that show a strong diurnal dynamics (e.g., radiation).In that case, using the last valid value would lead to erroneous diurnal forcing data when data gaps of a few hours occur, which can be quite often the case when using, e.g., FLUXNET data.
climatology: a climatological mean annual cycle i used for the calculations.
Interpolation methods can be easily changed by the user in configuration files.

The special case of radiation data
No direct interpolation is performed for the radiation data, as a linear approximation might not be sufficient to capture the diurnal cycle of the surface solar radiation flux.Instead, the clear-sky index (k) is interpolated in time and then used to calculate R g using Eq.(B3).

Figure 1 .
Figure 1.HOLAPS drivers, estimated surface fluxes and surface fluxes and modules.

Figure 2 .
Figure 2. Forcing data and variable interdependencies in the HOLAPS model.Only major output variables are illustrated.Details for model formulations can be found in Appendix B.

Figure 3 .
Figure 3. Distribution of FLUXNET stations used in this study.Light green corresponds to latitudes between 50 • N and 50 • S, which corresponds to the coverage of the TMPA precipitation data (see text).Stations in red cannot be used when forced with TMPA data.Light orange indicates approximate coverage of Meteosat data.

Figure 4 .Figure 5 .Figure 6 .
Figure 4. Comparison of surface net radiation flux (R N ) between FLUXNET measurements and HOLAPS estimates for the CTRL experiment: (a) hourly and (b) daily timescales.Colors indicate the frequency of occurrence of values (data density).Units in W m −2 .

Figure 7 .
Figure 7.Comparison of HOLAPS latent heat flux for (a) hourly and (b) daily timescale for the CTRL experiment using results from all stations and years.Units in W m −2 .

Fig. 7 .
The RMSD is 51.2, 30.7 and 26.3 W m −2 for the hourly, daily and monthly flux estimates, respectively, for the CTRL_G simulations.The correlation coefficient is 0.87 for hourly data, 0.79 for daily and 0.81 for monthly data.Error statistics for all experiments are provided in Fig.8.The increased uncertainty in the surface solar radiation and thus R N has a direct effect on the accuracy of the latent heat flux estimates.Correlation coefficients are the smallest for the experiments that use satellite surface solar radiation data.However, the correlations are still high with r > 0.74 for most of the stations and experiments.The RMSD for the CTRL simulations ranges between 35 and 52 W m −2 for the majority of the cases.The largest RMSD is observed for the METEOSAT and GRIDSAT experiments.However, results from the experiments when replacing the air temperature and wind speed with reanalysis data show that this also introduces uncertainties in the latent heat flux estimates.The RMSD ranges between 40 and 62 W m −2 for these experiments.Corresponding results for daily and monthly timescales are provided in Figs.D3 and D4.The overall error statistics for the sensible heat flux in the CTRL_G simulations are shown in Fig. 9.The RMSD ranges from 79.1 W m −2 (hourly) to 36.0 W m −2 (daily).The error statistics for all experiments are shown in Fig. 10 and show a result similar to the latent heat flux error statistics with worse statistics for the experiments with satellite radiation data as a forcing.The daily and monthly comparison results are shown in Figs.D5 and D6.

Figure 9 .Figure 10 .
Figure 9.Comparison of HOLAPS sensible heat flux for (a) hourly and (b) daily timescale for the CTRL experiment using results from all stations and years.Units in W m −2 .

A
. Loew et al.: High-resolution land surface fluxes from satellite and reanalysis data Appendix A: Acronyms Acronyms used throughout the text are summarized in the following table.

g
Shortwave

Figure E3 .
Figure E3.Error statistics for HOLAPS latent heat flux over daytime: (a) comparison using results from all stations and years, (b-d) box plots of validation statistics that are calculated at each station.

Figure E4 .
Figure E4.Error statistics for HOLAPS latent heat flux over nighttime: (a) comparison using results from all stations and years, (b-d) box plots of validation statistics that are calculated at each station.

Figure E5 .
Figure E5.Error statistics for HOLAPS latent heat flux over clear-sky condition: (a) comparison using results from all stations and years, (b-d) box plots of validation statistics that are calculated at each station.

Figure E7 .
Figure E7.Error statistics for HOLAPS latent heat flux over cloudy-sky condition: (a) comparison using results from all stations and years, (b-d) box plots of validation statistics that are calculated at each station.

Table 1 .
Overview of data sets used as drivers for HOLAPS.

Table 2 .
List of performed model experiments.Includes the number of stations and station years as well as the data source: F is FLUXNET data; S is satellite data for precipitation and radiation; additional data from satellites for albedo and LAI, and from ECMWF reanalyses for temperature, total column water vapor, and wind speed.
Control simulations are conducted for all of these different spatial domains.As a consequence a total of four different control simulations with a different number of stations are conducted.All the other experiments were also performed for these different spatial subsets where applicable.The differences between the same experiment type, at different spatial domains provides additional information on the variability of the error metrics as a function of the number of FLUXNET stations used.Table2summarizes all experiments conducted and the number of stations and simulation years.

Table 3 .
Overall HOLAPS accuracies for R N , LE and R g , at hourly (h), daily (d) and monthly (m) timescales for the CTRL, GRIDSAT and METEOSAT experiments.

Table A1 .
Acronyms used throughout the manuscript.

Table B1 .
Land cover specific parameters.

Table B2 .
Summary of spatial and temporal interpolation methods used within the HOLAPS framework for different driver variables.

Table C1 .
Loew et al.:High-resolution land surface fluxes from satellite and reanalysis data Appendix C: FLUXNET stations List of FLUXNET stations investigated.The coverage term specifies the location of each FLUXNET station.±50•refers to the station being within the latitudes 50 • N, 50 • S, while Meteosat indicates the station is within the coverage of Meteosat.For details on the spatial coverage see Fig.3.