Modeling Land Surface Fluxes from Uncertain Rainfall: A Case Study in the Sahel with Field-Driven Stochastic Rainfields

Abstract

In distributed land surface modeling (LSM) studies, uncertainty in the rainfields that are used to force models is a major source of error in predicted land surface response variables. This is particularly true for applications in the African Sahel region, where weak knowledge of highly time/space-variable convective rainfall in a poorly monitored region is a considerable obstacle to such developments. In this study, we used a field-based stochastic rainfield generator to analyze the propagation of the rainfall uncertainty through a distributed land surface model simulating water and energy fluxes in Sahelian ecosystems. Ensemble time/space rainfields were generated from field observations of the local AMMA-CATCH-Niger recording raingauge network. The rainfields were then used to force the SEtHyS-Savannah LSM, yielding an ensemble of time/space simulated fluxes. Through informative graphical representations and innovative diagnosis metrics, these outputs were analyzed to separate the different components of flux variability, among which was the uncertainty represented by ensemble-wise variability. Scale dependence was analyzed for each flux type in the water and energy budgets, producing a comprehensive picture of uncertainty propagation for the various flux types, with its relationship to intrinsic space/time flux variability. The study was performed over a 2530 km² domain over six months, covering an entire monsoon season and the subsequent dry-down, using a kilometer/daily base resolution of analysis. The newly introduced dimensionless uncertainty measure, called the uncertainty coefficient, proved to be more effective in describing uncertainty patterns and relationships than a more classical measure based on variance fractions. Results show a clear scaling relationship in uncertainty coefficients between rainfall and the dependent fluxes, specific to each flux type. These results suggest a higher sensitivity to rainfall uncertainty for hydrological than for agro-ecological or meteorological applications, even though eddy fluxes do receive a substantial part of that source uncertainty.

1. Introduction

In the populated and expanding semiarid regions of tropical Africa, such as the Sahel, vulnerable ecosystems and livelihoods are expected to be exposed to increasingly frequent and intense climate extremes [1]. Yet our knowledge of how they will respond to this climate stress is poor, making these regions hot spots for global change studies [2,3]. Monitoring soil and vegetation conditions, including soil–vegetation–atmosphere exchanges, is especially important in the Sahel [4,5,6] to (i) evaluate crop conditions and future production, to better anticipate and prevent food shortages, (ii) estimate space–time distributions of mass and energy transfers from the land surface to the atmosphere, to improve weather and climate predictions, (iii) improve the reliability of water resources and risk forecasts, (iv) help improve performance of irrigation systems, where they exist, and (v) understand the dynamics of fragile ecosystems. These are key issues in a region characterized by a semiarid, monsoonal climate [7], with convective rainfall displaying strong variability over a wide range of time and space scales [8], a high sensitivity of monsoon processes to surface conditions [9,10,11,12,13,14], and essentially rainfed agriculture [15]. Documenting energy and water budgets and their variability at various scales is crucial if we are to improve our understanding of the West African monsoon system and its intricate relationships with local ecosystems and societies [16].

Due in particular to the strong contrasts in Sahelian rainfall, surface exchanges and conditions are subject to rapid changes in time and space. This makes them impossible to estimate at any significant scale from measurements alone and calls for the development of reliable land surface models (LSMs) that can be used in combination with ground and space data. Mostly for lack of appropriate parameterization and evaluation data, few specific land surface modeling studies have focused on West African ecosystems. Some surface–vegetation–atmosphere transfer schemes were tested locally at single land-cover sites (e.g., [17,18,19,20,21,22,23]). The African Monsoon Multidisciplinary Analyses (AMMA) land model intercomparison project (ALMIP) intercompared a set of global LSMs applied at a regional scale [5] or over contrasting mesoscale sites of the AMMA-CATCH regional observatory [24,25] (please see Appendix A for definitions of abbreviations and acronyms). The Suivi de l’Etat Hydrique des Sols (SEtHyS)-Savannah model, which is used hereafter, was successfully adapted, calibrated, and validated for the local land-cover types from the AMMA-CATCH-Niger (ACN) Sahelian site near Niamey in South-West Niger [26], and was then gridded at the mesoscale over a 2530 km² domain [27].

For model predictions to be useful, their uncertainties and possible biases must be investigated. Given the importance of rainfall forcing as a driver of surface processes in the water-limited Sahel, and the difficulty in reliably estimating its time–space distribution, a critical issue is an evaluation of uncertainties and biases in the estimated rainfields and their propagation through the LSMs. Significant work has been done elsewhere in this regard with rainfall–runoff models, showing for instance that rainfall uncertainties can often, and under various conditions, largely outweigh other sources of uncertainty in modeled streamflow (e.g., [28,29]). The spatiotemporal representation of rainfall is pivotal for flood runoff estimation under various conditions (e.g., [30,31,32,33,34]). More contrasting results regarding error propagation from precipitation to soil moisture were obtained with LSMs for Oklahoma [35,36]. In Northern Germany, Bormann [37] found the highest sensitivity to rainfall resolution for LSM outputs at monthly timescale, but warned against directly transferring their conclusions to other model types and climate.

Rainfall uncertainty first depends on the data source(s) used to derive the space–time distributions of precipitation intensity, at a resolution compatible with the land processes of interest. Several studies evaluated various global, satellite-based, or model-derived rainfall products indirectly through their hydrologic predictive capabilities (e.g., [35,38,39,40]. Pierre et al. [41] compared such satellite products and their effects when applied to a seasonal vegetation model over the Sahelian belt. Nijssen and Lettenmaier [42] and Nikolopoulos et al. [43] investigated the effects of scale on the propagation of modeled satellite rainfall error through hydrological models. Global rainfall products appear quite promising especially for applications at large scales, but they still significantly lack the accuracy and resolution needed for smaller-scale studies, especially under convective tropical rainfall conditions. Radar-derived rainfall is attractive especially for flash-flood forecasting (e.g., [44,45]), but cannot be achieved over vast portions of the globe, including nearly all of Africa. Relatively little effort seems to have been put into evaluating the impact of gauge-based rainfall uncertainty/error or estimation methods on the full, water/energy cycling response of LSMs, and virtually none under the tropical, West African type of conditions. For a reliable global uncertainty analysis, the structure of forcing uncertainty needs to be expressed precisely and completely.

In this paper, we provide an analysis of the propagation of uncertainty from high-resolution time–space distributions of rainfall inputs to the dependent surface fluxes of water and energy, through the SEtHyS-Savannah LSM applied over a mesoscale domain in South-Western Niger. Rain distribution is derived from a dense raingauge network using a stochastic rainfield generator calibrated from this data [46]. SEtHyS-Savannah is a spatially-distributed column model, with multiple sub-cell land-use/ecosystem simulation. The study was performed over a 2530 km² window within the ACN area, and through a 6-month period covering a Sahelian monsoon season and the subsequent dry-down. Effects on the main water and energy output variables were analyzed in terms of their predictive uncertainty, expressed through a 50-member ensemble driven by the rainfield generator. As scales of interest may strongly vary with application fields (e.g., fine spatial resolution for agronomic studies, but fine time resolution for weather/climate applications), a range of spatial and temporal scales were explored in output analysis, from 1 day and 1 km² up to the entire space/time domain. This spanning of scales allowed us to investigate the scale dependence in flux intrinsic variability and propagated uncertainty as well as the interactions taking place between these properties, including the specific role played by ecosystem type. To our knowledge, this study is the first of its kind in the vast region of West Africa, and probably also in the entire continent.

As uncertainty analysis is a developing field, techniques of analysis are still required to address the variety of situations encountered and questions raised. In its broadest sense, uncertainty in a variable taken at a given time/space location can be defined as the distribution of departures of all possible values of the variable from their mean—the variable’s expectation—at that location. Hence, uncertainty adds yet another dimension to an already largely multidimensional modeling problem (time, 2D space, scales, and variables), and must be analyzed in light of the entire distribution of the mother variable itself. Such analysis calls for innovative techniques to explore and synthesize information from large multidimensional datasets, and a case study like this one is an opportunity to contribute to that development. Therefore, as an extension of the initial objective of assessing uncertainty arising in modeled surface fluxes from state-of-the-art field-based rainfall simulation in the Sahel, a complementary objective was to investigate new means of capturing and expressing uncertainty meaningfully in a high-dimensional space, through innovative graphical representation and synthetic metrics.

Section 2 describes the study case, the modeling protocol, and the methodology of output analysis—including the metrics and naming conventions used. Results are presented and interpreted in Section 3. Section 4 discusses and concludes on these results.

2. Materials and Methods

2.1. Study Area and Period

The study area is a 2530 km² rectangle (0.54 × 0.39 degree in longitude/latitude) located east of Niamey in the Fakara district of the Republic of Niger (Figure 1a). This area is part of the Niger mesosite (ACN) of the AMMA-CATCH long-term regional eco-hydro-meteorological observatory [4,6,47] and includes the Wankama supersite (~2 km²) at which detailed field monitoring of surface and subsurface processes has been carried out since 2005 (e.g., [22,48,49,50,51,52,53]. The area represents the Central Sahelian conditions in the ALMIP-Phase 2 multi-site regional land surface model intercomparison experiment [54]. It has been the topic of several other modeling exercises at different scales (e.g., [11,20,22,26,27,55,56,57,58,59]), and the support for development or evaluation of various remote-sensing based surface flux estimation methods [17,60,61,62,63,64].

Cappelaere et al. [6] give a detailed description of the study area. Large sandy valleys with gentle slopes are interspersed with small flat lateritic plateaus. Runoff is Hortonian and limited to small endorheic catchments, of at most a few km² area, where it collects in ponds. Land cover is an intricate mosaic of millet crop, fallow bush (temporarily uncultivated land), and bare soil in the valleys, while tiger bush relics and large bare soil patches occupy the plateaus. Climate is tropical semiarid, with a short rainy season from June to September which is followed by a long, almost totally dry season [22].

This study covers the 6-month period from 15 June to 14 December 2005 (183 days, i.e., exactly a half-year duration), consisting of the 2005 Sahelian rainy season plus the first two months of the subsequent dry season, and running roughly from the so-called African monsoon jump at the Northern Hemisphere summer solstice [65] to the following winter solstice. The period was part of the Enhanced Observation Period of the AMMA international field campaign [66], giving the study the benefit of an extensive set of field observations, for both model forcing and calibration.

2.2. Generating Stochastic Rainfields

The ACN observatory includes a permanent network of high-resolution recording raingauges [6], recognized as being unique in West Africa for its density and duration. Although gauge density has varied over the years, a core set of 32 tipping-bucket gauges, evenly distributed over a total area of 160 × 120 km² that includes the study domain (Figure 2 in Cappelaere et al. [6]), has been operating continuously since 1990. Times series of bucket tips were processed into continuous rainfall series at 5-min resolution [67].

Based on this time- and space-rich record, methods were developed to generate space–time rainfields through raingauge interpolation and stochastic simulation over the ACN observatory [46,68,69]. Despite the network density, estimating rainfall distribution in space and time inherently comes with substantial uncertainty, which may vary strongly with input and output resolution in both dimensions, due to the strong variability of convective Sahelian rainfall. To reflect this uncertainty faithfully, ensemble scenarios of equally possible rainfields, given all available observations, were generated using a two-step process. First, rainfields were simulated at the event timescale over a ~1 km (more precisely 0.1 degrees) resolution grid covering the entire domain. This was achieved using meta-gaussian geostatistical simulation [68], with point conditioning by gauge event rainfall [46]. All 32 stations—including the 22 located outside the study domain—were used to condition the simulated rainfields. Fifty stochastic replications per event were simulated. These ensemble event rainfields were then time-disaggregated at each grid node by scaling to the stochastic event rainfall an optimum event hyetograph shape derived by deterministic spatial interpolation from observed 5-min rainfall intensities. The interpolation consists of either Eulerian or Lagrangian kriging, depending on their respective performance as assessed for each event through a cross-validation scheme [69]. The conditional simulation scheme from Vischel et al. [46] has been used by others [70] to evaluate the effects of rainfall uncertainty on the performance of calibrated rainfall–runoff models.

Figure 1c displays the spatial distribution of time-aggregated ensemble rainfall over the study period, as maps of the ensemble mean and standard deviation. A sample time-aggregated rainfield from the simulation ensemble is provided as Supplementary Material (Figure S1).

2.3. Producing LSM Simulations

The above rainfields were used to generate ensemble mesoscale land surface simulations with the SEtHyS-Savannah LSM. This model was developed to simulate the energy and water processes in conditions such as in the Sahel [26]. It was previously applied with deterministic rainfall to the same study domain and period [27]. It is a one-dimensional, column model that makes a suitable compromise between the level of process description and the number of input parameters for mesoscale implementation. It uses a three-source modeling approach, as recommended by Verhoef and Allen [23] for Sahelian ecosystems. The canopy is described as two superimposed layers, one for crops, bushes and/or trees, and one for the grass understory. The model solves for three distinct energy budgets, one for each canopy layer and one for the underlying soil. Fractions of incoming downward shortwave and longwave radiation absorbed within each canopy layer are computed using a shielding factor [71]. For turbulent fluxes, an aerodynamic resistance scheme is applied to the three-source modeling structure [26], using the Ball [72] and Collatz et al. [73,74] formulations of stomatal resistances for C3 and C4 plant types, respectively. The soil column is divided into two horizontal layers with prescribed depths. The top layer contains roots from both vegetation layers, while the lower layer has roots from the top vegetation layer only. During drying, the development of a mulch sub-layer is simulated within the topsoil layer, with a resistance to water vapor diffusivity proportional to mulch thickness. Two separate mulch sub-layers may be simulated in case of alternating wetting/drying sequences. Infiltration-excess runoff, taking into account a surface crust conductivity, is output directly, without subsequent routing at any larger scale. The model was parameterized and calibrated at the plot scale, based on field observations from the ACN Wankama catchment for major ecosystem types [22,53], and on the existing literature for the South-West Niger area (e.g., [6]). Compared to Saux-Picart et al.’s implementation [26,27], model parameters were further tuned thanks to new data including soil moisture profiles [75].

For the spatial distribution of model parameters and inputs other than rainfall, we kept the approach used by Saux-Picart et al. [27]. A land-cover map was derived through unsupervised classification from a 20 m resolution SPOT-HRV image of 12 October 2005 (end of the wet season, the best period for this purpose) and area aggregation (tiling) at the same ~1 km (0.1°) resolution (Figure 1a), yielding land cover proportions in each km² cell. Eight different ecosystem types were distinguished in the simulation. For simplicity of the following analysis, three of these classes—each below 5% surface area—were aggregated as ‘other’, and the next least represented one—degraded plateau slopes, 5.2% surface area—was integrated into the bare soil class (Figure 1a). The three other major types are crops (millet essentially), young fallow, and old fallow (these two were discriminated on a vegetation density basis from the SPOT image, as density increases with age [27]). Corresponding time-dependent maps of the leaf area index (LAI) were inverted successively from seven clear-sky SPOT-HRV images distributed along the study period, using a neural network trained against a radiative transfer model [76], and were linearly interpolated in time (Figure 1b). For lack of reliable space-distributed data, meteorological forcing (shortwave and longwave radiation, wind speed, air humidity, temperature, and pressure) was taken from the observations at a 30 min timestep in the Wankama catchment [22], and applied uniformly over the mesoscale domain. The model was run at that time resolution for each ecosystem tile in each kilometer cell.

2.4. Analyzing Simulation Outputs

Please note that mathematical formulations of the computational steps described hereafter are detailed in Appendix B.

2.4.1. Fluxes

Flux types analyzed in this study were (i) for the water cycle: rainfall P, runoff Q, direct soil evaporation Es, and plant transpiration Tp (henceforth abbreviated to evaporation and transpiration, and summing to evapotranspiration), as well as drainage Dr; and (ii) for the energy cycle: net radiation Rn, sensible heat flux H, and latent heat flux LE. Ground heat flux was considered as less significant at timescales analyzed here, hence it is not presented. Direct canopy evaporation is quite negligible in this environment [22]. Evapotranspiration and latent heat flux are equivalent modulo the unit latent heat of vaporization. Whenever possible in the figures, water and energy fluxes are displayed together on the same graphs using this conversion factor.

2.4.2. Ecosystem Types

To simplify the presentation, only the four main ecosystem classes (LC)—bare soil, crops, young fallow, and old fallow (Section 2.3; Figure 1a)—are shown when results are displayed for separate ecosystem types. These four classes cover >91% of the surface area. All eight classes are accounted for, through area-weighted tile aggregation in each cell (Appendix B), in the mixed-ecosystem outputs which are the main results presented here. Unless otherwise specified explicitly as referring to “separate ecosystems”, fluxes are taken for “mixed ecosystems”.

2.4.3. Aggregation Scales

A large range of scales is investigated in this study, both in time and space. Over space, both the ~1 km² cell scale and that of the entire 2530 km² simulated domain were considered, denoted as the local and the meso scales, respectively. Timescales were analyzed from the daily to the half-year scale, including intermediate scales, in particular the 10/11 day scale. Although diurnal cycles were explicitly resolved in the simulations, sub-daily variability was not considered in the following analysis; hence “full resolution” (or “fully-resolved” fluxes) refers to the daily/cell scale. Combinations of these time and space scales were covered to build a comprehensive picture of scale effects on fluxes and their uncertainties. Note that for all flux types, including runoff, space and/or time upscaling from the computational resolution were obtained by simple arithmetic averaging of raw computed fluxes (Appendix B).

Table 1. Glossary of flux aggregation levels ¹ over the various dimensions (xy, t, e, and LC) of the simulation set, as considered in the analysis of results.

Dimension		Aggregation Level
	«Non-Aggregated»	Partially Aggregated	«Aggregated»
space (xy)	cell (~1-km2)	-	meso (2530 km2)
time (t)	Day	10/11 day	half-year (183 days)
ensemble (e)	Member	-	mean (=expectation, expected flux)
ecosystem type (LC)	separate ecosystem	-	mixed ecosystems (cell area-weighted average)

¹ Please see Appendix B for mathematical expressions of these aggregation operations.

summarizes the different aggregation levels considered in the analysis along each dimension of the simulation set.

2.4.4. Uncertainty Measures

This section defines the mathematical criteria (“measures”) together with the terminology used in the analysis to characterize uncertainty. Mathematical developments are provided in Appendix B. In the following, the generic term “uncertainty” is synonymous with ensemble member variability, and flux “expectations” designate the means over the ensemble member dimension. Describing uncertainty propagation properties in a high-dimensionality space (including, e.g., multiple flux types and time and space scales for both uncertainty input and outputs) requires compact measures that efficiently characterize uncertainties across these scales. Several synthetic uncertainty measures, both dimensional and dimensionless, are defined to summarize ensemble member anomalies (departures from flux expectations), as well as corresponding propagation properties from rainfall to any other flux type.

The measures may be taken either at each point in the time/space domain at any given scale (“elemental” measures), or lumped over one (“partial lumping”, over either time or space) or both dimensions (“full lumping”) of the domain for that scale. Hence, they are distributed over the non-lumped dimension(s), if any. Note that lumping a dimension is of no effect when the scale is the largest in that dimension. While both flux scale aggregation and uncertainty lumping represent a (different) type of integration over time and/or space, only the former amounts to a change in scale as referred to in this paper, the second corresponding rather to a synthesis of multiple uncertainty measures. Unless specified as lumped, a measure is elemental.

The dimensional measure used is the uncertainty’s first moment, i.e., the standard deviation over the ensemble dimension (standard deviation of ensemble members), which is referred to as “standard uncertainty” (USD). Lumping of USD is performed by simple quadratic (root mean square) averaging, amounting to ensemble variance aggregation (Appendix B). Two dimensionless uncertainty measures are considered, compared and evaluated to determine which most efficiently describes the uncertainty generation and propagation properties. These two measures are obtained by normalizing the standard uncertainty USD for a given time/space scale and a given level of lumping (none, or partial over one dimension, or full over both dimensions), with either the corresponding quadratically-averaged flux expectation or the corresponding total standard deviation (generated by all variability sources including uncertainty) (Appendix B). These two dimensionless measures are referred to as “uncertainty coefficient” (UC) and “uncertainty fraction” (UF), respectively. UC is commensurate with a coefficient of variation and numerically comes down to a quadratic average of elemental coefficients of variation (elemental UCs) weighted with corresponding expected fluxes. UF corresponds to a square-root variance fraction, which, when squared, sums to unity with the index of ensemble similarity (Ω) introduced by Koster et al. [77] and analyzed by Yamada et al. [78].

Finally, for any flux type at any scale and any time-lumped uncertainty measure, a propagation factor—from rainfall to that flux type—is defined as the ratio of that measure to the same measure for rainfall (see formulae in Appendix B).

2.4.5. Analysis and Presentation of Results

Analysis of simulation outputs is presented in Section 3 as three successive steps. In the first step (Section 3.1), the distributions and variability of plain flux intensities were examined directly at the various scales, including a global variance decomposition of the simulated ensemble set (see, e.g., [79]). The second step (Section 3.2), specifically analyzed the distributions of flux uncertainties, using the elemental or partially-lumped uncertainty measures defined in Section 2.4.4. Relationships to flux magnitude, rainfall uncertainty, and ecosystem type were investigated. In a third and final step (Section 3.3), uncertainty measures, including propagation factors, were taken fully-lumped to synthesize their scale behavior and scaling properties.

To support this presentation, original graphical representations were devised to display as much as possible of the multidimensional variations of fluxes and of their uncertainties, as well as relationships between them.

3. Results

3.1. Distributions and Variability of Flux Intensities

3.1.1. Distributions of Flux Ensemble Members and Means

The study period included 80 rainy days (i.e., with rainfall recorded at one raingauge at least), distributed from 15 June to 13 October 2005. Days with runoff (at least in a part of the domain and members) were a subset of 53 of these rainy days. While soil evaporation started with the first rainfall and essentially vanished after mid-October, substantial transpiration and latent heat flux persisted until early November. Drainage first appeared on 24 June, peaking mid-August, and remained significant until late October, then slowly decayed.

A summary description of marginal (i.e., unconditional) distributions of these variables over the entire space–time–ensemble simulation domain is provided in

Table 2. Summary statistics ¹ (top to bottom: range, mean, coefficient of variation (CV), skewness, kurtosis ²) for marginal distributions of analyzed variables (inner rows: mbr. = ensemble members, exp. = ensemble means, USD = elemental standard uncertainties) for the different flux types (left to right: P to LE) and scales (inner columns: c,d = cell/day, m,d = meso/day, c,p = cell/half-year).

Variable		Rainfall, P			Runoff, Q			Evaporation, Es			Transpiration, Tp			Drainage, Dr			Net Radiation, Rn			Sensible Heat, H			Latent Heat, LE
Scale		c,d	m,d	c,p	c,d	m,d	c,p	c,d	m,d	c,p	c,d	m,d	c,p	c,d	m,d	c,p	c,d	m,d	c,p	c,d	m,d	c,p	c,d	m,d	c,p
	mbr.	0	0	1.0	0	0	0	0	0	0.8	0	0	0	0	0	0	−0.8	−0.4	2.1	−6.9	−2.9	0.7	0	0.1	0.9
		121	41.6	3.5	82.5	16.4	1.4	10.8	7.2	1.7	3.2	1.4	1.0	6.6	0.2	0.4	6.0	5.3	3.2	4.3	2.8	1.5	12.5	8.2	2.6
Range (mm/d)	exp.	0	0	1.6	0	0	0.1	0	0.1	0.9	0	0	0	0	0	0	−0.5	−0.4	2.2	−6.1	−2.8	0.9	0	0.1	1.2
		82.9	33.4	2.6	32.0	10.6	1.1	9.8	6.9	1.6	2.7	1.2	0.9	0.7	0.1	0.1	5.8	5.2	3.0	3.6	2.7	1.4	11.4	8.0	2.1
	USD	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
		19.6	3.5	0.3	15.3	2.3	0.2	3.3	0.5	0.1	1.0	0.1	0.1	1.0	0.0	0.1	1.0	0.1	0.1	3.0	0.3	0.1	3.7	0.5	0.2
Mean (mm/d)	mbr.	2.19			0.38			1.26			0.44			0.02			2.74 (78.4 W/m2)			1.12 (31.9 W/m2)			1.70 (48.5 W/m2)
	exp.
	USD	1.22	0.33	0.22	0.42	0.12	0.11	0.29	0.07	0.07	0.06	0.01	0.03	0.017	0.004	0.016	0.11	0.02	0.04	0.19	0.04	0.05	0.34	0.08	0.09
	mbr.	2.95	2.39	0.16	5.73	3.41	0.50	1.15	1.04	0.10	1.03	0.83	0.39	3.71	1.38	1.77	0.53	0.52	0.07	0.77	0.63	0.12	1.03	0.94	0.12
CV	exp.	2.61	2.37	0.12	4.36	3.29	0.42	1.09	1.03	0.09	1.00	0.83	0.38	2.35	1.29	1.13	0.52	0.52	0.06	0.71	0.63	0.11	0.98	0.93	0.11
	USD	2.24	2.05	0.21	3.28	2.67	0.28	1.30	1.20	0.17	1.38	1.03	0.41	2.50	1.55	0.97	1.05	0.88	0.24	1.40	1.27	0.22	1.23	1.13	0.22
	mbr.	4.34	3.36	−0.02	9.02	5.41	0.83	1.47	1.29	0.02	1.04	0.29	−0.22	13.2	2.01	3.64	−0.05	−0.12	−0.60	−1.40	−1.31	0.17	1.16	0.93	0.04
Skew.	exp.	3.64	3.30	−0.53	6.72	5.06	0.77	1.38	1.28	−0.01	0.96	0.28	−0.29	4.74	1.46	1.90	−0.08	−0.12	−0.69	−1.41	−1.31	0.20	1.08	0.92	−0.06
	USD	2.47	2.14	−0.35	4.38	3.67	0.44	2.25	2.16	−0.61	2.44	1.15	0.22	5.53	2.21	1.34	1.61	1.02	0.20	2.94	2.60	0.03	2.11	2.03	0.04
	mbr.	26.6	15.9	2.63	109.	39.1	3.74	5.12	4.65	2.82	3.41	1.74	2.22	376.	7.21	22.0	1.80	1.75	2.49	8.72	7.68	2.47	4.22	3.62	2.58
Kurt.	exp.	18.6	15.3	2.06	57.9	33.1	3.41	4.95	4.64	2.64	3.15	1.73	2.18	34.5	3.68	7.93	1.80	1.75	2.46	9.18	7.74	2.35	4.10	3.61	2.21
	USD	8.73	6.93	3.04	24.6	18.6	3.09	9.51	9.11	4.00	11.4	3.67	2.44	48.9	7.60	4.48	5.98	3.62	2.98	15.3	11.8	2.84	8.89	8.53	3.09

¹ Unit for ranges and means is mm·day⁻¹ (energy fluxes expressed in equivalent evaporated water, also in W·m⁻² for means); other statistics are dimensionless. ² Kurtosis expresses tailedness and varies from 1 for a finite, uniform distribution, to 3 for a normal distribution, and above that value for a more narrow-peaking, fatter-tailed distribution.

through their ranges and first four statistical moments (see rows “mbr.” and “exp.” for ensemble members and means, respectively). Flux distribution at various time/space scales (from day to half-year and from cell to meso) is depicted through a set of three figures (Figure 2, Figure 3 and Figure 4), showing selected frequency quantiles for both ensemble members and ensemble means, conditioned to time in the season for the first figure and globally over the whole time/space domain for the last two.

Figure 2 maps over time the fully-resolved (cell/day) flux population—successively for each flux type—as the temporal courses of quantiles taken within a running 11-day window and the whole spatial domain, and over the ensemble members (dashed lines) or the ensemble means (plain lines). Hence three different variability modes are highlighted: seasonality (see time variability), combined spatial/short-timescale variability (<11 days; see quantiles for ensemble means), as well as propagated uncertainty (see quantile differences between ensemble members and means). The running mean, representing the meso/11 day-scale expected flux, is also represented (black dotted line).

For the same fully-resolved fluxes now taken globally over the whole space/time domain, Figure 3 displays the distributions of ensemble members conditioned to the value of the ensemble mean (inner 2-D crossplots) as well as corresponding marginal distributions of ensemble members or means (top and right 1-D whiskers, respectively). Hence the figure decomposes the overall distribution of ensemble members into the direct (first-order) effect of spatiotemporal variability (distribution of ensemble means) and the rainfall uncertainty effect (conditional distributions), highlighting the relationship of uncertainty to flux magnitude.

Similarly to Figure 3 but for the half-year timescale, Figure 4 shows the frequency distributions of time-integrated, cell water (Figure 4a) and energy (Figure 4b) fluxes for ensemble members and means (color plots in inner boxes, including both marginal distributions of members and means as well as conditional distributions of members to means, as in Figure 3). Figure 4 also shows distributions of these members and means conditioned to ecosystem type (1-D whiskers in outer plot areas). Finally, small vertical grayscale bars at the center of each 2-D crossplot (inner boxes in Figure 4) depict fully space/time-integrated member and expected fluxes (vertical and horizontal axes, respectively).

From these three figures, it can be seen that the effect of uncertainty (ensemble variability) is generally quite significant, with respect both to flux magnitudes and to the other variability modes, especially for the coarser resolutions. Considering changes in variability and skewness resulting from integration over one or several of the three dimensions –time, space, and uncertainty—(Table 2, Figure 2, Figure 3 and Figure 4) suggests that the time dimension predominates in the variability and shape of the fully-resolved ensemble fluxes for most flux types. Relative weights of these three variability sources and of their interactions are assessed quantitatively in Section 3.1.2, while space/time-scale effects on uncertainty (ensemble variability) are specifically investigated in Section 3.2 and Section 3.3. Interpretation of Figure 2, Figure 3 and Figure 4 in terms of patterns of first-order functioning of water and energy cycles in the study domain, is developed in Appendix C.

3.1.2. Global Variance Analysis of the Simulated Flux Ensemble: Confronting Variability Sources

In this section, flux variability in the entire, full-resolution, mixed-ecosystem simulated set is analyzed by decomposing total flux variance into its three orthogonal dimensions, namely time (t), space (xy), and ensemble (e). Flux variability in these different dimensions originates from corresponding variability in model forcings, namely: (i) time dependencies in rainfall, meteorology, and vegetation phenology, (ii) spatial variability in rainfall and ecosystem distribution, and (iii) rainfall uncertainty, for t, xy, and e respectively. This variance decomposition allows us to compare the respective effects of these different sources and to highlight those due to their interactions, through lumped measures over the whole time/space simulation domain. Analyzing temporal and spatial effects informs on how fluxes scale over time and space. It also highlights the time/space scaling of the third factor’s effects, namely flux uncertainty (taken here as the fully-lumped squared standard uncertainty, USD²).

Figure 5 displays this global variance decomposition for the different water and energy fluxes. Interpretation of the various component terms is detailed in the figure legend (right side). Taking, for instance, the uncertainty factor “e” for each flux type, the 1st-order uncertainty term amounts to the time/space-independent variance, i.e., uncertainty on the fully-integrated, meso/half-year flux. Adding either the «t-e» or the «xy-e» 2nd-order interaction term leads to the aggregated uncertainty on time or space-distributed fluxes, i.e., to uncertainty at meso/day or cell/half-year scales, respectively. Including both interaction terms as well as the 3rd-order «t-xy-e» (entire «e» bars in Figure 5) provides the finest, cell/day scale aggregated uncertainty. Similarly, the 1st-order term for the time or space factor represents the variance of the time-distributed (meso/day) or the space-distributed (cell/half-year) expected fluxes, respectively (decomposition of the latter into precipitation and land-use contributions is also reported in the figure). Adding these two 1st-order terms and the «t-xy» interaction produces the variance of expected fully-distributed (cell/day) fluxes. Entire time (left) or space (middle) bars represent total temporal or spatial variance in fully-resolved ensemble fluxes, aggregated over all ensemble members, and cells or days, respectively. Other combinations of terms provide further meaningful expressions of partial variances in the three-factor decomposition space (see legend).

As suggested in Section 3.1.1, the time factor (time-wise variability in forcings) dominates total variance for all variables—albeit less sharply for Dr (Figure 5). Its 1st-order component alone (time factor without interactions with space and uncertainty, i.e., ensemble- and space-averaged fluxes) even exceeds the other factors’ total variance except for Q and Dr. The relative effect of uncertainty on time variability is comparatively lower than that on space variability (see next paragraph). The space factor comes next in importance in the variance budget, but with a lesser 1st-order contribution (cell/half-year expectations), and mostly 2nd-order interaction with time (cell/day expectations) and 3rd-order interaction (full resolution ensemble). Direct interaction with uncertainty is less substantial, yielding moderate added spatial variability from expectations to ensemble at the cell/half-year scale (max. factor 2.4, for Dr). The uncertainty-induced spatial variability increase is always higher at daily scale (factor up to 8.8, for P), due to the 3rd-order term. Decomposition of spatial variability alone (1st-order term, i.e., for cell/half-year expectations) between precipitation and land-use largely varies with flux type (Figure 5). While the land-use factor dominates strongly for Tp and Rn, and to a lesser degree for Q, the two factors are rather balanced for Dr and H, and precipitation dominates only for Es and LE, not to speak of P itself. Note that this effect of land use (and hence the whole spatial variability) is largely tempered by the mixed-ecosystem composition of cells.

For the uncertainty factor, 1st-order terms (meso/half-year uncertainty) are much smaller than the interaction terms. Second-order interaction—hence downscaling—with time is stronger than with space for P, Q, and to a lesser extent Es and H, weaker for Dr, and similar for Tp and Rn. Further separating the spatial factor into its precipitation and land-use induced components, confirms (not in figure) that this interaction of space with uncertainty arises from the first component (precipitation) alone essentially, albeit also from the 3rd-order interaction with the two spatial components for Dr. Most important, however, is the 3rd-order interaction with both time and space, which dominates by far the composition of ensemble variability. These successive jumps in the magnitude of component uncertainty terms when moving from first to second and to third order induce very strong scale dependence of total uncertainty for all flux types.

Variance decomposition highlights the relationships existing between the variance components induced by the different factors. For instance, the decrease in variable uncertainty with spatial upscaling from cell/half-year to meso/half-year scales is identical to the decrease in spatial variability between ensemble members and their expectations at cell/half-year scale (see Figure 5 legend, common 2nd-order term to the uncertainty and space factors, in green). The same holds between the time and uncertainty factors (light-blue 2nd-order term), or between time and space (dark blue 2nd-order term), or still between the three for the 3rd-order term (gray). Also, considering the different types of variability simultaneously allows us to compare their respective magnitudes. For instance, it can be seen that uncertainty comes relatively close to spatial variance in the ensemble set (ratio ≥ 30%) for P, Q, Dr, and Es at cell/day or cell/half-year scale, also for daily LE, H, and Rn. It is smaller for Tp at both timescales and for the energy fluxes at a half-year scale. This ratio generally decreases with time upscaling and falls below 4% for Tp and Rn at cell/half-year scale. If comparing now with spatial variability of expected fluxes: (i) at cell/day scale, uncertainty comes out larger for P, Q, Dr, Es, and similar for LE and H; (ii) at cell/half-year scale, it is still larger for Dr and similar (ratio of 47%–72%) for P, Q, and Es.

3.2. Patterns of Distributed Uncertainties

3.2.1. Distributions of Elemental Uncertainty

A first picture of the distribution of ensemble flux anomalies arising from rainfall uncertainty—the departures of ensemble members from expected fluxes—is provided by the deviations of frequency quantiles from the 1:1 line in Figure 3 and Figure 4, in relation to ensemble means, for the various flux types at cell/day and cell/half-year scales, respectively. The time-dependence of uncertainty-generated anomalies at cell/daily-scale is depicted by the deviations of frequency quantiles from means in Figure 2. Ensemble anomalies can be condensed to the elemental (unlumped) standard uncertainties USD, to characterize the distribution of uncertainties over the time/space simulation domain for the various time and space scales. Figure 6 displays the frequency distributions at full, cell/day resolution of this elemental USD, both conditioned to ensemble means and as marginal distributions over the whole domain. Figure 7 shows, for the cell/half-year scale, the direct scatter of elemental USD against means, with color intensity mapping by the rainfall uncertainty coefficient UC (Section 2.4), as well as their marginal distributions both for the ecosystem mosaic and for individual ecosystems. Table 2 gives summary statistics of elemental USD at various scales.

3.2.2. Relating Flux Uncertainty to Flux Magnitude

The figures displaying ensemble member variability (Figure 3 and Figure 4) or standard uncertainty (Figure 6 and Figure 7) versus expected fluxes, suggest that space/time distributions of uncertainty are strongly imprinted by those of the underlying flux magnitudes. This is also hinted at by the general relative similarity in dimensionless statistics (especially coefficient of variation and skewness, see Table 2) between the fluxes—ensemble members and means—and their uncertainties, for any given flux type and scale.

At both timescales—daily or half-year—uncertainties generally grow with flux magnitudes, more or less clearly depending on flux type. On a daily scale however (Figure 6), for most flux types, they tend to level off for large magnitudes or even to decrease for the largest ones. At least part of this originates from this effect being present in P, meaning that rainfields are likely to be better constrained by the raingauge network for the strong events than for those with lower daily intensities. This results in a saturation-type effect on extreme member values for most flux types, which limits the differences in the range between the marginal distributions of means and members (Figure 2).

A particularly strong relationship—from flux magnitude to standard uncertainty—is obtained for P, Q, and also Dr when only the time dimension is expounded (i.e., meso/day or space-lumped cell/day scales—not shown). Conversely, when the space dimension is made explicit (cell resolution, no spatial lumping), the geometrical relationship to the gauge network layout also comes into play, as it controls the degree of source uncertainty and hence blurs the relationship to flux magnitude.

Except for Rn, uncertainties in daily fluxes (Figure 3 and Figure 6) generally tend to be comparatively lower for the most frequent flux magnitudes (see, e.g., the case of H).

When intercomparing flux types, Figure 7 highlights at cell/half-year scale the stronger response to rainfall uncertainty in relative uncertainty (uncertainty coefficient UC, represented by the slope from the origin in Figure 7) for Q and especially Dr compared with Es and Tp. Uncertainties are amplified relative to P for Q and Dr but attenuated for Es and Tp. Similarly, for energy budget components, LE and H display much higher relative sensitivity to rainfall uncertainty than Rn. It is shown in Section 3.3 that, taken relative to flux magnitude (i.e., as UC) as well as to its space/time variability (i.e., as UF), uncertainty is comparatively highest for Dr and lowest for Rn at all scales. It is also shown that UC makes the most powerful dimensionless uncertainty measure, and is therefore used intensively in the following.

3.2.3. Relationship to Rainfall Uncertainty: Uncertainty Propagation Patterns

The color darkness patterns in Figure 7 provide evidence of a relationship in flux USD, and even more so in flux UC (slope from the origin), to rainfall UC, for most flux types at the cell/half-year scale. The flux uncertainty-to-expectation relationship at cell/daily scale is also well structured by the cell/half-year rainfall UC for Q and Dr but less so for the other flux types (not shown). This section further analyzes the relationship in UCs between rainfall and each dependent flux type, when projected either in the temporal or the spatial dimension. In time-wise analysis, time-distributed fluxes are considered either integrated over space (mesoscale) or at cell scale with the uncertainty measure being lumped over space. Conversely, the space-wise analysis considers both daily cell fluxes with time-lumped uncertainty and time-integrated cell fluxes.

• Time-wise analysis (temporal patterns).

Figure 8 displays the seasonal courses of space-independent UCs for 10-day fluxes taken either at cell (space-lumped UCs, Figure 8a) or meso (space-integrated fluxes, Figure 8b) spatial resolution. Despite the different space scales, note the similarity in shapes between the two graphs. The source, rainfall uncertainty follows a time pattern that is roughly opposite to that of the flux itself (Figure 2a), which is a consequence of P uncertainty scaling tightly with magnitude (Section 3.2.2) according to a concave (roughly square-root) relationship. Q closely follows the P pattern, with a relatively constant amplification factor of around 2 and 3 for cell and mesoscale fluxes (Figure 8a,b), respectively. The P pattern can be seen as propagating with more or less distortion/smoothing and lag to most of the other flux types. Only Dr appears to be less connected to the P pattern, with a much smoother course presumably due to a more time-integrative response of this flux type to rainfall. Logically, Es and LE show a high correlation, as do, less expectedly, Tp, and Rn (with uncertainty in Rn essentially coming from the upwelling thermal radiation). The abrupt hollow in the curves for several flux types after Julian day 280 is due to some light late-season rainfall with no injected uncertainty (as can be seen in Figure 2a). Note that despite the different behaviors, flux type ranking remains overall relatively time-invariant. Amplitudes of ranges in UC largely vary with flux type, e.g., for cell-scale fluxes (Figure 8a), from a factor of ~3 for P and Q to as much as ~7 for Tp.

• Space-wise analysis (spatial patterns).

Using again the UC dimensionless measure, Figure 9 displays the spatially-variable (one point per cell) relationship to rainfall uncertainty of all water and energy flux uncertainties, for daily-scale fluxes (time-lumped measure, Figure 9a) or half-year fluxes (Figure 9b). On both graphs, diagonal projection to the right or top graduated axes provides the ratio of the two confronted UC values, i.e., the corresponding UC propagation factor. For time-integrated fluxes (Figure 9b), quantile distributions of the propagation factor are also drawn on the right side of the graph (with diagonal reading), for the mixed-ecosystem cells as well as for individualized ecosystems.

It can be seen that for both timescales the relationships of flux-to-rainfall uncertainties are fairly to very strong, depending on flux type. Dr is the least strong, including a few points scattered away from the main cloud: this is likely due to this flux acting as a residual from all other water cycling processes, hence integrating all variability sources that affect the various other fluxes. Because of spatial independence in simulated flux processes (no interaction between cells, unlike over time; spatial relationships stem from forcings only), flux uncertainty in a given cell depends on the rainfall ensemble for that cell only. Figure 9 shows that this dependence is rather well summarized using the UC uncertainty measure applied to rainfall and to the various dependent fluxes.

With daily scale fluxes, even for flux types for which the relationship is very strong (e.g., Es), the corresponding propagation ratio generally varies somewhat with the level of P uncertainty, e.g., decreasing from ~0.9 to ~0.5 (i.e., increased attenuation) for Es when rainfall UC sweeps the recorded range (0.15–0.7). Only for Q is the trend slighter, and opposite (propagation ratio increasing from slightly below to slightly above 2). Relationships generally appear even stronger and propagation ratios more invariant (linear relationship) when fluxes are considered at the time-integrated scale.

Given the generally moderate variations in propagation ratios for any given individual flux type in Figure 9 (especially 9b), compared to differences between flux types, a clear hierarchy in the uncertainty measure emerges over space between flux types. This hierarchy, which varies somewhat with the timescale, is hence well summarized by the fully-lumped UC, which is mapped as big circles in Figure 9 for both timescales. These are well centered within the corresponding space-resolving clouds—except for Dr for reasons explained above. Hence the analysis of fully-lumped measures as developed in Section 3.3 can be considered to largely represent the properties of space-distributed variables as well, especially in terms of propagation effects.

Note that for both timescales, UC propagation for spatially-distributed fluxes come out nearly everywhere as amplification for Q and Dr but attenuation for all other fluxes, with only very rare exceptions for Dr at daily scale and Tp at half-year scale. This quasi-invariance also holds over time for the 10-day scale (Figure 8).

The maps in Figure 10 illustrate how the spatial patterns in time-lumped UC of daily-scale P propagate to the other flux types. The marked pothole structure of P uncertainty is largely preserved for Es and LE as well as for H, Rn, and Q and to a lesser extent for Tp, but is considerably blurred in Dr.

3.2.4. Effects of Ecosystem Type on Uncertainty

The discussion of distributed uncertainties has so far considered a cell as an integrated mix of multiple ecosystems. This subsection discusses the information characterizing the different ecosystem types for time-integrated fluxes in Figure 7 and Figure 9b.

In terms of USD (Figure 7), significantly higher values among ecosystem types are obtained for the following combinations: bare soil for Q, crops for Dr, old fallow for LE and H. All of these four combinations but the last one—old fallow with H—also correspond to the highest flux magnitudes of all ecosystem types, for the given flux type (Figure 7, top). Conversely, statistically lower USD among ecosystem types is obtained from bare soil for all energy components (Rn, H, and LE), all of which also correspond to lower flux magnitudes. This property of bare soil could also be extended to Tp and Dr, but in this case, flux magnitudes are irrelevant or insignificant for those flux types. Differences in evaporation USD between cover types are not significant, despite larger differences in flux magnitude.

In terms of the relative uncertainty UC, it can be seen in Figure 9b that distributions of UC propagation factors are lowest from bare soil for Q, Es, and Rn, from crop fields for Dr and Tp (values for these latter flux types are not quite meaningful for bare soil, see above), and from both cover types for H. Es from bare soil even displays, in statistical terms, both the lowest USDs (slightly) and the highest flux magnitudes of all ecosystem types (Figure 7). Conversely, relative sensitivity is statistically highest from the old fallow type for Dr (but with very low flux magnitude) and H, from both young and old fallow for Tp and Rn, and the crop type for Q (lowest USD but also lowest magnitude). Propagation factor distributions are hardly distinguishable between these three vegetated ecosystem types for Es, and to a lesser extent for Rn. For Q, differences between all types are much smaller in the middle of the wet spells than in the drier spells (not shown).

3.3. Scaling of Fully-Lumped Uncertainty Measures and Propagation Factors

The global variance analysis of Section 3.1.2 (more precisely, the combinations of variance component terms for the “e” factor in Figure 5) as well as the uncertainty statistics in Table 2, indicate the effects of scale on the uncertainty measures. They suggest for instance a predominance of timescale over space scale effects for some flux types. This section analyzes quantitatively these scale effects, for both standard and non-dimensional fully-lumped uncertainty measures as well as for corresponding propagation factors. For dimensionless uncertainty, both UC and the uncertainty fraction UF are considered and compared.

• Uncertainty measures.

Figure 11a–c (top three graphs) display the scale dependency for the three uncertainty measures USD, UC, and UF, respectively. Time and space scales are combined multiplicatively in a single scale measure—of dimension [T]. [L]²—along the horizontal axis, and include intermediate timescales. This combination of dimensions stems from considering the count of additive aggregation operations applied to full-resolution flux values when upscaling these fluxes.

While USDs logically decrease steadily with any upscaling for all flux types (Figure 11a), this is also true with UCs (although very mildly for time upscaling of Dr, Figure 11b), but not with UFs for which patterns of variation strongly vary depending on time and space scales (Figure 11c).

Rankings of flux types in USD are invariant with scales, namely in decreasing order: P (the source uncertainty), Q, Es, H, Rn, Tp, and Dr. These rankings change when considering dimensionless uncertainty UC and UF, however, as with USD they remain relatively stable across scales, and also between the two dimensionless measures (Figure 11b,c; high: Dr and Q; low: Rn, Es, and Tp; P intermediate). Only H moves from rather low in UC to intermediate in UF, and even high for some large timescales. This relative similarity in flux type rankings for the two dimensionless uncertainty measures UC and UF, and the much more stable, regular behavior of UC across scales that contrasts with the loss of significance of UF for the largest scales, justify that only UC was used as the dimensionless measure in Section 3.2.3. Also note that for all scales, the spread in UC across flux types is larger than that in UF.

Figure 11a,b suggests that, for any given flux type, the time scaling effect (slope of plain lines in the log space) on USD or UC is nearly independent of the space scale (two plain lines, for cell or meso), and vice-versa. If these scale effects are approximated as power relationships of (standard or dimensionless) uncertainty ratio r_u to scale ratio s over the scale ranges, then combined time and space scaling can be written as:

r_u ≈ (s_t) ^α_t (s_xy)^α_xy

(1)

where s_t and s_xy denote given time and space scaling ratios, r_u is the corresponding ratio for the uncertainty measure, and α_t and α_xy are the mean time and space scaling exponents for the considered uncertainty measure and flux variable.

Table 3. Time and space scaling exponents (α_t and α_xy in Equation (1)) for the different flux types and uncertainty measures (standard uncertainty and uncertainty coefficient). For each flux type and measure, are displayed: the arithmetic mean and, in brackets, ±the half-range-width between the two relevant lines in Figure 11a or Figure 11b (each exponent being computed over the full line range).

	Standard Uncertainty, USD				Uncertainty Coefficient, UC
	Space (αxy)		Time (αt)		Space (αxy)		Time (αt)
rainfall, P	−0.186	(±0.007)	−0.507	(±0.010)	−0.180	(±0.012)	−0.319	(±0.018)
runoff, Q	−0.186	(±0.002)	−0.491	(±0.002)	−0.163	(±0.010)	−0.236	(±0.015)
evaporation, Es	−0.206	(±0.006)	−0.382	(±0.009)	−0.204	(±0.008)	−0.310	(±0.012)
transpiration, Tp	−0.218	(±0.007)	−0.224	(±0.011)	−0.208	(±0.006)	−0.172	(±0.009)
drainage, Dr	−0.230	(±0.003)	−0.140	(±0.004)	−0.174	(±0.001)	−0.043	(±0.001)
net radiation, Rn	−0.207	(±0.002)	−0.231	(±0.003)	−0.207	(±0.002)	−0.208	(±0.003)
sensible heat, H	−0.207	(±0.008)	−0.367	(±0.012)	−0.204	(±0.010)	−0.332	(±0.015)
latent heat, LE	−0.205	(±0.005)	−0.339	(±0.007)	−0.203	(±0.006)	−0.277	(±0.009)

displays the values of these exponents (as the mean ± the half range-width), estimated from the line vertices (full-scale ranges) in Figure 11a,b. The generally small range-widths testify to the independence suggested above.

The space scaling effect is quite invariant across the different flux types for both USD and UC (Figure 11a,b, respectively), and rather similar between these two uncertainty measures (Table 3). For USD, the time scaling effect is stronger than space scaling for the flux types with the largest USDs (P, Q, Es, LE, and H), and vice-versa for Dr. Fluxes with intermediate USD (Rn, Tp) display comparable time and space effects. Except for Dr and to a lesser extent Tp, differences between time and space scaling are smaller for UC than for USD. Note the very small time scaling effect for drainage UC.

The departure from a value of −0.5 of exponents for USD (first two columns in Table 3 and Figure 11a) reflects the degree of autocorrelation in flux anomalies over space or time (a −0.5 exponent would correspond to no autocorrelation and a null exponent to perfect correlation with homoscedasticity; see Appendix D). Thus it can be seen that ensemble anomalies in P appear fully independent over time for scales from daily upwards, due to independent simulation of rain events (Section 2.2). This is also the case to a slightly smaller degree for Q, presumably due to a slight memory effect through antecedent moisture conditions for runoff generation (see, e.g., [80]). Tp and even more so Dr anomalies are logically the most time-autocorrelated, because of the much longer response times associated with these flux types. Hence it is the system’s processes that induce time-correlation in anomalies (through correlation in fluxes), as rainfall is devoid of such correlation over the range of timescales considered. Conversely, rainfall appears to be the main source of spatial autocorrelation for the simulated system, resulting in inverse flux type rankings (the only other possible source of spatial flux correlation being land use, which appears much less effective in this respect than does rainfall).

While it could be seen in Figure 9 that Es and H have very close UCs at cell scale, Figure 11b shows that this is true at all scales, and Figure 8 that it holds over a large part of the season when the measure is resolved over time. Transpiration UC is not very different in magnitude, with a change in ranking versus the above when timescale changes, but not when space scale changes.

• Propagation factors.

The bottom graphs in Figure 11d–f display similarly, for each of the above three uncertainty measures, the scaling with time and space of the uncertainty propagation factor for each flux type, i.e., the ratio of the corresponding uncertainty measure value to that of rainfall.

When considering USD (Figure 11d), propagation factors are as expected all below unity, the highest being for Q at a nearly scale-independent value of ~0.5: USDs in Q are about half those in P. For the other flux types, USD propagation factors slightly decrease with space scale but sharply increase with the timescale, ranging 0.01–0.16 for daily fluxes versus 0.07–0.30 for half-year fluxes (lower bounds for Dr, upper bounds for Es). This reflects the generally faster USD decrease with the timescale for rainfall than for these other fluxes (see also Table 3).

With UC, Figure 11e shows that the propagation effect (either amplification or attenuation) is consistent over the whole range of time and space scales for all flux types: Dr and Q produce amplification at all scales, while all others produce attenuation. Amplification increases with scale for Dr (from a factor ~2 to ~10 over timescales) and to a lesser extent for Q (from <2 to <3). Attenuation is stable and similar for Es and H (~0.4–0.5), but decreases (rising <1 propagation factors) with time upscaling for Tp (from ~0.3 to >0.6) and Rn (in the ~0.07–0.16 range).

When comparing propagation of the two dimensionless measures UC and UF (Figure 11e,f), the same observations made in the previous sub-section for the comparison of UC and UF logically hold for their respective propagation factors. Additionally, besides the already-noted general similarity in rankings for fine scales and its decay with increasing scales due to UF scale instability (especially over time), this behavior applies when considering the nature (attenuation versus amplification) of the propagation (Figure 11f). This property further promotes UC as a better dimensionless measure of uncertainty patterns across the range of time and space scales.

4. Discussion and Conclusions

Using (i) a stochastic time/space rainfield generator, conditioned to raingauge records and reflecting the induced rainfall uncertainty, and (ii) the SEtHyS-Savannah land surface model—both calibrated against AMMA-CATCH field-data—a fine-resolution, ensemble simulation of water and energy surface fluxes was produced over a 2530-km² domain in South-West Niger for the 2005 monsoon season and subsequent dry-down (six-month total duration).

A detailed analysis of water—including rainfall P—and energy flux variabilities in the time, space, and ensemble dimensions was performed over a wide range of time and space scales, from the day and kilometer up to the entire simulated time/space domain. Variability in the ensemble dimension reflected the propagation of P uncertainty. Sources of spatial flux variability were rainfall and land use (ecosystem types), those of temporal variability were rainfall, meteorology, and vegetation phenology. Variability was described in the form of both distributions of simulated fluxes and their first-order statistical moments (variances/standard deviations), including global variance analysis of the simulated ensemble set.

The global variance analysis allowed characterizing the respective weights of, and interactions between, the three main sources of variability, namely spatial (from cell to meso) and temporal (daily to seasonal) variability in model forcings (rainfall, meteorology, and ecosystem type) as well as rainfall uncertainty. Interaction terms were found to be especially important, making, for example, spatial variability substantially higher when accounting for full flux uncertainty than when considering ensemble means—i.e., expected fluxes—only, especially at the daily scale. The effect of ecosystem type on flux variability was found to be significant for all flux types, but not generally exceeding the rainfall uncertainty effect. Respective weights of land use versus rainfall in first-order spatial flux variability varied greatly with flux type.

Uncertainty patterns were further investigated using synthetic measures of ensemble variability, taken either at each node of the time–space domain for a given scale or lumped over one or both time–space dimensions at that scale. In addition to a dimensional measure represented by the ensemble standard deviation—denoted as standard uncertainty (USD)—two dimensionless measures were used. Denoted as uncertainty coefficient (UC) and uncertainty fraction (UF), they were obtained by normalizing USD with the quadratic mean flux or with the total standard deviation in the simulated set for that flux type, respectively.

For most flux types, space/time distributions of uncertainty were found to be highly imprinted by the distributions of flux magnitudes, with a generally positive relationship between the two variables. This relationship tended to level off for the largest magnitudes at the daily scale, starting with P.

The relationship of flux uncertainty to P uncertainty was investigated by analyzing the distributions over time or space of UC propagation patterns—from rainfall to the various flux types—using the UC measure lumped over the non-analyzed—time or space—dimension (if applicable). This analysis highlighted general invariance over space in propagation direction (amplification versus attenuation) for all flux types at any timescale, with propagation factors showing less variation within a flux type than between flux types, for whichever daily-scale or time-integrated fluxes. Hence, with the relative exception of Dr that shows some more dispersion, the spatially-distributed UC for any dependent flux type appears to scale smoothly with that of P, irrespective of other local properties (land use and rainfall magnitude). In the time dimension, rainfall uncertainty at a given time logically impacts the various fluxes over a period of time, dependent on flux type. Interpretation of UC time series generally suggests that propagation includes a strong time-lag type of effect with a lag ranging from <1 day for runoff Q to several days for Tp. Only for Dr does the P uncertainty propagate into a completely transformed signal, suggesting large diffusion over time of instant source uncertainty.

Finally, the strong general consistency in distributed uncertainty patterns over time and space for the various flux types allowed us to summarize the effects of time and space scales using the fully domain-lumped, dimensional, and dimensionless uncertainty measures. For USD as well as UC, the ranking of the flux types was more or less scale-independent, with: (i) highest values of USD for P and Q, and lowest for Rn, Tp, and Dr in decreasing order, and (ii) highest values of UC for Dr and Q, lowest for Rn. Both measures display a highly regular scaling behavior for all variables, characterized by a steady decrease with both time and space scales. The scaling in one dimension was suggested to be nearly independent of the other dimension and to closely follow a power law. Dimensionless UC also displays a high degree of scaling invariance between flux types and between the two—time and space—dimensions, with exponents in the range −0.163 to −0.204 over space and −0.172 to −0.332 over time (except for Dr). Exponent values for USD-scaling give information on the degree of time or space autocorrelation in uncertainty on the various flux types, showing a general inverse ranking of flux types for space and time, with P at one end (highest autocorrelation in space, lowest in time) and Dr at the other end (vice-versa). Expectedly, many of the above regularities in scaling patterns were preserved in the uncertainty propagation factors. For the uncertainty coefficient UC, propagation shows amplification for Q (in the ~2–3 range) and—even more so—Dr (~2–10), but attenuation for all other flux types (0.3–0.6 for Tr, Es, LE, and H, and only ~0.1 for Rn). Hence, hydrological applications of land surface modeling appear to be the most exposed to Sahelian rainfall uncertainty, but, even though to a lesser extent, agronomic, ecological, and meteorological applications are also quite sensitive through the still substantial response in the eddy fluxes (Tr, Es, LE, and H). Uncertainty in the energy budget turns out much smaller for the available energy than for its partitioning into the turbulent heat fluxes. The highly-sensitive response of the area’s Hortonian runoff to P uncertainty was previously illustrated using a rainfall–runoff model [46].

Compared to the USD and UC measures, the other dimensionless measure UF, which is directly related to the Ω coefficient of Koster et al. [77], shows a much less regular behavior, in relation to the loss of significance of this measure when scales increase towards the time/space domain size. This suggests that the newly introduced lumped uncertainty coefficient UC could be a more efficient metric across the full range of scales, making a powerful measure of relative uncertainty for the various flux types.

Compared to a more standard sensitivity analysis of rainfall error propagation, the advantage of the stochastic ensemble method used here is that it explicitly accounts for the time and space dependence structure of P uncertainty, thus allowing us to derive global uncertainty estimates and uncertainty properties at any time/space scale included in the study domain. Limitations of this study include the single year of analysis which does not permit evaluating how results may vary over years or decades, as well as the disconnection between the meteorological and rainfall forcings that may induce some inconsistencies during rain events (unusually high evaporation conditions leading in particular to very strong negative H by artefact oasis effect). However, these particular situations are seen as too rare (<1% of time–space domain) to significantly affect our conclusions on uncertainty propagation.

Although there might be some risk in extrapolating the uncertainty propagation and scaling relationships obtained here beyond the levels of rainfall uncertainty and scales encountered in this study or outside the environment explored, the plots obtained do provide a reference for these relationships in the Sahelian context. As a spatial location with respect to the raingauge network stands as a major driver of time-lumped uncertainty of P and then of most dependent fluxes (Figure 10), a relationship could be sought between UCs for these fluxes and a synthetic indicator of this relative location, such as the geometrical component in kriging variance, that would allow generalizing the mapping of UC for the various flux types.

Finally, besides the quantification of uncertainties, an additional, valuable advantage of uncertainty-reflecting ensemble simulation is that it allows for unbiased estimation of flux expectations, by ensemble averaging. Conversely, as illustrated by Vischel et al. [46], using deterministic—even though unbiased—rainfields when these are uncertain, may lead to substantial biases in simulated fluxes, that strongly affect their inferred time/space distributions. These biases arising from the simulation with the unbiased ensemble mean rainfield, as well as those produced with presumably biased alternative, deterministic rainfall interpolation methods, are being quantified, using the ensemble mean fluxes as the reference.