Ph.D Nathalie Christian Saint-Geours

Christian Lavergne

Jean-Stéphane Bailly

Frédéric Grelot

Professor Jean-Stéphane Lavergne

Professor Frédéric Bailly

Grelot

Change of support effects in spatial variance-based sensitivity analysis

Keywords: sensitivity analysis, Sobol' indices, model upscaling, change of support, regularization theory, spatial model Sensitivity analysis, Sobol' indices, Model upscaling, Change of support, Regularization theory, Spatial model

come

Introduction

Variance-based global sensitivity analysis (GSA) is used to study how the variance of the output of a model can be apportioned to different sources of uncertainty in its inputs. Here, the term model denotes any computer code in which a response variable is calculated as a deterministic function of input variables. Originally developed in the 1990s [START_REF] Sobol | Sensitivity analysis for non-linear mathematical model[END_REF], GSA is now recognized as an essential component of model building (European Commission 2009; US Environmental Protection Agency 2009) and is widely used in different fields [START_REF] Cariboni | The role of sensitivity analysis in ecological modelling[END_REF][START_REF] Tarantola | Can global sensitivity analysis steer the implementation of models for environmental assessments and decision-making?[END_REF]. GSA is based on the decomposition of a model output variance into conditional variances. So-called first-order sensitivity indices measure the main effect contribution of each uncertain model input to the model output variance.

Based on these sensitivity indices, ranking the model inputs helps to identify inputs that should be better scrutinized first. Reducing the uncertainty on the inputs with the largest sensitivity indices (e.g., by collecting additional data or changing the geographical pattern of data locations) will often result in a reduction in the variance of the model output. More generally, GSA helps to explore the response surface of a black box computer code and to prioritize the possibly numerous processes that are involved in it.

Although GSA was initially designed for models where both inputs and output can be described as real valued random variables, some recent work has extended GSA to environmental models for which both the inputs and output are spatially distributed over a two-dimensional domain and can be described as random fields [START_REF] Lilburne | Sensitivity analysis of spatial models[END_REF], for a review). In these works, the computer code under study uses maps derived from field data (e.g., digital elevation models and land use maps). These maps are uncertain due to measurement errors, lack of knowledge or aleatory variability [START_REF] Brown | The Data Uncertainty Engine (DUE): a software tool for assessing and simulating uncertain environmental variables[END_REF][START_REF] Refsgaard | Uncertainty in the environmental modelling process: a framework and guidance[END_REF]). The uncertainty of these spatial inputs is usually modeled using random fields. Model output is also spatially distributed (e.g., a flood map or a pollution map). Authors use geostatistical simulation to incorporate spatially distributed model inputs into the GSA approach [START_REF] Ruffo | Hydrocarbon exploration risk evaluation through uncertainty and sensitivity analyses techniques[END_REF][START_REF] Saint-Geours | Is there room to optimise the use of geostatistical simulations for sensitivity analysis of spatially distributed models?[END_REF]) and they display estimation procedures to compute sensitivity indices in a spatial context, either with respect to the spatial average of the model output [START_REF] Lilburne | Sensitivity analysis of spatial models[END_REF] or with respect to the values of the model output at each site of a study area (Heuvelink et al. 2010;[START_REF] Marrel | Global sensitivity analysis for models with spatially dependent outputs[END_REF][START_REF] Pettit | Full-field sensitivity analysis through dimension reduction and probabilistic surrogate models[END_REF].

Nevertheless, to date, none of these studies has reported on a key issue: the link between uncertainty propagation and model upscaling/downscaling.

Model upscaling is the problem of translating knowledge from smaller scales to larger [START_REF] Heuvelink | Uncertainty analysis in environmental modelling under a change of spatial scale[END_REF]. In many environmental models, the physical quantities considered are spatially additive (e.g., porosity or evapotranspiration), i.e., their large-scale properties derive from small-scale properties by simple averaging (Chilès and Delfiner 1999, p.593). In this case, the model end user is usually interested in the spatial linear average or the sum of spatial output over a given spatial unit (e.g., the average porosity of a block or the total evapotranspiration over a plot of land) and model upscaling is thus reduced to a change of support problem (namely, a change of support of the end user's output of interest). [START_REF] Heuvelink | Uncertainty analysis in environmental modelling under a change of spatial scale[END_REF] observed that under a change of spatial support of the model output, the relative contribution of uncertain model inputs to the variance of the aggregated model output may change. Exploring how sensitivity analysis results interact with such a change of support is thus of great importance. It would allow the modeler to check the robustness of model-based environmental impact assessment studies and better assess the confidence of their results. Knowledge of this interaction would also allow the modeler to answer the following questions:

What are the model inputs that explain the largest fraction of the variance of the output over a given spatial support? For which output support size does a given spatially distributed model input contribute to the largest fraction of the variance of the model output? How does the contribution of a spatially distributed input to the variance of the model output depend on the parameters of its covariance function?

The change of support effect has been extensively discussed in geostatistics in the context of regularization theory (Journel and Huijbregts 1978, p.77). Hence, we attempt in this paper to integrate regularization theory with variance-based GSA framework. Our idea is to define site sensitivity indices and block sensitivity indices to i) provide a simple formalism that extends variance-based GSA to spatial models when the modeler's interest is in the spatial average or the sum of model output over a given spatial support (Sect. 2) and ii) discuss how the relative contribution of uncertain model inputs to the variance of model output changes under model upscaling (Sect. 3). We limit our study to point-based models, i.e., models for which the computation of the model output at some location uses the values of spatial inputs at that same location only [START_REF] Heuvelink | Accounting for spatial sampling effects in regional uncertainty propagation analysis[END_REF].

An example is used throughout this paper to illustrate formal definitions and properties. Finally, we discuss the limits of our approach and its connections to related works in Sect. 4.

C(h) =        C(0) if h = 0 (1 -η) • C(0) • ρ a ( h ) if h = 0 (1)
where ρ a (•) is some valid correlogram (Cressie 1993, p.67). The model output is a 2D random field {Y (x) : x ∈ D} that we will simply denote by Y (x). We assume that the first two moments of Y (x) exist. Finally, as discussed in the introduction, we limit our study to point-based models; hence, we assume that there exists a mapping ψ : R n × R → R such that:

∀x ∈ D , Y (x) = ψ [U, Z(x)] (2)
S i = Var [E(Y |X i )] Var (Y ) (3)
S U = Var [E(Y (x * ) | U)] Var [Y (x * )] ; S Z = Var [E(Y (x * ) | Z(x))] Var [Y (x * )] (4)
First-order sensitivity indices of model inputs with respect to the block average Y v are called block sensitivity indices and are denoted by S U (v) and

S Z (v): S U (v) = Var [E(Y v | U)] Var [Y v ] ; S Z (v) = Var [E(Y v | Z(x))] Var [Y v ] (5)
The ratio S Z (v)/S U (v) gives the relative contribution of model inputs Z(x)

and U to the variance of the output of interest

Y v . When S Z (v)/S U (v) is
greater than 1, the variance of Y v is mainly explained by the variability of the 2D input field Z(x); when S Z (v)/S U (v) is less than 1, it is the non spatial input U that accounts for most of Var (Y v ).

Illustrative example

The proposed formalism for spatial GSA is illustrated by the following example. A model Y = M (U, Z) is used for the economic assessment of flood risk over a given floodplain D. Z(x) is the map of maximal water levels (m) reached during a flood event. Z(x) is assumed to be a Gaussian random field with mean µ = 50 and exponential covariance C(h) with C(0) = 100, correlation length a = 5 and nugget parameter η = 0.1. U is a set of three economic parameters U 1 , U 2 and U 3 that determine a so-called damage function that links water levels to monetary costs. U 1 , U 2 and U 3 are assumed to be independent random variables following Gaussian distributions N (1.5, 0.5), N (55, 5) and N (10, 10), respectively. Random field Z(x) and random vector U are supposed to be independent. Model output Y (x) is the map of expected economic damages due to the flood over the area; these damages depend on U and Z(x) through the mapping ψ:

∀x ∈ D, Y (x) = ψ [U, Z(x)] = U 1 • Z(x) -U 2 • e -0.036•Z(x) -U 3 (6)
Stakeholders are interested in two outputs: the flood damage Y (x * ) on a specific building x * ∈ D and the total damage |v| • Y v over a district v (here, a disc of radius r = 50). Here, the expression of mapping ψ and the statistical characterization of model inputs may be simple enough that exact values of sensitivity indices could be derived, but this is usually not the case in real applications in which the model is very complex. A usual alternative is to consider model M as a black box and estimate sensivity indices with Monte-Carlo simulation. We chose to use the estimators and the computational procedure described by Lilburne and Tarantola (2009, Sect. 3.2), based on a quasi-random sampling design, using N = 4096 model runs (Table 1). It appears that on a given site x * , the variability of the water level map explains most of the variance of Y (x * ): S Z = 0.89. On a larger spatial support, the variance of the total flood damage |v| • Y v is mainly due to the economic parameters U 1 , U 2 and U 3 : S U (v) = 0.86. Thus, to improve the accuracy of damage estimation for a specific building, the uncertainty should first be reduced on the water level map Z(x); however, to improve the accuracy of total damage estimation over a large district v, the modeler should focus on reducing the uncertainty of economic parameters U 1 , U 2 and U 3 .

3 Change of support effect on block sensitivity indices

In this section, we assess how the ranking of uncertain model inputs based on their block sensitivity indices vary under a change of support v of model output.

Relation between site sensitivity indices and block sensitivity indices

Site sensitivity indices and block sensitivity indices are related. Let E Z Y (x)

denote the conditional expectation of Y (x) given Z(x), that is: 2)) where f U (•) is the multivariate pdf of random vector U. Under our assumptions concerning Y (x), E Z Y (x) is a second-order SRF. Let C * (•) denote its covariance function, σ 2 = C * (0) its variance and σ 2 v its block variance over support v, that is, the variance of block average

∀x ∈ D, E Z Y (x) = E [Y (x) | Z(x)] (7)
E Z Y (x) is the transform of the input SRF Z(x) via the function ψ(z) = R n ψ(u, z)f U (u)du (Eq. (
1/|v| v E Z Y (x)dx. Block variance σ 2 v is equal to the mean value of C * (h)
when the two extremities of lag vector h describe support v, which we denote by C * (v, v) (Journel and Huijbregts 1978, p.78). Using these notations, it follows from Eq. ( 4) and ( 5) that site sensitivity indices and block sensitivity indices are related by (see Appendix A for a proof):

S Z (v) S U (v) = S Z S U • σ 2 v σ 2 = S Z S U • C * (v, v) C * (0) (8)
3.2 Change of support effect Consider now that model M was initially developed to study the spatial average Y v over the support v, and that after model upscaling the modeler is interested in the spatial average Y V over the support V , where V ≫ v.

We know from Krige's relation (Journel and Huijbregts 1978, p.67) that the block variance σ 2 v decreases with increasing size of support: σ 2 V ≤ σ 2 v . It follows from Eq. ( 8) that:

S Z (V ) S U (V ) ≤ S Z (v) S U (v) (9)
The fraction of the variance of the aggregated model output explained by the input random field Z(x)-compared to the fraction explained by Uis thus smaller on support V than on support v. More specifically, let us suppose that the covariance function C * (•) of the random field E Z Y (x) has a finite effective range and that the support v is large with respect to this range. To a first approximation, the block variance σ 2 v is of the form σ 2 v ≃ σ 2 A/|v|, where A is the so-called integral range of C * (•) and is defined by (Chilès and Delfiner 1999, p.73). It follows from Eq. ( 8) that:

A = 1/σ 2 C * (h)dh
S Z (v) S U (v) ≃ |v| c |v| with |v| c = A • S Z S U (
C * (h) = ∞ k=0 λ 2 k • [C(h)] k (11)
For most of the usual transition covariance functions (e.g., spherical, exponential and Gaussian models), the covariance C(h) is a monotically increasing function of correlation length a. In this case, it follows from Eq. ( 11) that the integral range A = 1/σ 2 C * (h)dh also increases with correlation length a. An increase in correlation length a thus leads to an increase in the critical size |v| c , and the ratio of block sensitivity indices S Z (v) and S U (v) satisfies (Eq. ( 10)):

∂ ∂a S Z (v) S U (v) ≥ 0 (12)
The relative contribution of the uncertain model input Z(x) to the variance of the output of interest Y v increases when the correlation length of Z(x)

increases. Indeed, when correlation length a increases, the averaging-out effect that occurs when the model output is aggregated over spatial support v weakens; thus, the fraction of the output variance Var(Y v ) which is explained by the input random field Z(x) increases.

Nugget parameter's impact on the block sensitivity indices can be interpreted in the same manner. The nugget parameter η controls the relative part of pure noise in the input random field Z(x) (Eq. ( 1)). The smaller η is, the weaker the averaging-out effect will be when the block average Y v is computed over the support v, and the larger the part of output variance Var(Y v ) will be that is explained by Z(x). The critical size |v| c is thus a decreasing function of nugget parameter η, and the ratio of block sensitivity indices S Z (v) and S U (v) satisfies (Eq. ( 1), ( 8), ( 11)):

∂ ∂η S Z (v) S U (v) ≤ 0 (13)

Illustrative example

To illustrate the change of support effects on sensitivity analysis results, we performed spatial GSA on our numerical example in the following settings:

varying disc-shaped support v of increasing size (Fig. 2); varying correlation length from a = 1 to a = 10 (Fig. 3); varying nugget parameter from η = 0 to η = 0.9 (Fig. 4). For each setting, we computed estimates of the output variance Var(Y v ), the block sensitivity indices S U (v), S Z (v) and the ratio

S Z (v)/S U (v) over N = 4096 model runs.
Mean values with a 95% confidence interval were then computed for each estimate using bootstrapping (100 replicas). In accordance with Eq. ( 9), ( 12) and ( 13), it appears that the block sensitivity index S Z (v) (i) decreases when the support v increases It should be noted that our approach is based on conditions that may not be met in some practical cases. First, we considered a model M with a single spatially distributed input Z(x). In real applications, modelers may have to deal with several spatial inputs Z 1 (x), . . . , Z m (x), with different covariance functions C i (•), correlation lengths a i and nugget parameters η i . In this case, it can be shown that Eq. ( 8) still holds separately for each spatial input

Z i (x)
. However, no conclusion can be drawn a priori regarding how a change of support affects the relative ranking of two spatial inputs Z i (x) and Z j (x);

the ratio of their block sensitivity indices S Z i (v)/S Z j (v) will depend on the ratio of block variances σ 2 v,i /σ 2 v,j . Second, some environmental models are not point-based and involve spatial interactions (e.g., erosion and groundwater flow models). In this case, it still may be possible to build a point-based surrogate model as a coarse approximation of the original model; if not, then the change of support properties discussed in Sect. 3 may not hold. Third, we assumed the input random field Z(x) to be stationary; if it is not, site sensitivity indices depend on site x * (Eq. ( 4)). It is then possible to compute maps of these indices [START_REF] Marrel | Global sensitivity analysis for models with spatially dependent outputs[END_REF][START_REF] Pettit | Full-field sensitivity analysis through dimension reduction and probabilistic surrogate models[END_REF] to discuss the spatial variability of model inputs sensitivities.

Finally, we focused on the case in which the modeler's interest is in the spatial linear average or the sum of model output Y (x) over the support v. As discussed by [START_REF] Lilburne | Sensitivity analysis of spatial models[END_REF], other outputs of interest may be considered, such as the maximum value of Y (x) over v (e.g., maximal pollutant concentration over a zone), some quantile of Y (x) over v (Heuvelink et al. 2010), or the percentage of v for which Y (x) exceeds a certain threshold.

To our knowledge, no study has investigated the properties of sensitivity indices computed with respect to such outputs of interest.

Conclusion

This paper provides a formalism to apply variance-based global sensitivity analysis to spatial models when the modeler's interest is in the average or the As mentioned in Sect. 2, we assume that the first two moments of Y (x) exist.

The ratio of block sensitivity indices gives (Eq. ( 5)):

S Z (v) S U (v) = Var (E [Y v | Z(x)]) Var (E [Y v | U]) (14)
The conditional expectation of block average Y v given Z(x) gives:

E [Y v | Z] = E 1/|v| v Y (x)dx | Z(x) (definition of Y v ) = 1/|v| v E [Y (x) | Z(x)] dx (for a point-based model) = 1/|v| v E Z Y (x)dx (definition of E Z Y (x)) Thus we have Var (E [Y v | Z]) = Var 1/|v| v E Z Y (x)dx = σ 2 v (definition of σ 2 v )
. Moreover, the conditional expectation of block average Y v given input

U gives: E [Y v | U] = E 1/|v| v Y (x) dx | U (definition of Y v ) = 1/|v| v E [Y (x) | U] dx (Fubini's theorem) E [Y (x) | U]
does not depend on site x under the stationarity of SRF Z(x);

thus, we have in particular E [Y v | U] = E [Y (x * ) | U], and Var (E [Y v | U]) = Var (E [Y (x * ) | U]
). Combining these expressions with Eq. ( 14) yields:

S Z (v) S U (v) = σ 2 v Var(E [Y (x * ) | U]) (15)
The ratio of site sensitivity indices gives (Eq. ( 4)):

. S Z S U = Var(E [Y (x * ) | Z(x)]) Var(E [Y (x * ) | U]) (16)
We notice that for point-based models

Var [E(Y (x * ) | Z(x)] = Var [E Z Y (x * )] =
σ 2 (definition of E Z Y (x) (Eq. ( 7))). Finally, it follows from Eq. ( 15) and ( 16) that:

S Z (v) S U (v) = S Z S U • σ 2 v σ 2 Appendix B: Hermitian expansion of random field E Z Y (x)
The random field E Z Y (x) can be written (Eq. ( 2), ( 7)) as a transformation of the Gaussian random field Z(x) through the function ψ :

z → R n ψ(u, z) • f U (u) du: E Z Y = ψ (Z)
where f U (•) is the multivariate pdf of random vector U. Under the hypothesis that the first two moments of Y (x) exist, random field E Z Y (x) has finite expected value and finite variance. Thus, ψ belongs to the Hilbert space L 2 (G) of functions φ : R → R, which are square-integrable with respect to Gaussian density g(.). Hence, ψ can be expanded on the sequence of Hermite polynomials (χ k ) k∈N , which forms an orthonormal basis of L 2 (G) (Chilès and Delfiner 1999, p.399):

ψ = ∞ k=0 α k • χ k with χ k (z) = 1 √ k! • 1 g(z) • ∂ k ∂z k g(z)
where coefficients α k are given by:

α k = R χ k (z) ψ(z)g(z) dz. It follows that E Z Y (x)
can be written as an infinite expansion of polynomials of Z(x):

∀x ∈ D, E Z Y (x) = ∞ k=0 α k • χ k [Z(x)]
Its covariance then gives (Chilès and Delfiner 1999, p.396, Eq. (6.23) and p.399, Eq. (6.25)):

Cov (E Z Y (x), E Z Y (x + h)) = ∞ k=0 α 2 k • C(h) C(0) k = ∞ k=0 λ 2 k • [C(h)] k
where C(h) is the covariance function of GRF Z(x) and

λ k = α k • C(0) -k/2 .

Captions of tables

Tab

Response to reviewer's comments

We would like to thank the reviewer for his comments and suggestions. Our comments and answers are provided in the response below (in blue).

-------------

Conditional Expectation

Let f X,Y(x,y) be the joint probability density function for random variables X, Y. To compute E[X |Y ] we must first consider the conditional probability density function for X (given Y = y). This is

f X | Y(x | y) = f X,Y(x,y)/ f Y(y) for f Y(y) ≠ 0 = 0 for f Y(y) = 0 Where f Y(y) = ∫ f X,Y(x,y)dx (the marginal probability density for Y) Then E[X | Y = y] = ∫ xf X | Y(x |y)dx provided ∫| x|f X | Y(x |y)dx is finite.
While some properties of conditional expectation can be presented/discussed without explicitly mentioning the joint probability density it is always implicitly involved. In particular a proof of the "Total Law of Variance", i.e. the decomposition of the variance of X does not require a particular joint probability density function. But to compute Var{E[X |Y]} for a particular pair of random variables does require knowing the joint density function and both marginals.

Var{E[X | Y = y]} = EY{(E[X | Y])2 -(EY{E[X | Y]})2
The subscripting "EY" is a reminder that to compute this expectation the marginal density for Y must be used,

While (EY{E[X | Y]})2 = (E[X])2 , EY{(E[X | Y])2 ≠ E[X2] (unless X, Y are independent)
--> Thank you for this clarification. Some notations and explanations have been modified in the revised manuscript. U and Z(x) are supposed to be independent thoughout the manuscript (we forgot to make it clear in the initial submission). See the answers below for more details.

VARIOUS CONCERNS

• Page 2, lines 9-14 "Variance-based global sensitivity analysis" is an inanimate object or concept, hence it can't "aim" to do anything. The user of this methodology might have a particular goal or object in mind but the method itself does not. While variance is somewhat related to "uncertainty" it is not the same. "Uncertainty" would have to be quantified in terms of probabilities.

--> Thank you. We replaced "aim" by "is used to" in the revised manuscript. We also changed most occurrences of the term "uncertainty" --that was indeed too vague --by "variance" when it was appropriate.

• Page 3, line 29 The referenced document from EPA is no longer available at that URL. Elsewhere on the EPA site it indicates that this was a document put out to solicit public comment. It did not represent adopted policy or guidance.

--> Thank you. On April, 2012, 12th, we found the referenced document still available on the Web, at the URL which is given in the manuscript: http://www.epa.gov/crem/library/cred_guidance_0309.pdf • Page 3, lines 45-47 Since the authors are using the terms "model", "input", "output" in only an abstract, vague way there is no assurance that collecting additional data will reduce the fraction of the variance of the output that is attributable to a particular input. It may also be a question of the nature of the"data". Since the authors are specifically concerned with spatial models, the geographical pattern of the data locations may be very important and changing that pattern may result in a reduction in the variance of the output.

--> Thank you. Indeed, there is no assurance that collecting additional data will reduce the fraction of the variance of the output that is attributable to a particular uncertain input --even if it is the case in many applications. The question of the spatial pattern of the data locations is also very important, and exploring how sensitivity indices vary under a change of these data locations would be interesting.

Yet our goal in this paragraph was just to give a general --and probably too vague as you mentioned it -overview of what sensitivity analysis is and what it is usually used for. In most cases, sensitivity analysis is used to identify the main processes involved in a complex computer code, and to select the inputs that contribute the most to the variance of the model outuput. Then the modeller will try to get a more accurate estimation of the 'real' value of these inputs. In most cases (yet not always) reducing the variance of these uncertain inputs will result in a reduction of the variance of the model output. We changed this sentence in the revised manuscript, to try and be more precise:

"Based on these sensitivity indices, ranking the model inputs helps to identify inputs that should be better scrutinized first. Reducing the uncertainty on the inputs with the largest sensitivity indices (e.g, by collecting additional data or changing the geographical pattern of data locations) will often result in a reduction in the variance of the model output."

• Page 3, line 50 "influence" is a vague term. Presumably the authors mean that the presence/absence of the inputs in question don't cause the variance of the output to increase or decrease. They should be much more precise.

--> Thank you. Here we wanted to explain that if a model input has a very small sensitivity index compared to the other inputs, then trying to reduce the uncertainty on this input is of small interest : in any case it wouldn't cause the variance of the output to decrease much. As it was not clear enough and not really important, we removed this sentence in the revised manuscript.

• Page 3, line 52 "Model behavior" is a very vague term and can have a variety of meanings. "Understanding behavior" is at least partially a characteristic of the person(s) using the model.

--> Thank you. Our sentence was too vague. What we tried to say briefly is the following: global sensitivity analysis (GSA) is useful to explore the "response surface" of a complex computer code with many inputs, no simple analytical expression, that may represent various physical or anthropic processes. By performing GSA, the modeller will usually get a better idea of what are the main processes in his model, which processes strongly interact with others, which process may be simplified or even removed from the model, etc. We changed this sentence in the revised manuscript :

"More generally, GSA helps to explore the response surface of a black box computer code and to prioritize the possibly numerous processes that are involved in it."

• Page 3, line 56. Presumably in using the term "scalar inputs" the authors are contrasting those with "vector inputs". In the example on page 7 the inputs are either random variables or a random field hence the term "scalar inputs" is not appropriate. However the inputs in this example are presumably "scalar valued", i.e. each of the random variables is real valued and the random field is real valued. Also see line 9 on page 7 --> Thank you. What we wanted to oppose is: on one side, uncertain inputs that are modeled as real valued random variables; on the other side, uncertain inputs that are modeled as random functions or random fields (real valued or not). We changed the sentence in the revised manuscript:

"Although GSA was initially designed for models where both inputs and output can be described as real valued random variables, some recent work has extended GSA to environmental models for which both the inputs and output are spatially distributed over a two-dimensional domain and can be described as random fields."

We also changed the description of inputs U and Z(x) on section 2.1.

• Page 4, line 14 Delete "sampled". The field data is likely already a sample. Note that digital elevation models and land use maps might refer to physical presentations, e.g. printed on paper or they might refer to a visual representation in electronic form. In either case the maps are not the actual inputs.

To actually determine/choose a random field model will likely require collecting data and then estimating parameters. Moreover there are different types of random field models and those different types might result in different variances for the output, this could result in a significant amount of "uncertainty" that is not quantified.

-->Thank you. We removed the word "sampled" in the revised manuscript, and rephrased the whole sentence:

"In these works, the computer code under study uses maps derived from field data (e.g., digital elevation models and land use maps)."

We couldn't make a clear distinction in the manuscript (to try keep things simple and short enough) between the actual computer code inputs (which may be numbers or maps), and the description of the uncertainty on these inputs (using random variables or random fields).

From a very practical point of view, maps (digital elevation model, landuse map) are stored as layers in a GIS software. These layers are the actual inputs that the computer code uses to calculate its output. Then, for real valued inputs, uncertainty is described by representing the input with a real valued random variable with a given pdf. For spatially distributed inputs (maps), the uncertainty can be described by representing the input with a given random field model. In this case, describing properly the random field model will indeed require collecting data and estimating parameters. The choice of the random field model will also result in additional "uncertainty" that is not quantified.

We rephrased section section 2.1 in the revised manuscript to make a clearer distinction between these two different notions.

• Page 4, line 35 Does "efficient estimation procedures" pertain to the meaning used in the statistical literature or does it have some other meaning? In the statistical literature one might refer to an "efficient estimator", i.e., one which has a smaller variance than other estimators but that would not make sense with respect to "procedures".

--> Thank you. We misused the word "efficient" in this sentence. We didn't intend to talk about the "efficiency" of estimators with the meaning used in the statistical litterature. Thus we just removed the word "efficient" from this sentence in the revised manuscript.

• Page 4, line 54 Note that the term "scale" has multiple (and sometimes contradictory) meanings.

That is one of the reasons why the term "support" is used in the geostatistical literature as opposed to scale.

--> Thank you. We replaced most occurrences of the term "scale" by the term "support" in the revised manuscript. We kept the terms "scale" and "model upscaling/dowscaling" on p.4, line 54, because we are citing a paper (Heuvelink 1998, p.1) and we don't want to modify the expression used by the original author. We also kept the terms "small-scale properties" and "large-scale properties" on p.4, line 55, because we are citing a definition from Chilès & Delfiner (1999, section 8.1, p.593) that we don't want to modify either.

• Page 5, lines 43 & 53. What do "influential" and "influence" mean in this context? Also see line 22 on page 6 and many other places. "Influence" could mean deterministic, i.e. actual cause-effect or simply correlated or some other vague idea. The authors are using the wrong word --> Thank you. Manuscript was modified according to this recommendation: we removed all occurrences of the terms "influential" and "influence" in the revised manuscript.

• Page 5, lines 48-50. Change"most to model output uncertainty" to "largest fraction of the variance of the output of the model"

• Page 5, line 50 What is "spatial structure of uncertainty of a model input"? Moreover the authors never actually consider "uncertainty" they only consider the variance of the model output. All references to "uncertainty" should be omitted or appropriately changed.

--> Thank you. Our sentence was too vague. What we wanted to denote by the expression "spatial structure of uncertainty of a model input" is more precisely the shape and parameters of the covariance function C(.) of the random field model Z(x) that is used to describe an uncertain spatial input. We modified this expression in the revised manuscript to be more precise:

"How does the contribution of a spatial input to the variance of the model output depend on the parameters of its covariance function?"

Moreover, we replaced all occurences of the term "uncertainty" by "variance" when it was appropriate to do so (not in the very begining of the introduction which we want to be very general).

• Page 6, lines 27-33. The description of "point-based" models is incorrect.

--> Thank you. We took the term "point-based model" from the following paper (p.1): Heuvelink GBM, Brus DJ, Reinds G ( 2010) Accounting for spatial sampling effects in regional uncertainty propagation analysis. In: Tate NJ, Fisher PF (eds) Proceedings of the 9th international symposium on spatial accuracy assessment in natural resources and environmental sciences (Accuracy2010), pp 85--88

where the authors say: "Many environmental models involve spatial interactions. Examples are erosion, groundwater flow and plant dispersal models. However, there are also many environmental models that are essentially point-based. For instance, models that predict crop growth, greenhouse gas emission, soil acidification or evapotranspiration at some location typically use soil, landuse, management and climate input data at that same location only."

We want to limit our study to models that can be described by Eq.( 2) in our manuscript:

for all x in D, Y(x) = f [ U, Z(x) ]
We don't want to discuss the case where Y(x*) (model output at some location x*) would be a function of Z(x*) (model input at that same location) and Z(y) with y ≠ x* (model input at some other location). Thus we want to limit our study to the kind of models where "the computation of the model output at some location x* uses the values of spatial inputs at that same location x* only".

We don't know if that family of models has a name in the geostatistical litterature. We would be happy to change the term "point-based" to a more appropriate term if it exists. As we don't know which other term to use so far, we didn't modify this term in the revised manuscript.

• Page 7, lines 14-17. Although there is some variation in the statistical literature of the definition of strict stationarity, it always includes the condition that the joint distribution of Z(x1), …, Z(xn), for any finite collection of points in D, is translation invariant. It may or may not include the requirement that Z(x) have finite expected value and finite variance. If it does then it is also second order stationary. In most applications the statistical characteristics of the random field will not be known, e.g. the multivariate probability distribution and its parameters. Since it is unlikely that one can collect data from multiple realizations of the random field, one can't use the data to determine characteristics of the multivariate probability distribution.

--> Thank you for this clarification. We want to use the notion of "second-order stationary random field" as it is defined in [START_REF] Chilès | Geostatistics, modeling spatial uncertainty[END_REF].17 section 1.1.4.). If we understood properly:

-"strict stationarity" with finite first-order and second-order moments implies "second-order stationarity" -"second-order stationarity" does not necessarily implies "strict stationarity" under the hypothesis of second-order stationarity, it is possible to use the data of a single realization of the random field to determine its mean and the characteristics of its covariance function. Without this hypothesis, one can't use the data to determine characteristics of the multivariate probability distribution. We rephrased the whole section 2.1 in the revised manuscript.

• Page 7, lines 17-49. Since Y(x) is purportedly a function of U and Z(x), in order to consider statistical properties of Y(x) one must assume that U and Z(x) are defined on the same probability space. Are the authors assuming that U and Z(x) are independent or is there some unstated joint probability distribution? In any case, if Y(x) has finite expected value then both of the conditional expectations must exist. To actually compute any of these expectations one needs to know the probability densities and the conditional densities. All of this material needs to be re-written in a more careful way (and statistically correct).

--> Thank you for this question. Random vector U and random field Z(x) are supposed to be independent throughout the paper, but we unfortunately forgot to make it appear in the initial manuscript. We added the following sentence in the revised manuscript (section 2.1):

"U and Z(x) are supposed to be independent."

Under this hypothesis, the conditional density of U given Z(x) is the marginal density for U and the conditional density of Z(x) given U is the marginal density for Z(x).

Also, when we wrote the initial manuscript, we were not sure that when Y(x) has finite expected value and variance, then E[ Y(x) | Z(x)] and E[Y(x) | U] (which are both random variables) necessarily both have finite expected value and finite variance. After some more careful examination, it appears that this property simply derives from conditional Cauchy-Schwarz inequality. Thus we removed the additionnal conditions on E[ Y(x) | Z(x)] and E[Y(x) | U] in the revised manuscript.

Finally, we are not really interested in discussing the expression of the conditional densities of Y(x) knowing U or Z(x), so we don't think there is a need to explicitly mention these conditional densities in the manuscript.

• Page 7, line 17. Change "covariogram" to "covariance function", see page 70 of Chiles and Delfiner for the definition of "covariogram" which is not the same as "covariance function". Also see the same problem on page 12, line 28 and page 13, line 37, page 14, line 54.

--> Thank you. Manuscript was modified according to this recommendation: all occurrences of the term "covariogram" were replaced by "covariance function".

For information, we took the definition of "covariogram" from Cressie (1993, p.53, eqn. 2.3.3) : " The function C(.) is called a covariogram or a stationary covariance function: cov(Z(s1),Z(s2)) = C(s1 -s2) for all s1, s2 in D" and also p.67: " Call the function C(.) given by (2.3.3) (provided it is well defined), a covariogram. (Notice that it has also been called an autocovariance function by time-series analysts.)"

• Page 8, line 31. It is essential to restrict the possible choices for the function f, for example if f(X1, X2) = X1/X2 where X1, X2 are independent standard normal random variables then f(X1, X2) does NOT have either an expected value nor a variance. Thus all of the discussion that follows is erroneous.

--> Thank you. The discussion holds if Y has finite expected value and finite variance. The manuscript was modified according to this recommendation:

"Let us consider a model Y = f(X 1 , ..., X n ) where X i are independent random variables and where the first two moments of Y exist."

• Page 8, lines 43-45. The statement "It is the expected part of the output variance Var(Y) that could be reduced" is incorrect The sensitivity index is really based on the following Var(Y) = E(Var (Y | X)) + Var(E(Y |X)) The two terms on the right hand side are the "fraction of variance unexplained" and "fraction of variance explained". This relationship is sometimes called the "Total law of variance". See most text books on probability theory for more details.

--> Thank you. We took most of the material of this paragraph from the following book, which is very often cited if the field of sensitivity analysis:

Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, Saisana M, Tarantola S, eds (2008) Global sensitivity analysis, the primer. Wiley, New York

In the literature of sensitivity analysis, the definition of sensivity indices is usually built from the decomposition of the function f( .) into what is often called "high dimensionnal model representation" (HDMR) (that was originally suggested by [START_REF] Hoeffding | A class of statistics with asymptotically normal distribution[END_REF], see also [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF][START_REF] Sobol | Theorems and examples on high dimensional model representation Reliability Engineering & System Safety[END_REF], references below):

f(x 1 ,...x n ) = f 0 + f 1 (x 1 ) + ... + f n (x n ) + f 1,2 (x 1 ,x 2 ) + ... + f 1,...n (x 1 ,....,x n )
where the integral of each function f i1,...,ik () is equal to 0. The decomposition of the variance of f(X) that derives from this functional decomposition is very close to the ANOVA framework:

Var[ f(X) ]= Var[ f 1 (X 1 ) ] + ... + Var[ f n (X n ) ] + Var[ f 1,2 (X 1 ,X 2 ) ] + ... + Var[ f 1,...n (X 1 ,....,X n ) ].
First-order sensitivity index of input X i is then defined as the ratio Var[ f i (X i ) ] / Var [ f(X) ], which is the same as the Equation (3) in our manuscript, written in terms of conditional expectations.

Higher-order sensitivity indices can also be defined from this decomposition, to characterize the contribution of the interactions between several uncertain inputs X i1 ,...,X ik to the variance of the output f(X).

The statement "It is the expected part of the output variance Var(Y) that could be reduced" may be incorrect.

In the revised manuscript, we just removed it.

• Page 9, lines 47-50. See the comment just above.

--> OK.

• Page 9, eqs 4 & 5. Delete the ":x' ε D", that is logically incorrect --> Thank you, the manuscript was modified according to this recommendation.

• Page 10, lines 24-26. The exponential covariance function does not have a true range in the sense that there is no spatial dependence beyond that distance. However in the geostatistical literature it is common to refer to the "effective range", i.e. the distance at which the value of covariance function has decreased to .05C(0). If a is the parameter in the functional form of the covariance function then it is not the effective range, however it is what is called the correlation length. The authors need to clarify which it is, "effective range" or "correlation length/integral range". In any case it appears that the authors have chosen a range that is quite small in comparison to the scale of the "block" used later.

--> Thank you. We misused the term "range", denoted by 'a', to refer to the parameter in the functional form of the covariance function C(.). In our illustrative example, the "effective range" is thus close to 3*a (as we use an exponential covariance function:

C(h) = C(0) . [ η δ(|h|) + (1-η)exp(-|h|/a) )]
We changed the manuscript according to the recommendation, replacing all occurrences of the term "range" by "correlation length" or "effective range" depending on the intended meaning.

The term "integral range", denoted by A in our manuscript, was taken from Chilès and Delfiner (1999, p.73):

A = 1/σ² ∫ C(h) dh
where C(.) is a covariance function. If we understood properly this definition, for a given covariance function C(.), the integral range A is not necessarily equal to the parameter a in the functional form of C(.) (that we now call "correlation length"). If C(.) is an exponential covariance function with correlation length a, then the integral range is equal to A = π a². We didn't change the occurrences of the term "integral range" in the revised manuscript.

In the illustrative example, we chose a correlation length varying from a=1 to a=10, corresponding to an effective range varying from 3 to 30. The radius of the disc-shaped "block" was varying from r = 1 to r = 50.

The ratio between the range and the size of the block is of the same order of magnitude than in the real case-study that we work on (not described in the manuscript), where Z(x) is a map of water levels and M is a model that computes total flood damage over a large study area.

• Page 10, section 2.3. Are the random variables and the random field dependent or independent? They have not sufficiently characterized the random field Y(x)

--> Thank you for this question. Random vector U and random field Z(x) are supposed to be independent throughout the paper, but we unfortunately forgot to make it appear in the manuscript. We added the following sentence in the revised manuscript (section 2.1):

"U and Z(x) are supposed to be independent."

• Page 10, line 29 Change "describe" to "determine", very different meanings --> Thank you. Manuscript was modified according to this recommendation.

• Page 10, line 31 Change "expected monetized costs" to "monetary costs". The "expected" term is clearly incorrect here.

--> Thank you. Manuscript was modified according to this recommendation.

• Page 10, lines 53-55. What does "analytical expression" mean here? Does it only mean the formula given in eq(6) or does it mean a complete statistical characterization? In any case the authors have left out the key question of whether the two inputs are dependent or independent? If they are dependent then it is essential to give the actual joint probability distribution. This would be necessary even if Monte Carlo simulation is used.

--> Thank you for this question. By "analytical expression" we mean eq(6) + the knowledge of the pdf of random variables U1, U2, and U3 + the knowledge of random field model of Z(x). Random vector U and random field Z(x) are supposed to be independent throughout the paper. In the revised manuscript we changed this sentence to:

"Here, the expression of mapping ψ and the statistical characterization of model inputs may be simple enough that exact values of sensitivity indices could be derived"

• Page 10, line 60. Change "model complexity is high" to "model is very complex" --> Thank you. Manuscript was modified according to this recommendation.

• Page 11. line 9. The authors should provide more information about how they used a "Monte Carlo approach" and also the "sampling-based method". [START_REF] Lilburne | Sensitivity analysis of spatial models[END_REF] do not use the terminology "sampling-based method", they do consider various "sampling designs".

--> Thank you. We estimated sensitivity indices using the computational procedure described in section 3.2 of [START_REF] Lilburne | Sensitivity analysis of spatial models[END_REF], with a quasi-random sampling design, and with the estimators given in formulas ( 8), (9) (11)-( 16). We changed the sentence in the revised manuscript:

"A usual alternative is to consider model M as a black box and estimate sensivity indices with Monte-Carlo simulation. We chose to use the estimators and the computational procedure described by Lilburne and Tarantola (2009, Sect. 3.2), based on a quasi-random sampling design, using N = 4096 model runs."

• Page 12, line 26 The expression apparently refers back to page 7 but "p(u)" is not identified either on page 7 or on this page (nor on page 25). The integral should be ∫y f Y | Z (y |z) dy where f Y | Z (y |z) = f Y , Z (y , z)/ f Z (z), i.e. the conditional density of Y given Z is the quotient of the joint density of Y, Z and the marginal density of Z. To determine the joint density of Z, Y one must know the joint density of U and Z. The authors have not given that density nor said anything about the joint distribution of U and Z. It is also possible to use a change of variables in the integral to write the above integral in terms of the conditional density of U given Z. Perhaps the authors are assuming (but without ever saying so) that U and Z are independent. In that case the joint density would be the product and then the conditional density of U given Z would be the marginal density for U.

--> Thank you for this clarification. p(u) is supposed to denote the multivariate pdf of the random vector U, as it was (not very clearly) said on page 7, line 12 in the initial manuscript.

Random vector U and random field Z(x) are supposed to be independent throughout the paper, but we unfortunately forgot to make it appear in the initial manuscript. We added the following sentence in the revised manuscript (section 2.1):

"U and Z(x) are supposed to be independent."

Under this hypothesis, the conditional density of U given Z is the marginal density for U, i.e. the multivariate pdf p(u). To make notations more explicit, we replaced the notation p(u) by f U (u) in the revised manuscript.

We also put the definition of f U (u) on p.12 in the revised manuscript: 2)] ---where f U (.) is the multivariate pdf of random vector U."

"E Z Y(x) is the transform of the input SRF Z(x) via the function Ψ(z)= ∫ Ψ(u,z) f U (u) du [Eq. (
• Page 13, line 32 Again the term "influence", also the comparative word should be "lesser", Also see page 14, line 22., page 14, line 30, page 14, line 54. Also see the title of section 3.3, page 15 line 39 --> Thank you. Manuscript was modified according to this recommendation: we replaced the term "influence" by more appropriate terms throughout the manuscript. We also replaced "lower" by "lesser".

• Page 13, line 35 Again "uncertainty" when the authors really mean variance of the output, also see line 9 page 14., line 19 page 14 --> Thank you. Manuscript was modified according to this recommendation: we replaced the term "uncertainty" by "variance" when it was more appropriate.

• Page 13, line 55 Change "drives" to "determines" --> Thank you. Manuscript was modified according to this recommendation.

• Page 14, change "low" to "small" --> Thank you. Manuscript was modified according to this recommendation.

• Page 14, lines 49-54. The real reason that EZ[Y(x] can be represented as an infinite series of Hermite polynomials is that EZ[Y(x] is square integrable with respect to the standard normal density function. Note this is "mean-square convergence, not "point-wise" --> Thank you for this clarification. This is the proof that we try to develop in more details in Appendix B. The sentence page 14, lines 49-54 has been modified to make this reason clearer:

"EZ[Y(x)] is then square-integrable with respect to the standard normal density function. It can be decomposed into an Hermitian decomposition and its covariogram can be written as: (...)"

• Page 15, line 22 Change "verifies" to "satisfies", also see page 16, line 9 --> Thank you. Manuscript was modified according to this recommendation.

• Page 15, line 50 Change "lower" to "smaller" --> Thank you. Manuscript was modified according to this recommendation.

• Page 16, line 40. See comment with respect to page 29 concerning the confidence intervals --> Thank you. Confidence intervals were computed using bootstraping, with n=100 replicas. A clarification was added on this point in the revised manuscript: "Mean values with a 95% confidence interval were then computed for each estimate using bootstrapping (100 replicas)."

• Page 17, line 12 "uncertainty" again --> Thank you. Manuscript was modified according to this recommendation: we replaced the term "uncertainty" by "variance".

• Page 17, line 14 Change "illustrates" to "shows" --> Thank you. Manuscript was modified according to this recommendation.

• Page 17, line 25. The R software system incorporates a great many "packages", some are part of the standard download but others are not. The authors should be more explicit about how they used R, i.e., which "packages" and which functions within those packages..

--> Thank you for this question. We used the GaussRF() function from the RandomFields package to generate realizations of random field Z(x). We also modified the function sobol() from package sensitivity and the S3 method "tell" for objects of class "sobol", to estimate variance-based sensitivity indices with the estimators given by Lilburne and Tarantola (2009, section 3.2, eq. ( 8), ( 9), ( 11)-( 16)). To be more explicit, we added a sentence in the revised manuscript:

"All calculations and figures were realized in R: random realizations of Z(x) were generated with the 'GaussRF()' function from the 'RandomFields' package [START_REF] Schlather | Simulation of stationary and isotropic random fields[END_REF], while computation of sensitivity indices was based on a modified version of the 'sobol()' function from the 'sensitivity' package."

We also added a reference on the RandomFields package [START_REF] Schlather | Simulation of stationary and isotropic random fields[END_REF] in the list of references.

• Page 17, line 42 Change "previous related works" to "various prior publications" --> Thank you. Manuscript was modified according to this recommendation.

• Page 18, line 14 "influence", also line 45 --> Thank you. Manuscript was modified according to this recommendation: we replaced the term "influence" by more appropriate terms throughout the manuscript.

• Page 18, lines 53 & 58 How is "accuracy" quantified? In what sense is accuracy "appropriate"or not appropriate?

--> Thank you for these questions. By "accuracy", we mean a measure of how much the input maps (e.g. a digital elevation model, a landuse map) that are used by the computer code under study, differ from the "real" spatially distributed physical variable they represent. Measuring this "accuracy" is part of the quality assessment of the input map. It could be measured for example from an extra set of validation field points, by computing the root mean square error between the value of the physical variable at validation field points and the value of the variable taken from the input map at the same locations.

The precise meaning of the term "appropriate" depends on the goal of the modeller. In a given context, the modeller will usually fix a limit on the variance of the output of his computer code: under this threshold, he will consider that the model output is "precise enough" for what it will be used for (decision making for example). To reach such a variance, the modeller may have to improve his knowledge on the 'real' value of some the model inputs. In that case, he'd better start working on the model inputs with the largest sensitivity indices, because reducing the variance of these inputs will likely reduce the variance of the output by a large fraction.

We rephrased this sentence in the revised manuscript:

"Our contribution also promotes an increased awareness of the issue of sharing out efficiently, among the various inputs used by a complex computer code, the cost of collecting field data. At some point of the model building process, the modeller will usually aim at reducing the variance of the output below a given threshold, that will depend on the model use. To do so, the modeller may have to improve his knowledge on the real value of some of the model inputs, usually by collecting extra data. In this case, gathering extra field data on inputs maps that have small sensitivity indices (S Z (v) < 0.1) would be unefficient, as it would be costly but could not reduce the variance of the model output by a large fraction."

• Page 24, lines 31 & 34 Delete the ":x' ε D", that is logically incorrect. The block average is not conditioned on the random field but rather on a random variable, i.e. Z at a particular location. Also see eqs 4 & 5 on page 9

--> Thank you, the manuscript was modified according to this recommendation.

• Page 29 In the captions for Figures 3a, 4a it refers to "Yv" but on the figures themselves the titles seem to use a different identification --> Both Figures 3a and4a show the total variance of Y v on the y-axis. The label on the figures was "Variance of Y v ", while the two captions said "total variance of Y v" . To make things clearer, we changed the labels of the y-axis on the figures to "Var[Y v ]".

• Page 29 In the captions for Figures 2,3,4 the authors claim that the error bars show 95% confidence intervals. How are the confidence intervals computed, did the authors use "bootstrapping" or did they assume a normal distribution to do so but with no justification. In any case they should provide clarification.

--> Confidence intervals were computed using bootstrapping, with n=100 replicas. A clarification was added on this point in the captions for Figures 2,3,4: "Error bars show 95 % confidence interval computed by bootstrapping (100 replicas)"

2

  Variance-based sensitivity indices for a spatial model 2.1 Description of spatial model M We want to study a computer code M whose output is a map and whose inputs are a map and a set of n real valued variables. Both inputs and output are uncertain and are described as random variables or random fields. More precisely, we use the following notations: let D ⊂ R 2 denote a 2D spatial domain, x ∈ D a site, h the lag vector between two sites x and x ′ , and v ⊂ D some spatial support (block) of area |v|. We consider the model Y = M(U, Z) where U = (U 1 , . . . , U n ) is a random vector and {Z(x) : x ∈ D} is a second-order stationary random field (SRF) -that we will often simply denote by Z(x). U and Z(x) are supposed to be independent. Covariance function C(•) of Z(x) is assumed to be isotropic, characterized by correlation length a ∈ R, nugget parameter η ∈ [0; 1[ and of the form:

A

  sensitivity analysis of the model M must be performed with respect to a scalar quantity of interest derived from spatially distributed model output Y (x). Here, we consider two different outputs of interest: the value Y (x * ) at some specific site x * ∈ D and the aggregated value Y v = 1/|v| v Y (x) dx over support v. Because model inputs U and Z(x) are uncertain, Y (x * ) and Y v are both random variables; the sensitivity analysis will describe the relative contribution of uncertain model inputs U and Z(x) to their respective variances. 2.2 Site sensitivity indices and block sensitivity indices Before defining sensitivity indices for spatial model M, we briefly review the mathematical basis of variance-based GSA. Let us consider a model Y = G(X 1 , . . . , X n ), where X i are independent random variables and where the first two moments of Y exist. The first-order sensitivity index S i of model input X i is defined by:

S

  i ∈ [0; 1] measures the main effect contribution of the uncertain model input X i to the variance of model output Y . Sensitivity indices can be used to identify the model inputs that account for most of the variance of the model output (model inputs X i with high first-order indices S i ). Sum of S i is always less than 1 and the difference 1i S i accounts for the contribution of the interactions between model inputs X i to model output variance Var(Y ). Please refer to Saltelli et al. (2008) for more details on GSA theory and on the estimation of sensitivity indices. To extend GSA to spatial model M, we propose to use different types of sensitivity indices to describe the relative contribution of the uncertain model inputs U and Z(x) to the variance of the model output: an index on a point support (i.e., with respect to output of interest Y (x * )) and an index on a larger support (i.e., with respect to output of interest Y v ). First-order sensitivity indices of model inputs with respect to Y (x * ) are called site sensitivity indices. Under the stationary hypothesis on SRF Z(x), these indices do not depend on site x * and thus will simply be denoted by S U and S Z :

(

  Fig. 2(b)), (ii) increases with the correlation length a (Fig. 3(b)), and (iii) decreases with the nugget parameter η (Fig. 4(b)). The opposite trends are observed for sensitivity index S U (v). The change of support effect is clearly highlighted in Fig. 2(b): the contribution of the economic parameters U 1 , U 2 and U 3 to the variance of total flood damage |v| • Y v exceeds the contribution of the water level map Z(x) when the radius r of v is greater than r c ≃ 18; for radius r < r c , the variance of total flood damage over the support v is mainly explained by the variability of the water levels Z(x). Finally, Figure 2(c) shows that the ratio S Z (v)/S U (v) is proportional to 1/|v| when the support v is large enough. The theoretical curve S Z (v)/S U (v) = |v| c /|v| (Eq. (10)) was fitted (least squares -R 2 = 0.99) on data points (for r ≥ 20 only), yielding an estimate of the critical size |v| c ≃ 1, 068. All calculations and figures were realized in R (R Development Core Team 2009): random realizations of Z(x) were generated with the GaussRF() function from the RandomFields package (Schlather 2001), while computation of sensitivity indices was based on a modified version of the sobol() function from the sensitivity package.4 DiscussionOur first goal was to provide a formalism that extends the variance-based GSA approach to spatial models when the modeler is mainly interested in the linear average or the sum of a point-based model output Y (x) over some spatial unit v. Our approach is strongly motivated by various prior publications. Other authors had already computed site sensitivity indices[START_REF] Marrel | Global sensitivity analysis for models with spatially dependent outputs[END_REF][START_REF] Pettit | Full-field sensitivity analysis through dimension reduction and probabilistic surrogate models[END_REF] and block sensitivity indices[START_REF] Lilburne | Sensitivity analysis of spatial models[END_REF], but did so without naming them or exploring their analytical properties or their relationship. Our work is an attempt to do so. Equation (8) provides an exact relation between the site and block sensitivity indices, it may prove useful in the case of a model with a simple enough analytical expression.Our research also sought to account for the change of support effects in the propagation of uncertainty through spatial models, within a variance-based GSA framework. We proved that the fraction of the variance of the model output that is explained by a spatially distributed model input Z(x) decreases under model upscaling; when the support v is large enough, the ratio of the block sensitivity index of spatially distributed input to the block sensitivity index of non-spatial inputs is proportional to |v| c /|v|. The critical size |v| c depends on the covariance function of the input SRF Z(x); it usually increases with an increase of the correlation length a or a decrease of the nugget parameter η. These findings are a translation into GSA formalism of the averaging-out effect clearly exhibited by[START_REF] Journel | Mining geostatistics[END_REF] in the regularization theory. Our contribution is to discuss this issue from the point of view of GSA practitioners. Formalizing the effect of a change of support on sensitivity analysis results may help modelers when they consider model upscaling; it will orientate future data gathering by identifying model inputs that will explain the largest fraction of the variance of the model output over a new spatial support. Our contribution also promotes an increased awareness of the issue of sharing out efficiently, among the various inputs used by a complex computer code, the cost of collecting field data. At some point of the model building process, the modeller will usually aim at reducing the variance of the output below a given threshold, that will depend on the model use. To do so, the modeller may have to improve his knowledge on the real value of some of the model inputs, usually by collecting extra data. In this case, gathering extra field data on inputs maps that have small sensitivity indices (S Z (v) < 0.1) would be unefficient, as it would be costly but could not reduce the variance of the model output by a large fraction.[START_REF] Saint-Geours | Analyse de sensibilité de Sobol d'un modèle spatialisé pour l'évaluation économique du risque d'inondation[END_REF] discuss this issue on a flood risk assessment case study.

  sum of the model output Y (x) over a given spatial unit v. Site sensitivity indices and block sensitivity indices allow us to discuss how a change of support modifies the relative contribution of uncertain model inputs to the variance of the output of interest. We demonstrate an analytical relationship between these two types of sensitivity indices. Our results show that the block sensitivity index of an input random field Z(x) increases with the ratio |v| c /|v|, where |v| is the area of the spatial support v and the critical size |v| c depends on the covariance function of Z(x). Our formalization is made with a view toward promoting the use of sensitivity analysis in model-based spatial decision support systems. Nevertheless, further research is needed to explore the case of non-point-based models and extend our study to outputs of interest other than the average value of model output over support v. Appendix A: Proof of the relation between site sensitivity indices and block sensitivity indices

Fig. 1

 1 Fig. 1 Spatial model with uncertain inputs U and Z(x) and spatial output

Figure 1

 1 Figure 1Click here to download high resolution image

Figure 2

 2 Figure 2 Click here to download high resolution image

  . 1 Sensitivity analysis results over N = 4096 model runs with respect to the outputs of interest Y (x * ) and |v| • Y v . Mean values with ±s.d. computed by bootstrapping (100 replicas).

		TABLE 1
	Support	Site x *	Block v
	Output of interest	Y (x * )	|v| • Y v
	Mean of output	66.5 ± 4.2	539 • 10 3 ± 3.6 • 10 3
	Variance of output	1393 ± 188	9 • 10 9 ± 0.2 • 10 9

  impose legally binding requirements; depending on the circumstances, it may not apply to a particular situation. The U.S. Environmental Protection Agency (EPA) retains the discretion to adopt, on a case-bycase basis, approaches that differ from this guidance." (extract from the document disclaimer)

	This document is a final version approved by EPA, after they collected public comments:

"A draft version of this document was reviewed by an independent panel of experts established by EPA's Science Advisory Board and revised by CREM in response to the panel's comments. This final document is available in printed and electronic form." (extract from the document preface) The referenced document does not represent any legal requirement, yet it provides official guidance to develop and evaluate environmental models: "This document provides guidance to those who develop, evaluate, and apply environmental models. It does *Response to Reviewer Comments not


Nathalie Saint-Geours

Linda Lilburne

Comparison of Three Spatial Sensitivity Analysis Techniques


Nathalie Saint-Geours

Christian Lavergne

Jean-Stéphane Bailly

Frédéric Grelot

Analyse de sensibilit é globale d'un mod èle d' évaluation économique du risque d'inondation

Keywords: analyse de sensibilité globale, indices de Sobol, spatial, inondation global sensitivity analysis, variance-based, Sobol' indices, spatial, flood

1 : L'analyse de sensibilité globale peine à se développer dans le champ de la modélisation environnementale. Dans sa formulation initiale, elle est limitée à l'étude de modèles Y = f (X 1 , . . . , X p ) où les variables d'entrée X j et la sortie Y sont scalaires, alors que nombre de modèles environnementaux incluent une dimension spatiale marquée, soit qu'ils fassent appel à des cartes comme variables d'entrée, soit que leurs sorties soient distribuées spatialement. Au travers d'une étude de cas détaillée, nous présentons dans cet article une extension de l'analyse de sensibilité globale à l'étude de modèles spatialisés. Le modèle étudié, nommé ACB-DE, est un outil d'évaluation économique du risque d'inondation. La réalisation de cartes d'indices de sensibilité permet d'étudier les sorties spatialisées du modèle ACB-DE et de rendre compte de la variabilité spatiale des indices de Sobol. L'influence relative des variables d'entrée à différentes échelles d'étude est analysée par la réalisation de cartes d'indices de sensibilité de résolution croissante. L'analyse réalisée permet d'identifier les variables d'entrée incertaines qui expliquent la plus grande part de la variabilité de l'indicateur économique fourni par le modèle ACB-DE ; elle apporte un éclairage nouveau sur le choix de l'échelle adéquate de représentation spatialisée de cet indicateur selon la précision des variables d'entrée.

Introduction

L'analyse de sensibilité d'un modèle Y = f (X 1 , . . . , X p ) vise à quantifier l'impact de l'incertitude pesant sur les p variables d'entrée X i du modèle sur la variabilité de sa sortie Y . L'analyse de sensibilité globale basée sur la variance proposée par [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF] peut être utilisée pour étudier des codes numériques de type boîte noire peu coûteux en temps de calcul. Elle s'appuie sur la décomposition de la variance de Y en variances conditionnelles ; on y définit des indices de sensibilité S i (parfois appelés indices de Sobol), compris entre 0 et 1, qui traduisent la part de chaque variable d'entrée X i dans la variance de Y :

S i = Var [E(Y | X i )]
Var(Y ) Dans sa formulation initiale, l'analyse de sensibilité globale basée sur la variance est limitée à l'étude de modèles où les variables d'entrée X i sont scalaires, tout comme la sortie Y . Or, dans le vaste champ de la recherche en environnement, notamment dans l'étude des risques naturels, nombre de modèles incluent une dimension spatiale marquée, soit qu'ils fassent appel à des cartes comme variables d'entrée, soit que leurs sorties soient spatialisées (cartes de risque par exemple). De ce fait, l'analyse de sensibilité peine à se développer dans ces champs thématiques, comme l'a montré [START_REF] Delgado | Sensitivity analysis in multicriteria spatial decisionmaking: a review[END_REF]. Nous présentons ici une application de l'analyse de sensibilité globale basée sur la variance à un modèle d'évaluation économique du risque d'inondation (modèle ACB-DE), dont les variables d'entrée et les sorties sont distribuées spatialement.

Modèle ACB-DE

Le modèle ACB-DE est un outil d'évaluation économique du risque d'inondation développé par [START_REF] Erdlenbruch | Une analyse coûtbénéfice spatialisée de la protection contre les inondations : Application de la méthode des dommages évités à la basse vallée de l'Orb[END_REF]. Il vise à caractériser l'exposition d'un territoire au risque d'inondation, au moyen d'un indicateur spatialisé décrivant les dommages liés aux crues potentielles. Les données utilisées dans ce travail sont issues d'une étude de cas sur la basse vallée de l'Orb, dans l'Hérault (superficie de 100 km 2 ). Le modèle ACB-DE fait appel à des modules de natures différentes (hydrologique, hydraulique, économique), et combine six variables d'entrée notées X 0 à X 5 : X 2 , X 3 et X 5 appartiennent respectivement à R 5 , R 20 et R. Les trois autres variables d'entrée sont distribuées spatialement et sont représentées sous forme de couches dans un logiciel de Système d'Information Géographique : X 0 est une carte décrivant l'étendue des plaines d'inondation, X 1 est un Modèle Numérique de Terrain et X 2 une carte d'occupation du sol. La résolution horizontale de ces cartes est de 5 m. La sortie Y du modèle ACB-DE est une carte donnant en tout point u du territoire les dommages évités Y (u) en ce point. Sa résolution horizontale est de 5 m. A partir de cette sortie Y , on crée quatre autres cartes plus grossières notées Y (200) , Y (400) , Y (800) et Y (1600) , de résolution respectives s = 200, 400, 800 et 1600 m, à destination des gestionnaires de bassins versants (Figure 1). Plus précisement, pour une résolution s donnée, Y (s) désigne la carte donnant en toute cellule C de côté s la somme Y C des dommages évités sur cette cellule :

Y C = u∈C Y (u) du.

Analyse de sensibilité globale spatialisée

Plusieurs travaux ont déjà apporté des éléments de réponse aux problèmes posés par l'analyse de sensibilité globale d'un modèle spatialisé ; ils se sont intéressés à l'intégration dans l'analyse de variables d'entrée spatialisées et à l'étude de sorties spatialisées. On renvoie à Lilburne et Tarantola (2009) pour une liste de références sur ce sujet.

On propose ici de combiner deux approches pour mener à bien l'analyse de sensibilité du modèle ACB-DE : on associe le traitement des variables d'entrée spatialisées selon la méthode de Lilburne et Tarantola (2009) et l'analyse de sorties spatialisées via des cartes de sensibilité comme proposé par Marrel et al (2009). On explore de plus un point non abordé à notre connaissance, qui est l'impact sur la valeur des indices de Sobol de la résolution s à laquelle est représentée la variable de sortie Y spatialisée.

Description des incertitudes Dans un premier temps, on spécifie l'incertitude pesant sur chaque variable d'entrée (Table 1). Des distributions de probabilité sont identifiées pour les variables X 2 , X 3 et X 5 . Afin de rendre compte de l'incertitude sur les deux variables spatialisées (X 1 , carte d'occupation du sol et X 4 , Modèle Numérique de Terrain), deux jeux de n 1 = 1000 et n 4 = 100 réalisations aléatoires de ces variables ont été générés : un champ d'erreur Gaussien présentant une structure spatiale de covariance modélisée par un variogramme est ajouté au Modèle Numérique de Terrain initial selon un algorithme de simulation Gaussienne séquentielle proposé par [START_REF] Castrignano | Accuracy assessment of Digital Elevation Model using stochastic simulation[END_REF] ; des confusions probabilistes entre classes sont introduites dans la carte d'occupation du sol. Estimation des indices de sensibilité Dans un second temps, un échantillon de taille N = 28 672 est constitué afin d'explorer l'espace des variables incertaines X 1 à X 5 , puis le modèle ACB-DE est simulé en chaque point de cet échantillon. La variable spatialisée X 1 (resp. X 4 ) est intégrée dans cet échantillon en associant un entier compris entre 1 et n 1 (resp. n 4 ) à chaque réalisation aléatoire préalablement générée de cette variable, ces réalisations étant considérées équiprobables. On estime alors les indices de sensibilité de premier ordre de chacune des variables X j par rapport à Y , notés S j (Y ), selon les expressions proposées par [START_REF] Saltelli | Global Sensitivity Analysis -The Primer[END_REF]. Des intervalles de confiance à 90 % sont estimés pour chacun des indices de sensibilité par bootstrap (nombre d'échantillons bootstrap n boot = 100).

Réalisation de cartes d'indices de sensibilité

Pour chacune des quatre sorties spatialisées Y (200) à Y (1600) , on établit des cartes de sensibilité pour chaque variable d'entrée X 1 à X 5 . Ces cartes sont définies sur la même grille régulière que la sortie Y (s) considérée. La valeur de la carte de sensibilité de la variable X j en une cellule C est égale à l'indice de sensibilité S j (Y C ) de cette variable par rapport à la somme Y C des dommages évités sur cette cellule. On obtient ainsi un jeu de 20 cartes de sensibilité qui traduisent l'influence de l'incertitude des différentes variables d'entrée X 1 à X 5 sur la variabilité de la carte Y des dommages évités à différentes résolutions spatiales.

Résultats

La Figure 2 présente les cartes d'indices de sensibilité de trois variables d'entrée (X 1 , X 2 et X 3 ) par rapport à la carte Y des dommages évités à deux résolutions différentes, s = 200 et s = 1600 m. Sur chacune de ces cartes, on observe une variabilité spatiale des indices de sensibilité. À une résolution s donnée, l'influence d'une variable incertaine X i sur la variance de la sortie Y C (somme des dommaes évités sur la cellule C) varie d'une zone à l'autre du territoire. On peut identifier deux éléments de structuration spatiale majeurs : les effets de bord et la présence de villes (bâti urbain).

Par ailleurs, on observe que les cartes de sensibilité par rapport à la sortie Y (200) sont différentes des cartes réalisées par rapport à la sortie Y (1600) . Pour permettre une comparaison entre les cartes d'indices de sensibilité aux différentes échelles, on a calculé pour chaque variable d'entrée X j et chaque résolution s = 200, 400, 800 et 1600 m la valeur moyenne µ (s) j de l'indice de sensibilité S j (Y C ) sur l'ensemble des cellules C de la grille.

On montre que pour chaque variable d'entrée X j , l'indice moyen µ (s) j dépend de manière monotone de la résolution s. L'indice de sensibilité des variables distribuées spatialement (ici X 1 et X 4 ) diminue à mesure qu'augmente le côté s des cellules sur lesquelles sont calculés les dommages évités. Un effet de compensation des incertitudes en est responsable : les incertitudes qui affectent les variables spatialisées sont locales ; lorsqu'une large zone est considérée, ces incertitudes locales se compensent et leur impact relatif par rapport aux autres variables diminue. Symétriquement, l'indice de sensibilité des variables non distribuées spatialement (ici X 2 , X 3 et X 5 ) augmente avec le côté s de la cellule C sur laquelle on calcule les dommages évités. Selon la résolution s à laquelle on représente le champ Y des dommages évités, les variables d'entrée qui expliquent sa variabilité ne sont donc pas les mêmes.

Conclusion

L'approche proposée est applicable à un large éventail de modèles spatialisés, à la condition qu'ils soient peu coûteux en temps de calcul. Elle offre de nouvelles perspectives dans la compréhension des liens entre effets d'échelle, incertitudes et sensibilité, questions qui sont propres à ce type de modèles. Des travaux complémentaires sont à mener pour préciser les effets de biais induit par la méthode de représentation de l'incertitude sur les variables spatialisées, pour optimiser l'échantillonnage de réalisations aléatoires de ces variables, ou encore pour mieux comprendre et prédire la structure spatiale des cartes de sensibilité.

Figure 1 :

 1 Figure 1: Y (200) , carte des Dommages Évités Moyens Annualisés sur une grille régulière constituée de cellules de côté s = 200 m. En chaque cellule C est représentée la somme Y C des dommages évités sur cette cellule.

Figure 2 :

 2 Figure2: Cartes de sensibilité. Pour chacune des trois variables X 1 à X 3 , la carte du haut (resp. du bas) présente la carte de sensibilité de la variable X j par rapport à la sortie spatialisée Y(200) (resp. Y (1600) ) : sur chaque cellule C de côté s = 200 m (resp. 1600 m), on représente S j (Y C ), indice de sensibilité de la variable X j par rapport à la somme Y C des dommages évités sur la cellule.

Un article présentant ce travail plus en détails a été soumis au Journal de la SFdS pour un numéro spécial Analyse de sensibilité.