Modelling sub-grid wetland in the ORCHIDEE global land surface model: evaluation against river discharges and remotely sensed data

. The quality of the global hydrological simulations performed by land surface models (LSMs) strongly depends on processes that occur at unresolved spatial scales. Approaches such as TOPMODEL have been developed, which allow soil moisture redistribution within each grid-cell, based upon sub-grid scale topography. Moreover, the coupling between TOPMODEL and a LSM appears as a potential way to simulate wetland extent dynamic and its sensitivity to climate, a recently identiﬁed research problem for biogeochem-ical modelling, including methane emissions. Global evaluation of the coupling between TOPMODEL and an LSM is difﬁcult, and prior attempts have been indirect, based on the evaluation of the simulated river ﬂow. This study presents a new way to evaluate this coupling, within the ORCHIDEE LSM, using remote sensing data of inundated areas. Because of differences in nature between the satellite derived information – inundation extent – and the variable diagnosed by TOPMODEL/ORCHIDEE – area at maximum soil water content, the evaluation focuses on the spatial distribution of these two quantities as well as on their temporal variation. Despite some difﬁculties in exactly matching observed lo-calized inundated events, we obtain a rather good agreement in the distribution of these two quantities at a global scale. Floodplains in the model, is The difﬁculty of reproducing the year-to-year variability of the observed inundated area (for instance, the trend by the end of 90s) is underlined. Classical indirect based on comparison between simulated and river and after

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.
Abstract. The quality of the global hydrological simulations performed by land surface models (LSMs) strongly depends on processes that occur at unresolved spatial scales. Approaches such as TOPMODEL have been developed, which allow soil moisture redistribution within each grid-cell, based upon sub-grid scale topography. Moreover, the coupling between TOPMODEL and a LSM appears as a potential way to simulate wetland extent dynamic and its sensitivity to climate, a recently identified research problem for biogeochemical modelling, including methane emissions. Global evaluation of the coupling between TOPMODEL and an LSM is difficult, and prior attempts have been indirect, based on the evaluation of the simulated river flow. This study presents a new way to evaluate this coupling, within the ORCHIDEE LSM, using remote sensing data of inundated areas. Because of differences in nature between the satellite derived information -inundation extent -and the variable diagnosed by TOPMODEL/ORCHIDEE -area at maximum soil water content, the evaluation focuses on the spatial distribution of these two quantities as well as on their temporal variation.
Despite some difficulties in exactly matching observed localized inundated events, we obtain a rather good agreement in the distribution of these two quantities at a global scale. Floodplains are not accounted for in the model, and this is a major limitation. The difficulty of reproducing the yearto-year variability of the observed inundated area (for instance, the decreasing trend by the end of 90s) is also underlined. Classical indirect evaluation based on comparison between simulated and observed river flow is also performed and underlines difficulties to simulate river flow after coupling with TOPMODEL. The relationship between inundation and river flow at the basin scale in the model is analyzed, using both methods (evaluation against remote sensing data and river flow). Finally, we discuss the potential of the TOP-MODEL/LSM coupling to simulate wetland areas. A major limitation of the coupling for this purpose is linked to its ability to simulate a global wetland coverage consistent with the commonly used datasets. However, it seems to be a good opportunity to account for the wetland areas sensitivity to the climate and thus to simulate its temporal variability.

Introduction
Land surface processes control the partition of incoming radiative energy into sensible and latent heat fluxes, and heat storage. This partition depends on the available net energy and on soil moisture, which limits transpiration in dry regions (e.g. Teuling et al., 2009). Many global and regional modelling studies have demonstrated the influence of soil moisture on climate variability and predictability (e.g. Douville, 2003;Seneviratne et al., 2006;Vautard et al., 2007). The spatiotemporal variability of land surface processes is usually represented as a boundary condition to the atmosphere by land surface models (LSMs). In this context, soil moisture is a key variable of LSMs, and has motivated intercomparison projects (GSWP-2; e.g. Guo et al., 2007), but the evaluation of soil moisture in models remains challenging. One issue is that in many LSMs, soil moisture is not calculated as a physical state variable, but rather as an index for predicting evapotranspiration and runoff, thus making it difficult to evaluate with observations. Then, there is a scale discrepancy between a model grid, and either in-situ or satellite measurements of soil moisture. Finally, remote sensing data do not provide a direct observation of soil moisture and have errors often poorly characterized (Schumann et al., 2009). Alongside with above issues regarding the evaluation of simulated soil moisture at the resolution of global climate models (typically from 50 to 300 km), an improvement in the treatment of subgrid-scale variations is also needed. Indeed, several studies showed a better simulation of interception, regional water budgets and streamflow when accounting for a subgrid description of precipitation, soil moisture, drainage and/or runoff (Decharme and Douville, 2007;Wood et al., 1998). One strategy for treating subgrid runoff in LSMs relies on the introduction of the simple hydrological model, TOPMODEL, which was initially developed to account for saturation excess runoff (Dunne, 1978) owing to the variable contributing area concept (Beven and Kirkby, 1979). The use of TOPMODEL in LSMs (Decharme et al., 2006;Decharme and Douville, 2005;Habets and Saulnier, 2001;Famiglietti and Wood, 1994;Koster et al., 2000;Gedney and Cox, 2003) relies on soil moisture redistribution within each grid-cell based upon a topographic index, what allows to determine the saturated fraction of each grid cell, given grid-cell average soil moisture computed by the LSM.
Wetlands cover a large diversity of ecosystems, including peatlands, marshes, and swamps (Reichardt, 1995), but have in common a strong impact on climate, both directly by humidifying the atmosphere (Krinner, 2003) or indirectly as the single largest source of atmospheric methane (CH 4 ) (Forster et al., 2007). Until recently, studies related to current wetland CH 4 emissions (e.g. Walter et al., 2001;Zhuang et al., 2004, Wania et al., 2010 considered fixed wetland extent and introduced a vertical hydrologic variability through the change of the water table depth. But recent findings have shown that the spatial and temporal variability of wetland areas may play an important role in controlling seasonal and interannual changes of CH 4 emissions from wetlands (Bloom et al., 2010;Ringeval et al., 2010) and thus changes in global CH 4 growth rate (Bousquet et al., 2006). This recent finding overturns much previous thinking that the source areas of CH 4 were unchanging, or slowly varying on millennial time scales (Kaplan et al., 2006). Thus, accounting for wetland dynamics in land surface modelling has emerged as a major research problem. This is particularly true in the context of anthropogenic climate change, as changes in wetland extents are expected to impact the future evolution of CH 4 emissions with feedbacks on global climate (Gedney et al., 2004;Ringeval et al., 2011). Because of the sensitivity of CH 4 emissions to wetland extent, it is particularly important to have a well-constrained estimate of the dynamic response of wetland extent to climate.
In this context, it is necessary to account for wetland dynamic through a process-based approach rather than empirically (Krinner, 2003;Riley et al., 2011). Given that TOP-MODEL enables to simulate the dynamic of saturated areas, its coupling with an LSM is more and more used to diagnose the wetland extent dynamic (Gedney and Cox, 2003;Ringeval et al., 2011) or the presence of peatlands . However, the interest of this coupling to diagnose wetland extent at global scale remains poorly evaluated. The evaluation of the coupling between TOPMODEL and an LSM is usually indirect, based on comparisons between simulated and observed river discharge (e.g. Decharme and Douville, 2007). Recently, a global, multi-year dataset quantifying the monthly distribution of flooded areas at ∼ 25 km resolution has been generated from multiple satellite observations optimized specifically for surface water detection Prigent et al., 2001Prigent et al., , 2007. This dataset represents an unprecedented source of information to evaluate directly the coupling between TOPMODEL and a LSM. Given that the Prigent et al. (2007) data have already been used successfully to approach wetland dynamics Melton et al., 2012, WETland and Wetland CH4 Intercomparison of Models Project, http://arve.epfl.ch/research/ wetchimp/), the evaluation of the TOPMODEL/LSM coupling against these data will allow to better characterize the ability of such physically-based approach to describe the subgrid area of wetland. Through the study of the year-to-year variability, this evaluation will also enable us a better characterizing the TOPMODEL/LSM ability to describe the wetland extent sensitivity to climate.
In this study, a subgrid soil moisture redistribution scheme based on TOPMODEL is implemented into the ORCHIDEE (ORganizing Carbon and Hydrology In Dynamic Ecosys-tEms) global LSM (Krinner et al., 2005), which is run globally at 1 degree resolution. The ORCHIDEE model is first refined to account for the soil water freezing and thawing cycles, a necessary step for modelling wetland dynamics in cold regions. These developments are described in Sect. 2. The experimental set-up is presented in Sect. 3. Then, the B. Ringeval et al.: Modelling sub-grid wetland in the ORCHIDEE global land surface model 943 evaluation of the current and new versions of ORCHIDEE is performed through the classic simulation of streamflow and comparisons with observations of river discharge at gauging stations (Sect. 4). The coupling of TOPMODEL and OR-CHIDEE is then further evaluated through the comparison of predicted saturated areas against the satellite products of Papa et al. (2010) (Sect. 5). The discussion has two objectives (Sect. 6). First, it focuses on the ability of the new version of ORCHIDEE to capture the relationship between streamflow and inundation extent at the scale of river basins. Second, the suitability and limitations of an approach like TOPMODEL to simulate wetland areas at large scales are discussed. Finally, Sect. 7 concludes the study.

Current soil hydrology model of ORCHIDEE
The ORCHIDEE model is a process-based dynamic global vegetation model developed for carbon cycle applications and as the land component of the IPSL coupled model (Krinner et al., 2005) and using the approach of Decharme et al. (2006). The model includes parameterizations of canopy physiology (photosynthesis and canopy conductance) that are intimately linked to energy and water fluxes and operated at a time scale of 30 min. Soil hydrology is modelled using a semi-empirical approach (Ducoudré et al., 1993).
The soil moisture reservoir has a fixed depth of 2 m, representing essentially the root zone, and a spatially uniform volumetric soil water holding capacity (SWC = 150 kg m −3 , fixed value for all grid-cells) that corresponds to the maximal amount of water that plants can extract. In ORCHIDEE, the evolution of soil moisture is computed daily with a twolayer soil model (Ducoudré et al., 1993), composed of a 2 m thick soil reservoir and a surface layer of variable thickness receiving incoming precipitation. The water content of each layer is updated by inputs from snowmelt and rainfall not intercepted by the canopy, and by losses from soil evaporation, transpiration, sublimation, deep drainage, and surface runoff. Drainage from the top layer depends on its volumetric soil moisture (θ ) in a nonlinear way. It remains small until θ reaches 0.75, and increases strongly above this threshold (De Rosnay and Polcher, 1998). Drainage from the lower layer (base flow) is parameterized to equal to 95 % of the overflow runoff, simply defined as the excess water when the soil is over its maximum capacity, as in a classic bucket model (Manabe, 1969).
ORCHIDEE represents heterogeneous vegetation using a "mosaic" in each grid cell, given a set of 12 Plant Functional Types (PFT) and bare soil. The total water flux from the land surface to the atmosphere in each grid cell is computed as the sum of snow sublimation, soil evaporation, transpiration by plants and evaporation of water intercepted by the canopy over all the PFTs. The soil hydrology module provides water limitations to the carbon module, and thus governs carbon allocation, litter and soil carbon decomposition. A separate water stress function acting on photosynthesis and plant transpiration is calculated from the convolution between the water content of the two soil layers and the root profile. This stress function varies between 0 and 1, respectively at the minimum (ω min ) and maximum soil water content (ω max ) (McMurtrie et al., 1990). (ω max -ω min ) defines the potential maximum soil water content (SWC). We made the following modifications to the current OR-CHIDEE hydrological processes (Krinner et al., 2005): -Transpiration does not stop when water is intercepted by the canopy, i.e. we assume that stomata can still emit water vapour (N. Viovy, personal communication, 2010).
-The water-holding capacity of the canopy is increased (0.1 to 0.5 mm of water per unit of LAI) to be in closer agreement with values reported in Crockford and Richardson (2000).
-The potential maximum soil water content (SWC) is assumed to vary spatially, as a function of the surface soil texture from the Food and Agricultural Organization dataset. To do so, we computed a soil saturation capacity (ω sat ) and a wilting point (ω wilt ) using regression fits to soil parameter values (Noilhan and Lacarrere, 1995) as well as to soil organic matter content (Lawrence and Slater, 2007). The ORCHIDEE SWC is then defined as 60 % ω sat − ω wilt so as to keep a global mean soil water holding capacity value of 150 kg m −3 , as in the current model version, and thus bring minimum perturbations (Ducharne and Laval, 2000). SWC used as input for ORCHIDEE is given in Fig. A1 (right panel).
This modified ORCHIDEE version is the starting point upon which freeze/thaw parameterization and TOPMODEL subgrid soil moisture redistribution are included, as summarized in Fig. 1. In order to account for saturated area extension during the spring thaw and their shrinking or disappearance during the autumn freeze, we add to ORCHIDEE a simple parameterization of frozen soil water. The modified version is called ORCHIDEE-FRZ (Fig. 1). This includes taking into account the soil heat conductivity and heat capacity of liquid and frozen phases, as in Poutou et al. (2004). An apparent heat capacity representing latent heat release and uptake during freezing and thawing respectively is introduced. The thermal and hydrological schemes of ORCHIDEE are respectively discretized vertically using 7 and 2-layers. The frozen fraction of each thermal layer that varies between 0 and 1 is diagnosed from the thermal scheme then multiplied by a water content profile obtained by interpolating the 2 hydrological layers onto the 7 thermal layers. This allows diagnosing the frozen water content in each of the 7 thermal layers, keeping in mind that vertical water distribution into ORCHIDEE is not explicitly represented.

Freeze and thaw processes
The second modification introduced is the implementation of a reduction of water infiltration in frozen soils (Farouki, 1981cited by Poutou et al., 2004. This process is represented by scaling the potential infiltration with the reduction of liquid water holding capacity within the upper 20 cm of soil caused by presence of ice in the soil. where ω froz is the frozen soil water content in the upper soil top 0.20 m and P g is the water input of the soil (rain + snowmelt). In Eq. (1) the soil liquid water content is reduced when a fraction of the pores is filled up by ice, which has the consequence of increasing meltwater runoff during spring (R froz see Table 1) (Takata and Kimoto, 2000). This is a typical case of Horton runoff occurring when production of meltwater exceeds the soil saturated hydraulic conductivity (see e.g. Dunne, 1978). Finally, frozen water is discarded for computing the plant water stress function and water availability for litter degradation.

Subgrid soil water content distribution from TOPMODEL scheme
We follow the approach of Decharme et al. (2006) and Habets and Saulnier (2001) for the ISBA model and Gedney and Cox (2003) for the MOSES model to describe a subgrid soil moisture distribution into ORCHIDEE using TOPMODEL concepts (Beven and Kirkby, 1979). Following Decharme and Douville (2007), we also incorporate the bias correction of Saulnier and Datin (2004). This defines the ORCHIDEE-TOP version, which we apply globally and evaluate in this study ( Fig. 1) and which allows us to compute at each time-step the fraction of each grid-cell which reaches ω max , further noted F max . TOPMODEL was initially developed at river catchment scale. It attempts to combine the distributed effects of channel network topology and dynamic contributing areas for runoff generation (Beven and Kirkby, 1979;Sivapalan at al., 1987). This formalism takes into account topographic heterogeneities explicitly based on the spatial distribution of a topographic index λ i (m), defined at the pixel scale as follows: where a i (m) is the drainage area per unit of contour of a local pixel, i, and tanβ i , the local topographic slope, approximates the local hydraulic gradient where β i is the local surface slope. If a given pixel in a catchment has a large drainage area and a low local slope, its topographic index will be large and thus, its ability to be saturated will be high. TOPMODEL gives a relationship between the mean water deficit in a catchment (D t ), the local deficit of a given pixel (d i,t ) into the considered catchment and the topographic index: where M (m) is a parameter describing the exponential decrease of the soil transmissivity with local deficit and λ m is the mean topographic index over the catchment. For a given mean deficit over the catchment (D t ), a threshold topographic index λ th can be diagnosed in such a way that all pixels with a local topographic index λ i > λ th have no local deficit (d i,t = 0). Then, a fraction of the catchment, noted F max , defined by the pixels with no water deficit can be estimated from the partial integration of the spatial distribution of the topographic index in the catchment, noted δ: The coupling between TOPMODEL and ORCHIDEE assumes that the relationship between local soil moisture, mean deficit and topography holds within each grid-cell at the LSM resolution (Gedney and Cox, 2003). It requires the estimation of the grid-cell mean deficit from variables computed by ORCHIDEE. Following Decharme et al. (2006), we consider that the grid-cell average deficit (D t ) and the soil moisture computed by ORCHIDEE (ω soil ) are proportional for each time-step, so that the grid-cell average deficit D t can be simply expressed as here, h soil is the ORCHIDEE soil depth and D t is computed as a deficit with respect to the maximum soil water content ω max . As a result, F max corresponds to the subgrid fraction at ω max , and it varies at each time-step, being inversely proportional to the grid-cell mean deficit D t deduced from the soil water content ω soil computed for each time-step by OR-CHIDEE.
Following Decharme and Douville (2007), we also incorporate the bias correction of Saulnier and Datin (2004). This correction leads to more complex relationships than given here by Eqs. (3) and (4) for the sake of simplicity. All details can be found in Decharme et al. (2006).
As in Kleinen et al. (2012), the maximum soil water content cannot be considered as a saturated state and thus, F max cannot be directly compared to saturated fractions simulated by LSM/TOPMODEL coupling as in Decharme et al. (2006).
Given the simple representation of soil hydrology in our LSM, the assumptions used in TOPMODEL to establish the relationship between mean deficit, local deficit and topography (namely hydraulic gradient equal to surface slope, steady-state conditions, and exponential decrease of hydraulic conductivity) are not explicitly accounted for into ORCHIDEE. The above coupling allows for the transition from the notion of vertical water flux, which is present in the LSM, to that of the horizontal water fluxes, on which TOP-MODEL is based.
In this new version, surface runoff is computed as the sum of Dunne runoff from precipitation falling on the fraction at maximum soil water content, F max , and of Horton runoff over frozen soil (Table 1). For partially frozen soils, we thus introduce an effective fraction at maximum soil water content (F eff max ) after Gedney and Cox (2003), defined by where ω liq and ω froz are the liquid and frozen soil water content in the top 0.75 m of the upper soil, respectively. The minimum 0.75 m depth is required so as to prevent surface wetland creation during early winter in boreal regions, when the ORCHIDEE soil still contains significant liquid water in the deep layer. In reality, this liquid water gets trapped at the bottom of the soil profile, and cannot contribute to wetland formation. Runoff from the lower layer (drainage) is then parameterized at each time step as the amount of water that cannot infiltrate into the soil column if the total water holding capacity is reached. We also introduce a subgrid variability of evapotranspiration. In partly saturated grid-cells, it is computed as the average of the evapotranspiration depending on the grid-cell average soil moisture weighted by (1-F eff max ) and potential evapotranspiration weighted by F eff max . The surface energy balance is modified in consequence.

Description of the simulations
The global meteorological forcing to drive ORCHIDEE is provided by the monthly NCEP climate forcing data corrected by the high-resolution gridded data sets of the Climatic Research Unit (N. Viovy, personal communication, 2009, http://dods.extra.cea.fr/data/p529viov/cruncep/ readme.htm) on a 6-hourly time step at 1 • resolution. This forcing is a combination of global observation-based datasets with the model reanalysis of the National Center for Environmental Prediction-National Center for Atmospheric Research (NCEP-NCAR). Precipitation is distributed homogeneously over each time step. The present-day vegetation map is prescribed from Loveland et al. (2000). Soil albedo is defined from Zobler (1999). For running ORCHIDEE-TOP, we use the topographic indices at the 1-km resolution given by the HYDRO1k dataset (http://www.edcdaac. usgs.gov/gtopo30/hydro). In each grid-cell, we calculate the mean, standard deviation, and skewness of the HYDRO1k topographic indices and we used them to parameterize a three parameter gamma distribution (Decharme et al., 2006;Sivapalan et al., 1987). Simulations were performed at a 30-min time step over the 1985-2004 period. All simulations start from initial condition of soil moisture at ω max . The 1985-1992 years are used to reach the soil water steady state equilibrium then the 1993-2004 period is used for evaluation against river runoff and satellite data. Three simulations are performed to compare the current version of ORCHIDEE with the new versions ORCHIDEE-FRZ and ORCHIDEE-TOP.

Evaluation against river discharge
The simulated river discharge is compared with monthly observations at gauging stations distributed over the largest river basins globally (see Fig. 2a the spatial distribution of these watersheds), in line with Decharme and Douville (2007). Simulated runoff and drainage on a global grid is converted into river discharge using the routing scheme of ORCHIDEE (Ngo-Duc et al., 2005). In this scheme, at each time step, surface runoff and drainage are temporarily stored in three reservoirs, each of which has different residence time constants. The water is then progressively routed

Fig. 2. (a)
Map of the river basins considered in this study and (b) global breakdown into 11 regions inspired from the TRANSCOM atmospheric transport models intercomparison study (Gurney et al., 2002).
to the oceans in the direction of greatest slope and taking into account the tortuous path of the river channels through a cascade of linear reservoirs along the river network, described at the 0.5 • resolution (Ngo-Duc et al., 2007). Mean annual cycles and monthly discharge anomalies are computed for the years 1985-2000 for which we have observations. Model results are evaluated against observed river discharge using the ratio between simulated and observed annual mean discharges Q sim /Q obs (hereafter named annual discharge ratio criterion), the model efficiency, Eff (see definition in Appendix, Nash and Sutcliff, 1970) as well as correlation (r 2 ) calculated between observed and modelled monthly discharge anomalies (these anomalies are obtained removing a mean seasonal cycle over the period). In addition, the root mean square error diagnostic brings information about observation-model co-variability (as r) but also about the conditional and unconditional model bias (Weglarczyk, 1998). The fact that the model does not include anthropogenic modifications of river hydrographs (e.g. dams over the Ob river, reservoirs, irrigation in the Mississippi basin) could lead to a potential limitation to the use of observed streamflow to evaluate the models.

Evaluation against satellite observations of inundated area
We compare inundated areas inferred from satellite observations , called hereafter P10, with the distribution of F max modelled by ORCHIDEE-TOP. The P10 product is an estimate of the extent of episodic and seasonal inundations, wetlands, rivers, lakes and irrigated agriculture at 0.25 • resolution at the equator. The P10 method uses a complementary suite of satellite observations including passive microwave observations (SSM/I emissivities), active microwave observations (ERS scatterometer), along with AVHRR-NDVI. An unsupervised classification of the three sources of satellite data is performed, and pixels with satellite signatures likely related to inundation are retained. For each inundated pixel, the monthly fractional coverage by open water is obtained using the passive microwave signal and a linear mixture model with end-members calibrated with scatterometer observations to account for the effects of vegetation cover. For the boreal regions, where microwave measurements are sensitive to the snow cover, snow masks were used to edit the results and avoid any confusion with snowcovered pixels. The weekly Northern Hemisphere and Southern Hemisphere snow mask from the National Snow and Ice Data Center (NSIDC) is adopted and averaged on a monthly basis (Armstrong and Brodzik, 2005). The methodology described above was initially applied to the 1993-2000 period  then updated to 2004 . We use here the dataset available at a monthly time scale for 1993-2004. More detailed information concerning the seasonal and interannual behavior of the surface water extent dataset for specific regions can be found in  and Papa et al. ( , 2008 respectively for the Indian subcontinent and the boreal regions. We aggregated the P10 inundated area data to 1 • resolution for comparison with the model outputs. The P10 dataset provides the distribution of inundated fraction (water-logged) whereas the ORCHIDEE-TOP model calculates the subgrid fraction at maximum soil water content. Hence, the two variables are not comparable in absolute value. In the aim to simulate wetland extent compatible with P10, we introduce a global parameterization in order to deduce calibrated wetland fractions (F wet ) from fractions at maximum soil water content (F max ). This parameterization is based on the assumption that the mean topographic index (λ m ) is prone to uncertainty, in particular linked to the resolution of the used digital elevation model (here HYDRO1k), what justifies some calibration (Ducharne, 2009). We performed a shift of the topographic index distribution by modifying the topographic index in all grid-cells: where c is a global constant. This leads to modifying the sub-grid topographic index distribution, called hereafter δ ′ .
In the idealized case in which the equations of the Sect. 2.3 are given, i.e. without the bias correction of Saulnier and Datin (2004), the wetland fraction, F wet , would be computed similarly to Eq. (4): c has been optimized to obtain an annual global F wet close to the global annual P10 +29 %, i.e. P10 + the estimate of drained wetland extent since the pre-industrial period (Sterling and Ducharne, 2008). This leads to a global "pristine" wetland fraction about ∼ 3.2 %. With c (unit in ln(m); see Ducharne, 2009), the yearly global F wet is equal to 3.4 % while the mean annual F max fraction over 1993-2004 is 9.7 %. However, the P10 absolute inundated area fractions themselves are prone to some uncertainties, e.g. the satellite observations have difficulties in catching small, isolated water saturated patches in dry areas, as well as small dry patches in largely wet areas (see also the discussion in Sect. 6). Moreover, the floodplains processes are not simulated by ORCHIDEE-TOP that could lead to mismatch between the two products. All these elements make that the inundated fraction given by P10 and the subgrid fractions given by ORCHIDEE-TOP are difficult to compare in absolute values even after the above-described parameterisation. Hence, in the following, we focus our evaluation on comparing (1) the spatial distributions of P10 and F wet (and F max ) at the scale of large regions with extensive naturally inundated areas, and (2) the seasonal and interannual variability of normalized P10 and F wet .
The natural hydrological unit is the river catchment, but some regions rich in surface water bodies in P10 do not belong to the largest rivers basins of Fig. 2a, e.g. over the northwest of Canada. Therefore, we also consider the approach which consists in dividing the global land area into 11 regions, inspired from the TRANSCOM atmospheric CO 2 transport models intercomparison study (Gurney et al., 2002) (Fig. 2b). This arbitrary breakdown is a reasonable compromise to evaluate continental scale processes that lead to saturated area, and it is compatible with the estimation of CH 4 fluxes from wetlands using atmospheric inversions (Bousquet et al., 2006). In addition, we also analyze individual river catchments in order to study the relationship between F wet and river discharge, including inter-catchment differences.

Results of the comparison with river discharge
The seasonal cycle and interannual discharge variability are shown in Fig. 3 for representative watersheds in boreal regions (the Yenisey, Ob, and Amur), in temperate regions (the Mississippi and Danube) and in tropical regions (the Ganges and Rio-Amazonas). The corresponding values of each statistical criterion are given in Table 2 for both seasonal cycle and year-to-year anomalies. Adding soil water freezing (ORCHIDEE-FRZ) leads to a seasonal peak in runoff whose intensity, except for the Ob river, tends to be in better agreement with observations (Yenisey in Fig. 3 and Lena, not shown). But a too-early runoff peak (one month earlier) is produced, as compared to the observations, and this limits the improvement of Eff over boreal regions. Possible explanations for the early discharge peak include the crude representation of snow, and thus of the snowmelt timing (as mentioned in Decharme and Douville, 2007), and the lack of glacier runoff in the model, both of which could contribute to the error. However, comparison between seasonality of snowmelt in ORCHIDEE and observed snow cover (NSIDC snow cover product provided by ISLSCP-II) (not shown) did not show any systematic bias in the snow-melt timing over the different watersheds. Cox et al. (1999) noted that on coarse spatial scales such as those of typical GCM resolutions, heterogeneities in frozen soil and soil freezing hydraulic conductivity would allow surface runoff from frozen soil surface to infiltrate into the soil elsewhere in the same grid box. To account for this, one option is to introduce a buffer reservoir with a time constant that enables the storage of water from melting snow and its delivery to the soil through delayed infiltration later in the season (Poutou et al., 2004). Concerning the Ob river where the runoff peak is coarsely estimated, the misfit could be related to the nonrepresentation of large floodplains as well as to the effects of dams (Coe, 2000).
The ORCHIDEE-TOP version leads to an important decrease of yearly river discharge (see Q sim /Q obs in Table 2). The coupling with TOPMODEL leads to a decrease in the soil water content and thus the runoff generated by excess water when the soil is at its maximum capacity. The intensity of this decrease varies throughout the year and from one river to another (left-hand panel in Fig. 3). This decrease globally leads to a degraded simulation of the river flow seasonal cycle, except for the Ganges (Table 2). This suggests that the decrease of the soil water content following the representation of fractions at maximum soil water content is too large. This could be attributed to the difficulty to simulate the absolute value of the contributing area with a two-layer soil model (see discussion in Sect. 6). Regarding interannual variability (right-hand panel in Fig. 3 and Table 2), the combined effect of the added frozen soil processes and subgrid water redistribution leads to an increase in r 2 and/or decrease in RMSE for many rivers (Yenisey, Amur, Danube, Mississippi, Ganges). In the Ob watershed, the decrease in the ability to simulate the year-to-year variability in the river flow is totally explained by the freeze/thaw processes (Table 2); the evolution in the permafrost melting or in the infiltration due to discontinuous permafrost could be explain a part of the disagreement. Only for the Rio-Amazonas, accounting for  Table 2. subgrid water distribution leads to decrease both the r 2 and the RMSE.
For the Mississippi at Vicksburg (Fig. 3), the results of the initial version of ORCHIDEE show a large soil water content during summer (not shown) that leads to an excess of river flow during the period July-October. Guimberteau et al. (2010) attributed this river flow overestimation to a too weak evapotranspiration and corrected it by increasing root depth of the present vegetation (natural and anthropogenic PFT "C3 grass"). Part of the difference in river flow between observations and our model could also be attributed to anthropogenic impacts, which cannot be taken into account in our approach. For instance, irrigation could modify the hydrological cycle too (10 % of the irrigated agriculture over the world is in the USA; Siebert et al., 2005). When coupling with TOPMODEL, this large summer soil water content leads to substantial contributing areas then Dunne runoff resulting in a negative feedback on the soil water content (i.e. Table 2. Skill score for simulated discharges by ORCHIDEE (black), ORCHIDEE-FRZ (green) and ORCHIDEE-TOP (red) over the 7 rivers given in Fig. 3. The ratio between simulated and observed annual mean discharges (Q sim /Q obs ) and the monthly efficiencies (Eff) (calculated over the whole discharge observation period given in Fig. 4, right panel) are given as well as the correlation (r) and RMSE (in m 3 s −1 ) between simulated and observed monthly anomalies. leading to a decrease in the soil water content). It seems that this decrease is too large (see Q sim /Q obs in Table 2) as mentioned above.
For the Rio-Amazonas at the Obidos gauging-station located ∼ 800 km upstream from the river mouth in the state of Para, Brazil, the application of ORCHIDEE-TOP shows a degradation compared to the current version (Fig. 4). Both ORCHIDEE versions have a too early peak as compared with observations. This basin is prone to large scale floodplain inundation that leads to the formation of a buffer reservoir during the wet season with its own time constant linked up to re-infiltration rate and aquifer recharge, which can explain the observed lag between precipitation and flow maximum. It seems necessary to account for this reservoir as in Decharme et al. (2008) to improve the Amazon discharge simulation. Another model shortcoming could also be an underestimation of the total soil depth over the Amazon basin as suggested by Kleidon and Heimann (1998).

Spatial distribution
In this section, we evaluate ORCHIDEE-TOP results with respects to remote sensing observations. Figure 4a compares the spatial distribution of the yearly mean fractional inundation from satellite (P10) and the mean annual modelled wetland fraction (F wet ). Figure 4b gives the annual maximum fraction for the same variables (P10, F wet ), all averaged over 12 yr regardless of the month of maximum. Given the caveats of comparing the absolute values of these different quantities, as explained in Sect. 3.2.2, we also ranked and separated each grid-cell into deciles in Fig. 4c. This allows us to compare the spatial distribution of the grid-cells that contribute most to the global signal. Maps are also given for F max (middle column). Table 3 also gives the Spearman correlation between the inundated fractions and F max (F wet ) at the annual maximum of all the grid-cells into each region defined in Fig. 2b and at global scale. This non-parametric correlation is preferred here to the Pearson correlation over the direct values of modelled fractions and P10 because it is compatible with the relative values, as expressed by percentiles. The histograms of P10 and F wet spatial distributions at the annual maximum are given in Fig. 4d. Figure 5 shows the relative contribution (in percent) of each region of the TRANSCOM breakdown to the total area of P10 or F wet , for the yearly mean (Fig. 5a), for the mean annual maximum regardless of month of maximum ( Fig. 5b), for January (Fig. 5c), and July (Fig. 5d) over 1993-2004. Globally, the distribution of P10 in Fig. 4a, b and c shows key regions: the Ob plain, the Ganges plain in northern India and the regions of north-eastern Canada. This results in a high spatial variance (Fig. 5) between the different regions of the TRANSCOM breakdown at any time periods of the year (annual, January or July). The distribution of F wet in Fig. 4a shows more key regions that have a large contribution to the global maximum saturated area than P10. Regions that contribute to global F wet are more evenly distributed than for global P10 in Fig. 5. Therefore, F wet has a lower betweenregions spatial variability than P10 (std P10 > std F wet ). For instance during July, the two biggest contributing regions to Table 3. Spearman correlation (ρ P10,F max ) between the inundated (P10) and fractions at maximum soil water content (F max ) at the annual maximum (regardless of the month of maximum) of the gridcells contained in each region defined in Fig. 2b and  P10 shows significant fine-grained structure with gridcells belonging to the highest percentiles-class (Fig. 4a, b and c). These structures are rather homogeneously distributed, and affect the river catchments of Volga, Mississippi and Amazon (Fig. 4a, b and c). The modelled F max distribution does not capture these fine-grained patterns (Fig. 4a, b and  c). As a result, the histograms of P10 and F max spatial distributions look different (Fig. 4d). The grid-cells between the 0th and 30th P10 percentiles have an inundated fraction equal to zero (Fig. 4d). The median value of the distribution is close to the mean for F max , but not for P10, because the P10 field is more skewed towards low values (Fig. 4d). These differences could be linked to local differences in soil water holding capacity and water routing that are not resolved in ORCHIDEE-TOP, as well as sub-grid rainfall distributions, such as preferential convective rainfall over mountain regions, that may covary with the sub-grid distribution of flooded areas. The introduced parameterization leads to a fine-grained structure ( Fig. 4a and b) and to a distribution (Fig. 4d) of the subgrid fraction more consistent with P10. This suggests that a part of the differences between P10 and F max could be linked to intrinsic differences between inundation and "full bucket" state of a soil in the model. For instance, inundation has more complex mechanisms than field capacity and saturation such as irrigation, damming, formation of ponds, flooding of impermeable soils, thus being spatially more heterogeneous over the globe. This can explain also the difference between regions observed for the Spearman correlation (Table 3). Figure 5 shows differences in the contribution of some regions to global P10 and F wet . Tropical Asia has a bigger contribution to P10 than to F wet (Figs. 4a and 5). This region    Fig. 2 for the regions definition) to the global inundated/saturated area, both for the yearly mean (a), for annual maximum (b), mean January (c) and mean July (d) over 1993-2004. For each part of the Fig. 5, the regions are ranged in descending order of contribution to global inundation (P10). The Pearson (Spearman) correlation between the global distribution among the different regions of P10 and F wet , r P10,F wet , (ρ P10,F wet ), is also given for each panel (a, b, c, and d).
The spatial standard deviation among the different regions for P10 and F wet (std P 10 and std Fwet respectively) is indicated for (a, b, c) and (d).
is strongly influenced by irrigation even after removing rice paddy areas, which is not accounted for in ORCHIDEE-TOP. Accounting for spatial variability in the soil water holding capacity (RU) into ORCHIDEE-TOP provides small holding capacities in this region that seems to prevent highenough F wet (Fig. A1). In contrast, grid-cells in the eastern United States have larger relative contributions to F wet than to P10 (Fig. 4a). This may be because, in our modelling approach, we do not account for human activities, especially drainage, that may have lead to an estimated loss of 53 % of the original wetlands in the United States (i.e. 200 yr ago) (Mistch and Gosselink, 2000). The greatest historical Papa et al. data  wetland losses occurred in the lower Mississippi alluvial plain and the prairie pothole region of the north central states (Mistch and Gosselink, 2000), just where F wet is found to be larger than the inundation extents in P10. Finally, we can see higher local F wet contributions than P10 over tropical South America in Fig. 4a. In this region, the flooding mechanism of floodplains and the routing of water are not accounted for in ORCHIDEE-TOP. The Spearman correlations between regional averages of P10 and F max (not shown) are close to the correlations between P10 and F wet given in Fig. 5. This means that the introduced parameterization to deduce F wet from F max has a small effect on the contribution of the different large regions to the global signal. However, it has a more significant effect on the distribution over the grid-cells in each TRANSCOM region: in the Table 3, the ρ P10,F wet correlation is lower than ρ P10,F max for most regions.
The contribution of each region to global F wet and P10 varies through the year (Fig. 5). The between-regions correlation r P10,F wet is higher for the different months and for the yearly mean (0.76, 0.80 and 0.82 for respectively the yearly mean, January and July) than for the yearly max (0.59). The between-region difference in spatial variance between P10 and F max is larger in July than in January, that is | std P10std F wet | July > | std P10 -std F wet | January . In January when soils are frozen in the Northern Hemisphere, the contribution of boreal Europe, Siberia and North American regions to the total is negligible both for F wet and P10. This, and the fact that this is the dry season over Southern Asia and the Amazon, makes it easier to reproduce the relative contribution of each region to the total, hence the large value of r P10,F wet in Fig. 5.

Seasonal cycle
In this section we compare satellite P10 and simulated F wet seasonal variability. Figure 6 shows the distribution of the most frequent month of maximum P10 and F wet occurrence (Fig. 6a) and the lag in months that produced a maximum correlation of both P10 and F wet with monthly precipitation (Fig. 6b). We find a good agreement between the timing of maximum values of P10 and F wet (Fig. 6a). Indeed, the latitudinal gradient is well reproduced, except for tropical Asia where maximum F wet lags maximum P10 by one month. For the northern Sahel region and central western Africa, F wet also reaches its maximum after the one observed in the P10 satellite data. There is also good agreement between F wet and observed P10 regarding the latitudinal gradient of the lag between inundation and precipitation (Fig. 6b). At high northern latitudes, this lag is likely explained by the lag between snowfall and snowmelt. We have large difficulties in representing fine-scale structures seen in the lag of P10 over the Amazon, North America and temperate Eurasian region, where the link between precipitation and F wet is stronger than the one observed between precipitation and P10.
The mean seasonal cycle for each region of the TRANSCOM breakdown (Fig. 1) and for the globe is shown in Fig. 7. The seasonality of snowfall and rainfall is also shown. Each curve is normalized to unity by dividing each  Fig. 7. Mean seasonal cycle for each region and at global scale for P10 and F wet . Seasonality of snowfall and rainfall (both in relative way compared to snowfall+rainfall) are also added. All the curves are divided by their maximum value to obtain the same upper limit for all the curves. The FMT between P10 and F wet (see Appendix for FMT definition) is given in each region's plot.
value by the yearly maximum averaged over the whole period except for the seasonal cycle of snowfall and the one of rainfall which are divided by the maximum of total precipitation (snowfall + rainfall). Precipitation, either in the form of rainfall or snow, is a forcing of ORCHIDEE-TOP. The agreement between the seasonal cycle of P10 and F wet is expressed by the value of the figure of merit in time (FMT) (Hourdin et al., 1999;Krinner et al., 2005 and see Appendix) given in each sub-panel of Fig. 7. Overall, we obtain relatively high FMT values, yet with some discrepancies among regions. The FMT goes from 64.2 % for Northern Europe up to 92.4 % for tropical Asia. Despite the probable influence of the rice paddies flooding and/or irrigation on the wetland extents in tropical Asia, the seasonal cycle is well captured by ORCHIDEE-TOP. The FMT at global scale is at the bottom range of regional values, which is explained by the discrepancy between the contribution of each different region in global P10 and F wet discussed above.
Over South America (tropical and temperate), we obtain a high value of FMT, due to the low amplitude of the seasonal cycle over these regions for both P10 and F wet even though there is a discrepancy between the timing of maximum P10 and the one of F wet , as seen in Fig. 7. In tropical South America, there is a lag of few months between P10 and rainfall, and a weak correlation between these two variables (r = 0.03). The correlation between P10 and rainfall increases to 0.60 at a 3 months lag. Using ORCHIDEE-TOP, we capture this lag between F wet and rainfall, but it is underestimated for tropical South America and over-estimated for temperate South America. The underlying mechanisms are complex. Over the Parana region of temperate South America, inundation seems to be driven by upstream precipitation events . Over the Amazon, the water residence time in floodplains and aquifer recharge may also play a role (Decharme et al., 2008).
Over boreal North America, northern Europe and boreal Eurasia, the match between the seasonal cycle of P10 and F wet is good, despite lower FMT values than for tropical regions due to larger seasonal amplitude (Fig. 7). The seasonal cycle is longer for F wet than for P10 (Fig. 7). The observed onset of the seasonal increase in satellite P10 is well captured by modelled F wet for boreal Eurasia, but it occurs ∼ 1 month too early for boreal North America and Europe. Year-to-year variability of P10, F max and of snowfall + rainfall. Each curve is obtained using moving average over 12 month on monthly anomalies. These anomalies were obtained by removing the mean seasonal cycle over 1993-2004 then dividing by the maximum value. Correlation between P10 and F wet is given between brackets in each region's plot.
The termination of the F wet seasonal cycle when soils freeze in autumn occurs one month too late in all three boreal regions. Some of the seasonal cycle phase differences between F wet and P10 could be attributed to the poor resolution of vertical water distribution in ORCHIDEE-TOP that prevents a realistic simulation of the liquid soil water content. Moreover, flooding is one of the different processes that can lead to inundation (P10) which is not accounted for in ORCHIDEE-TOP. The thickness criteria used to compute F eff wet from F wet (75 cm, see Eq. 6) have little influence on the length of the seasonal cycle of F wet , e.g. FMT changes to 2 % for boreal Eurasia if the depth criteria changes from 75 cm to 1.2 m. The FMT between P10 and F wet is higher than one between P10 and liquid soil water content (not given) precisely due to this depth criterion that excludes the contribution of liquid water located in the deep soil layers to the wetland formation. The high FMT (not shown) between rainfall (blue curve in Fig. 7) and P10 seasonal cycle (black curve in Fig. 7) for boreal regions is attributed to the fact that maximum rainfall occurs in summer, when the soils are not frozen at the surface, hence preventing complex lag effects of P10 to rainfall. Figure 8 compares the interannual variability of satellite P10, modelled F wet and snowfall + rainfall forcing. A 12 month moving average is applied to the monthly anomaly time series to remove the seasonal cycle. Anomalies are normalized first by dividing the time series of each region by its maximum value, and then removing the mean seasonal cycle over 1993-2004. The interannual correlation between P10 and F wet over the 12-yr period is given in each sub-panel of Fig. 8. There is a large variability for the inter-annual correlation (r) between P10 and F wet among the regions with r going from −0.30 in tropical Asia to 0.85 in temperate Eurasia. Except for boreal regions where snow can create a delay between precipitation and water availability to saturate soils, the interannual correlation between F wet and precipitation is always high (> 0.54) in all regions. This is observed even in the tropical and temperate regions where correlation between F wet and P10 is poor (< 0.1). For instance, the correlation between F wet and precipitation is equal to 0.81 over tropical Asia and to 0.944 over Southern Africa. Simulated F wet year-to-year variability is thus close to the one of precipitation variability, whereas the variability of P10 seems to be explained, at least in some regions, by other factors.

Interannual variability
Note that the P10 data show a decrease in P10 in some regions (temperate South America, North America, etc.) as well as at the global scale at the end of the 1990's. P10 captures well this decrease for regions where a coincident decrease of precipitation is also observed (North America), associated with the drought there (Hoerling and Kumar, 2003). However, in other regions (temperate South America for instance) and at the global scale, the decrease observed in P10 is not associated with a decrease in precipitation. We assume that such changes in P10 could be attributed in part to anthropogenic influence. But, given the lack of independent observations at large scale and over many years to evaluate P10, more investigations are needed to draw such conclusions. Therefore, at global scale, as for the seasonal cycle, the poor ability of ORCHIDEE-TOP to match the year-toyear variability of P10 (r = −0.47) can be explained, at least partially, by some discrepancy in the contribution of the different regions (cf. Sect. 5.1).

Discharge-to-flooding relationship for Siberian rivers
In the previous sections, we evaluated independently ORCHIDEE-TOP against both river discharge (an indirect validation of the entire hydrology) and satellite observations of flooded area (a more direct validation of saturation, despite differences between inundation and saturation). The seasonality of river flow and inundation are not independent. Hence, we focus the first point of this discussion on our potential ability to capture the relationship between river streamflow and inundation extent at the basin scale as well as the variability of this relationship from one basin to other. This will allow us to assess if the conclusions obtained through the two independent cross-validation datasets of discharge and satellite flooded area are consistent with each other. We focus on three large Siberian river basins where P10 was extensively evaluated against in-situ river discharges, as well as in-situ and satellite-derived snow water equivalent products. In particular, the relationships between snow, inundation, and runoff were studied by Papa et al. ( , 2008, showing significantly different behaviors between the Ob, Yenisey, and Lena basins. In order to examine the capacity of ORCHIDEE-TOP to reproduce the regional relationships between inundation and river flow over Siberian watersheds, the monthly climatology  of inundation extent and discharge is shown in Fig. 9 (left) for both observations and simulations. Concerning the Yenisey basin that shows a large downstream/upstream gradient in inundation extent, we consider only the inundated/saturated area downstream of the gauging station. All the curves are normalized by their maxima. For each basin, we analyse the relationship between the river discharge on the one hand, and the satellite-derived inundation and simulated saturation fractions on the other hand ( Fig. 9, right). Papa et al. (2008) observed an hysteresis between the seasonality of the discharge and flood extent, with variable degrees from one basin to other. While there is a large increase of both inundation and discharge when snow melting season starts over these three basins, the relationship between discharge and inundated area during summer is different. Contrary to the Ob, the Yenisey and the Lena basins present a sharp decrease in the discharge while flood extent only slowly declines. This could be explained by an interruption of the connection between the river channel and the surrounding floodplains during the flood recession from July. Water that is ponding in the lowlands and is disconnected from the river channel only slowly disappears through evaporation and percolation during the summer period. The situation is different for the Ob River that has a large hydraulic network with high storage capacity (especially in the northern part) due to floodplains and less extensive permafrost coverage than the other two Russian watersheds. Figure 9 shows that ORCHIDEE-TOP is able to capture one important trait of the observed relationship between discharge and inundated area for Lena and Yenisey, namely the lag between recession of inundated area and discharge decrease in June. For both basins, we simulate too much delay between the maxima of saturation area and discharge, in comparison with the observed values. The timing of the flooding recession in F wet area occurs also too late by ∼ 1 month in comparison with observations. This disagreement in the floodingto-discharge relationship can be explained, as mentioned before, by the intrinsic difference between what is observed (inundation) and what is modelled (wetland fraction). The flooding-to-discharge relationship is less well-captured for the Ob basin, however. The simulated discharge increase and inundated area decrease occur too early and too late, respectively (Fig. 10). That can be explained by the fact floodplains and delta formations are not simulated by the model or by a crude representation of snow, and thus of the snowmelt timing.  also explained the lags between maximum P10 and river discharge for the Ob River by ice jams in the river valley that contribute to delay between the peak of the discharge and the flooding maximum.

Suitability of TOPMODEL to simulate wetlands for CH 4 emissions modelling
Several studies used subgrid topographic information (Coe, 1998;Krinner, 2003) or TOPMODEL's concepts (Gedney and Cox, 2003) to diagnose wetland extents at global scale and their time variability. We focus the second point of this discussion on the suitability of TOPMODEL's concepts to simulate wetland areas at large spatial scale for modelling the associated CH 4 emissions. Because wetlands have a considerable range of hydrologic conditions and because of their great variation in size, location, seasonality and human impact (Mistch and Gosselink, 2000), their definition is controversial (Reichhardt, 1995). Nevertheless, they have many specific features, the most important being the presence of standing water for some period during the growing season either at the surface or within the root zone. The presence of unique chemical soil properties and organisms, especially vegetation, are consequences of saturated soils and resulting anaerobic conditions (Mistch and Gosselink, 2000). At the global scale, saturation seems mainly linked to geomorphology and climate (Mistch and Gosselink, 2000), corresponding well with the assumption used in TOPMODEL. Contrary to early subgrid parameterizations proposed by Coe (1998) and Krinner (2003), TOP-MODEL does not need prescribed residence time of water and pre-dimensioned reservoirs. As a consequence, it offers the possibility of a more mechanistic representation of wetland.

Difficulty to simulate reasonable global wetland coverage
The coupling between ORCHIDEE and TOPMODEL as described in this study allows diagnosing the fraction at maximum soil water content. The simulated fractions are much larger than from the P10 dataset , which provides the distribution of inundated fraction (waterlogged). These two variables are not comparable in absolute value. This problem does not arise from the coupling with TOPMODEL but underlines the necessity to have a more process-based model for the soil hydrology to simulate physically the wetland dynamics. This limitation is not linked to the size of the bucket but more to the physics used to represent the water fluxes in the soils. A two-layer bucket does not allow estimating the mean deficit to the saturation over each grid-cell. This problem in the absolute value of the contributing area explains also the degradation of the Q sim /Q obs criteria for river flow after coupling with TOPMODEL. This underlines the difficulty to simulate the absolute values of contributing areas and wetlands with a simple two-layer bucket even if it is coupling with a TOPMODEL approach. A possibility to improve the representation of the global water budget could be also to "play" with the computation of the evapotranspiration (ET) on the F max . The ET computation is not modified depending on the TOPMODEL/LSM coupling in Decharme et al. (2006) while it is the case in this study (see Sect. 2.3) and in Koster et al. (2000) and . The modification of ET computation feeds back on the soil water content and could lead to change both the river flow magnitude (Q sim /Q obs criteria) and the absolute value of F max . The mean ET over a given grid-cell could also be modified only depending on F wet and not on F max , in a similar way as Koster et al. (2000) and . These problems underline again the difficulties to correctly simulate the absolute value of the saturated/wetland areas (which contributes to Dunne runoff and on which ET is close to evapotranspiration) with a two-layer bucket model. However, both the spatial distribution of the wetland areas and its time variability are rather well-captured.
It seems that the coupling of TOPMODEL with a more process-based LSM (MOSES) (Gedney and Cox, 2003) allowing to reach the saturation leads also to saturated areas that are much larger than the accepted wetland inventories (Aselmann and Crutzen, 1989;Matthews and Fung, 1987;and now Prigent et al., 2007). As for F max , modelled saturated areas and inundated areas as given by Prigent et al. (2007) are not comparable in absolute value. Saturated areas do not necessarily appear to satellites as free-water inundated area, making saturated areas smaller than inundated ones. Inundated area could locally be larger than saturated areas, if flooding corresponds to the presence of water on impermeable soil with no evident link to a saturated soil in depth. Moreover, the absolute inundated area fraction in the dataset is prone to uncertainties ) (see below).
If these saturated areas simulated by a TOPMODEL/LSM are combined with wetland emissions models based on flux sites measurement (e.g. Walter et al., 2001), this will lead to overestimated global wetland CH 4 emissions as compared to the accepted range of the contribution of wetland to global sources (Denman et al., 2007). This raises the question of the required hydrologic conditions (saturation versus inundation) to have wetlands that emit CH 4 . To solve this issue, Gedney and Cox (2003) optimized a maximum critical topographic index to match the wetland inventory of Aselmann and Crutzen (1989). They assume that wetlands could be considered as areas of stagnant water, excluding all other areas where the water table could rise above the surface (i.e. where the local deficit > 0) and results in a significant flow. This strategy is possible because they used the classical approach of Beven and Kirkby (1979) where the local deficit could be negative. However, Saulnier and Datin (2004) have underlined an approximation in the Beven and Kirkby (1979) approach and suggested a formulation to correct the resulting bias. This corrected formulation is adopted in the present study: it consists in limiting the local deficit to positive values (see Appendix 1 of Decharme et al., 2006). Thus, the strategy of Gedney and Cox (2003) and the correct formulation of Saulnier and Datin (2004) are not consistent. Kleinen et al. (2012) also introduces thresholds on the mean topographic index. This allows reproducing fine-grained structure suggested in the Sect. 5.1.
Another solution adopted by Ringeval et al. (2011) is to calculate anomalies from the field capacity area given by TOPMODEL relatively to the Prigent et al. (2007) data. This approach is not totally satisfying either due to different natures of the two products as explained before. In the present study, a parameterization consisting in a shift of the topographic index distribution is introduced (see Sect. 3.2.2) to match the global wetland coverage. This leads to more consistent spatial patterns with respect to the data. It is based on the assumption that the mean topographic index is prone to uncertainties (Ducharne, 2009) and hence could appear as a good alternative to the parameterizations introduced in Kleinen et al. (2012) and Gedney and Cox (2003). Such a parameterization could also be tested with a more process based soil water model as in Gedney and Cox (2003). Moreover, the involved parameter would be optimized at regional scale (instead of global as in the present study) and would allow both to match in a better way the data and to account for variability in the drained wetland area between the different world regions (Sterling and Ducharne, 2008). In all the cases, these types of parameterizations have small effects on the simulated time variability.
Finally, a better coupling between carbon and hydrology could be an alternative way to better constrain the wetland extent. In fact, soil carbon accumulates under anaerobic conditions and is necessary as a substrate for methanogenesis bacteria. The combination of subgrid topography used in TOPMODEL and subgrid soil carbon (e.g. IGBP-DIS data at 5 min resolution, Global Soil Data Task Group, 2000) could allow limiting spatially the wetland extent to the area with high soil carbon content.

Evaluation of modelled wetland extents
Even if multiple satellite approaches as in  could be able to circumvent the problems mentioned by Frey and Smith (2007), they still have limitations. For instance, grid-cells with high soil carbon content representative to peatlands are not well-captured in P10 (R. Spahni, personal communication, 2010): peatlands are not necessarily free-surface waters and as consequence are not present in P10. Moreover, as mentioned above, the absolute inundated area fraction in the P10 dataset is prone to uncertainties . Satellite observations have difficulties in catching small, isolated water-saturated patches in largely dry areas, as well as small dry patches in largely wet areas. Along the coast, satellite products are also contaminated by ocean signals. Thus, regional inventories as mentioned in Peregon et al. (2008) closer to ground-truth should also be considered. Or, as done in the present study, the evaluation has to focus on the spatial and time variability.
Finally, evaluation of wetland extent should be coupled with the evaluation of wetland CH 4 emissions. Year-to-year variability in wetland extent seems to explain a large part of the variability in wetland CH 4 emissions . Thus, the evaluation of the variability in modelled emissions obtained using a TOPMODEL/LSM coupling is also a way to evaluate the wetland dynamic. Evaluation of year-to-year variability in global CH 4 emissions could be done against top-down approach results (Bousquet et al., 2011;Ringeval et al., 2011) but are themselves prone to lots of uncertainties (Denman et al., 2007).

Accounting for the water table depth
Improving the modelling of wetland CH 4 emissions using TOPMODEL/LSM is achieved through the representation of the water table depth (WTD). Current couplings between TOPMODEL and an LSM are restricted to simulation of areas where the WTD is at or above the soil surface. However, a wetland can emit CH 4 even if the WTD is under the soil surface.
The WTD is a key-variable by delimiting the extension of the anoxic soil zone (where CH 4 is produced) and of the overlying oxic one (where CH 4 is oxidized by methanotrophy) (Walter and Heimann, 2000). The WTD seems to act as an on-off switch (Christensen et al., 2003).
Although uncertainties remain about the sensibility of CH 4 emissions to WTD (Christensen et al., 2003;Updegraff et al., 2001), the value of the WTD at which no emissions occurs seems below the soil surface. Thus, it is necessary to account for both vertical and horizontal hydrological variability as underlined by Bohn et al. (2010). . Spatial variability in soil water holding capacity introduced in the initial version of ORCHIDEE (RU). Spatial soil water holding capacity computed using both soil texture and soil organic content (b) is used in all ORCHIDEE versions described in this paper (Fig. 1). Spatial soil water holding capacity computed using only soil texture is given in (a) for information. Improvements could be brought by associating not only a topographic index with saturation (λ sat , see Appendix 1 of Decharme et al., 2006) but also some topographic index with different soil water deficit into each grid-cell. Nevertheless, such improvements brought to Saulnier and Datin (2004) formalism could not give information about the value of the WTD as soon as the WTD is above the soil surface (cf. Sect. 6.2.1). The thickness of the water layer above the soil surface determines yet the CH 4 oxidation intensity before it reaches the atmosphere (Walter and Heimann, 2000) and can be used to discriminate between lakes and wetlands (Coe, 1998;Krinner, 2003;de Noblet-Ducoudré et al., 2002).

Accounting for wetland diversity
Finally, a global approach such as TOPMODEL does not give the possibility to model the hydrologic diversity in ecosystems covered by the term "wetland". First, TOP-MODEL accounts only for wetlands developed from saturated soils by beneath and not from floodplains mechanisms or ombrotrophic bogs (i.e. wetland receiving water exclusively from precipitation and not influenced by groundwater). Then, even among the class of wetland developed from saturation, the water recharge could be different and TOP-MODEL cannot account for this diversity. The local recharge processes will partly determine the wetland productivity (ombrotrophic, connected fens/unconnected bogs) (Mistch and Gosselink, 2000), the substrate amount for methanogenesis, and finally CH 4 fluxes (Updegraff et al., 2001). Different approaches are exploring the possibility to include a statistical representation of the grid cell micro-topography within TOPMODEL and analyse its impact on the simulated saturated area. Also, addition of several PFTs could allow accounting for this diversity (Wania et al., 2009).

Conclusions
This study shows the impact of a soil freeze/thaw and TOP-MODEL sub-grid scale parameterizations on global hydrological simulations performed with the ORCHIDEE LSM. First a classical comparison between modelled river discharge and observations at gauging station is conducted. Globally, the quality of the simulated discharge increases with accounting for freeze/thaw processes but is degraded after coupling with TOPMODEL in relationship to the difficulty to estimate the absolute value of the saturated areas. We suggest that over specific basins (Ob, Rio-Amazonas), the improvement is further limited because the LSM does not simulate floodplains. Second, an original evaluation of the coupling between TOPMODEL approach and a LSM is presented using remote sensing data . The satellite observations provide an estimate of the inundated extent, whereas the model diagnoses the areas at maximum soil water content (F max ). As a consequence, the model results cannot be validated in the absolute sense, but the consistency of the spatial and temporal variability of the two related estimates (inundation areas and F max ) is carefully evaluated. A parameterization is also introduced to estimate the wetland fractions from F max . Despite some difficulties in matching exact locations of individual inundated events, a good agreement is obtained in the spatial distribution of the inundation/saturation areas, into pre-defined world regions. The model reproduces the seasonality of the inundation correctly, but interannual variability is difficult to simulate, especially the decreasing trend observed in the inundated area in the late 90's. The two evaluations (against river flow and remote-sensing data) are complementary, imposing strong constraints on the simulation at basin scale.