Integrated hydrological modeling of domains with complex subsurface features requires many highly uncertain parameters. Performing a global uncertainty analysis using an ensemble of model runs can help bring clarity as to which of these parameters really influence system behavior and for which high parameter uncertainty does not result in similarly high uncertainty of model predictions. However, already creating a sufficiently large ensemble of model simulation for the global sensitivity analysis can be challenging, as many combinations of model parameters can lead to unrealistic model behavior. In this work we use the method of active subspaces to perform a global sensitivity analysis. While building up the ensemble, we use the already-existing ensemble members to construct low-order meta-models based on the first two active-subspace dimensions. The meta-models are used to pre-determine whether a random parameter combination in the stochastic sampling is likely to result in unrealistic behavior so that such a parameter combination is excluded without running the computationally expensive full model. An important reason for choosing the active-subspace method is that both the activity score of the global sensitivity analysis and the meta-models can easily be understood and visualized. We test the approach on a subsurface-flow model including uncertain hydraulic parameters, uncertain boundary conditions and uncertain geological structure. We show that sufficiently detailed active subspaces exist for most observations of interest. The pre-selection by the meta-model significantly reduces the number of full-model runs that must be rejected due to unrealistic behavior. An essential but difficult part in active-subspace sampling using complex models is approximating the gradient of the simulated observation with respect to all parameters. We show that this can effectively and meaningfully be done with second-order polynomials.

Water flow in the subsurface is an integral part of the water cycle. In recent years, integrated hydrological modeling based on the partial differential equation (pde), coupling flow in the subsurface and on the land surface, has become a rather standard tool

A global sensitivity approach that does not directly fit into any of the categories listed above but has recently gained increased attention is the active-subspace method

A general problem when performing any type of global sensitivity analysis is the choice of how to sample the parameters. Apart from defining ranges and distributions of single parameters, which are unique to the problem at hand and can often be addressed by experts in the field, questions related to unfavorable parameter combinations are harder to deal with a priori. Unfavorable parameter combinations may lead to model behavior that is not observed in reality, such as severe floods or strong droughts. In regionalized (also known as generalized) sensitivity analysis, such parameter sets are classified as non-behavioral, and statistical differences between the behavior and non-behavioral parameter sets are sought

Recognizing that a sufficiently large set of model runs is needed for a reliable stochastic analysis,

In this paper, we use a meta-model, derived by the active-subspace method, to pre-determine presumably behavioral parameter sets and perform the global sensitivity analysis with the full model using the pre-selected parameter sets. By this we aim at reducing the number of discarded simulations using the full model. We use two-dimensional active subspaces to derive both multiple meta-models and sensitivity patterns for an integrated surface–subsurface-flow model. The rest of the paper is structured as follows: Sect.

Flow in the subsurface is computed using the software HydroGeoSphere

In this section we consider a general function

The basic idea of the active-subspace decomposition is to find primary directions in the original parameter space, composed of linear combinations of the parameters, along which the solution

An active subspace of

For assessing the global sensitivity of each parameter in

A major issue with computing an active subspace for a subsurface-flow model is that it requires the gradient of the target quantity

Our standard fit is the second-order polynomial with cross terms. If the set of model runs is smaller than about twice the number of required coefficients, we use the gradient fit 2, excluding the second-order cross terms. The linear fit 3 implies that the gradient

In theory, also higher-order polynomials could be used to approximate the local gradient vectors. In practice, however, we are limited by the number of regression coefficients we need to estimate. With the 32 parameters considered in this work, we require a rather large ensemble of full-model runs to obtain the 561

With a functional active-subspace decomposition, we may construct a low-order approximation of the observation (

To this end, we construct the vector

In summary, the construction of an active subspace contains two strong approximations which both give rise to errors: (1) the active-subspace decomposition itself (dimension reduction) and (2) the gradient approximation. As these two errors can be strongly correlated, it is difficult to show the effect of the dimension reduction when the gradient approximation is still uncertain. However, in Sect. 4.3 we attempt to show the effect on the total error when altering the accuracy of the gradient approximation.

A key difficulty in running complex models with random parameters drawn from wide prior distributions is that a significant number of the resulting model simulations may show behavior that is contradictory to the prior knowledge of the modeled system. Such non-behavioral runs should be discarded in subsequent analyses, which implies that running them was a waste of computational resources.

An approach to limit the number of non-behavioral model runs could be adaptive sampling, in which a meta-model (i.e., a simplified, fast-running low-order approximation of the true model) is used first to predict whether a parameter-set is behavioral. In this study we utilize the ability of an active subspace to construct low-order meta-models between our unknown parameters, represented as active variables, and the chosen observations whose behavior we wish to control. In our application, we use the third-order polynomials of the active subspaces, explicated in Eq. (

The setup of an adaptive sampler using active subspaces consists of the following seven steps.

Run a first set of flow models with random parameters drawn from a wide distribution of plausible values. In this work we use 500 as our initial sample size and apply Latin hypercube sampling.

For every unique observation type related to a behavioral target, construct a sufficiently detailed active subspace. If the gradient computation allows it, several subspace dimensions would be better than one. Here we use two active-subspace dimensions (i.e., two active variables), and, hence, our meta-model is a surface in the two-dimensional space of active variables. For each meta-model, an NSE value larger than 0.7 is required to be considered for pre-assessing the behavior of new parameter sets in the following steps.

A new candidate of the full parameter vector is now drawn from the same initial distribution as that used in step 1. This parameter vector is projected onto the active subspace(s), and the meta-model is used to approximate the behavior of all target predictions.

The new candidate will be accepted at stage one if any of the following criteria are met:

All approximated behavioral targets are on the permitted side of their limits.

While a target is on the non-behavioral side of the limit, it is within a reasonable distance of the limits. In practice, this is implemented as a linear decay function from 1 at the limit to 0 at an outer (user specified) point. If the decay function value is larger than a random number drawn from a uniform distribution, the candidate is accepted at stage one. This criterion is implemented to construct a soft region around the limits that accounts for the imperfection of the meta-model. The outlined approach is similar to the classical Metropolis–Hastings sampling

With a 10 % probability, a candidate is accepted at stage one independent of its predicted performance. This criterion is included to make sure that we maintain a sufficiently good sample of the full parameter space so that the recalculation of the active subspace (see below) still sees the unwanted regions.

For each stage-one-accepted candidate, we run the full flow model to obtain the prediction of the real model (stage two).

After performing a predefined number (in our case 100) of flow simulations, recalculate the active subspaces using all flow-simulation outputs obtained so far.

Repeat steps two through six until the sample size is large enough for the purpose of the stochastic modeling. This is a model-purpose-specific choice and can be done both on hard limits to the (stage-two) simulated data or on the number of flow-model runs. Here, we require 10 000 runs of the flow model (i.e., 9500 stage-one acceptances plus 500 initial samples, which are per se stage-one accepted).

It is important to note that all post-sampling sensitivity analyses performed in this work are done on the subset of the sampled parameter sets that are deemed behavioral after running the full HydroGeoSphere flow model (stage-two accepted). Hence, we use the meta-model only as a pre-selection tool to avoid sampling those regions of the parameter space that will clearly generate non-behavioral runs. As we aim to sample 10 000 parameter sets (which are stage-one accepted), the analyses will be performed on a notably smaller number of parameter sets.

In this work, we have chosen to construct the meta-model using two active variables. Although more active variables could potentially lead to higher accuracy of the meta-model, we saw no major improvement when increasing the number of subspace dimensions beyond two. Along the same line of thought, model outcomes in two subspace dimensions can easily be visualized, thus facilitating an intuitive judgement of the goodness of the meta-model. In this light, we refrain from going beyond two active-subspace dimensions in the current work. Other application may require considering more active dimensions.

It should be noted that applying the active-subspace sampling does not necessarily restrict the sensitivity analysis to calculating the activity score (Eq.

The test bed used in this paper is a steady-state flow model setup and run in HydroGeoSphere. It draws its main features from the catchment of the stream Käsbach in the Ammer valley in southwestern Germany

As illustrated in Fig.

Illustration of the modeled catchment, including important features and being surrounded by an explicit view on the different geological features.

The model domain measures about 4 km

The boundary conditions, set up to allow water to leave the domain both through the surface and the subsurface, are as follows. The bottom of the model features a Dirichlet boundary with values read in from a larger-scale model of the region

To avoid long runtimes and complications of complex topsoils (including plant–atmosphere interactions), which are unimportant once the steady state is reached, the top of the HydroGeoSphere model is 1 m below the land surface. Flow across the top boundary is only incoming and modeled as a Neumann boundary, corresponding to the steady-state groundwater recharge. The recharge varies with land use, split into three categories: cropland, grassland and forest, in which urban areas are treated as grassland. It should be noted that starting the model at 1 m below surface still allows for a notable unsaturated zone to develop in the model domain; only the uppermost meter is missing. Also, the outgoing drain fluxes described above are applied to this top boundary that is 1 m below the surface, hence requiring that the exit pressure head is 1 m larger than the ponding pressure head described above. Technically, this does not apply to stream nodes, which are considered to be the real top of the porous medium (that is, there is no unsaturated zone on top of a stream).

For the evaluation of the model sensitivities, we consider four observation quantities: (1) the discharge almost at the outlet of the catchment (gauge C in Fig.

We computed the geological-unit-specific residence time using the visualization software TecPlot360, considering 199 particles that each start in an observation well, 5 m below the groundwater table. We separated the discrete particle tracks into the segment spent in each geological layer and computed the total residence time per layer as the mean over all 199 particles. It should be noted that this way of computing travel times is slightly imprecise, as the velocity-field output of HydroGeoSphere is non-conforming and Tecplot particle tracking is not primarily designed for quantitative outputs. However, the purpose of the residence-time calculation is not an exact prediction of exposure times in the real Käsbach catchment but rather a qualitative transport indicator. As we used the same method with the same numerical parameters on the same grid in all model runs of our stochastic ensemble, we believe that the associated variability in computed residence times is good enough for the inter-comparison between different stochastic runs.

In the setup used in this paper, 32 parameters related directly to the flow model are randomized. All parameters are sampled from a uniform distribution. Table

Sampling ranges for of stochastic parameters considered.

mo: Upper Muschelkalk. ku: Lettenkeuper. km1: Gipskeuper. km1-w: weathered Gipskeuper. Q: Quaternary fillings.

Six examples of distinctively different realizations of the geological layer Lettenkeuper. Color shows the natural logarithm of the saturated hydraulic conductivity.

Besides these material properties of the lithostratigraphic units and boundary conditions, also the exact subsurface structure is uncertain. To address this uncertainty, we drew three parameters controlling the size of the main geological layers from uniform distributions: the vertical offset of the fault running north–south through the domain (see Fig.

The stochastic engine for HydroGeoSphere is set up in and controlled by MATLAB. The full stochastic suite is run on a midrange cluster with 20 nodes, each featuring two Intel Xeon L5530 eight-core 2.4 GHz processors with 72 GB RAM. The 10 000 samples discussed later took on this setup about 96 000 CPU hours (

In the present work, we define five behavioral targets that are all based on expert knowledge about the modeled catchment. In the following, we specify the behavioral target values and the point of maximum deviation from that target for which we may probabilistically accept a simulation run, denoted the “outer point”:

Flooding, here viewed as water leaving the domain through the top drain at the surface at places outside of the streams, occurs in the model at a hydraulic head of 0.2 m above the land surface. The total flooding in the domain should not exceed

At measurement gauge C (Fig.

Knowing that stream B produces flow, a minimum flow is set to

The stream residing on the steep eastern side of the hill is known to only produce flow under extreme conditions. Hence, at the steady state the flow in stream A should be minimal. The maximum accepted flow is therefore set to

It is known that the catchment in question loses a notable amount of its water to the subsurface. A rough estimate is that, in the real catchment, the stream flow amounts to

In a preliminary test of a model very similar to the one used here, and with the same randomized parameters, we performed Monte Carlo simulations of the full model without pre-selecting presumably behavioral parameter sets. Here about 75 % of a total of 10 000 runs had to be discarded only due to severe flooding. This highlights that it is highly beneficial if the sampling is targeted only to simulations that show a response that is, within reason, representative of the modeled domain. In the case of flooding in the domain, it is not a single parameter that controls this behavior but a complex relation between many parameters, deeming a priori decisions about behavioral parameter ranges unfeasible. As all targets are based on knowledge about the real catchment on which the current model is based, all model outputs produced by a stage-two-accepted model will be in line with what we would expect to be realistic in the catchment. However, it is important to note that even though it would be possible to point out a location in the active subspace which corresponds to a real observation, the active variables themselves are non-unique with respect to the flow parameters so that a simple back transformation from the active subspace to flow-model parameter space is not possible.

To allow the reader to better see the 3-D structure in the results presented here, the main results can be viewed in a plug-and-play app designed for MATLAB, denoted the “Active Subspace Pilot”, which is available as supplementary information to this publication.

The effect of using the active subspaces as a sampling strategy for the flow simulations can be seen in Fig.

Marginal posterior distributions of the parameters influenced by the adaptive sampling. Blue bars show the sampled posterior, and red bars show the constrained posterior sample used in the sensitivity analysis.

Two-dimensional active subspaces for three of the behavioral targets: top drain (limited flooding – target 1), ratio (target 5) and flow in stream A (target 4). The observed values are illustrated by the color.

The blue bars of Fig.

The red bars in Fig.

An example of the performance of the adaptive sampler for the top drain observation (limited flooding target).

Figure

Figure

In this section, we use the active-subspace method to analyze parameter sensitivities. In this analysis we only consider the behavioral parameter sets. That is, the red bars in Fig.

The aim of a sensitivity analysis is to identify how influential the individual parameters of a model are on (a set of) observations. In global sensitivity analysis, this is evaluated over the entire parameter space. We start the discussion with the least influential parameters in our application. None of the observations considered depend on the two van Genuchten parameters

We now consider the significant parameters. Figure

Active subspaces and square root of activity scores for the observations “stream discharge at gauge C” and “net flow across subsurface boundaries”. Please note that the plot is limited to the seven most important parameters.

Square root of activity score for 199 wells

Figure

Two-dimensional active subspaces

As a last observation, we consider the total residence time in the geological layers of the aquifer, which may be a relevant proxy for applications to reactive transport. Figure

In contrast to previous works with active subspaces in subsurface flow, we approximate the gradient

To test the consistency of the results, we drew 2500 samples in 1000 repetitions using classic bootstrap sampling without replacement from the original sample (4533 members). In each repetition, we computed the activity scores and NSE of the meta-models based on the three gradient approximations. Figure

Square root of activity score and corresponding NSE for flow at gauge C using different approximations to compute the gradients. Each approximation is used to compute 1000 active subspaces based on a bootstrap resampling with 2500 samples from the original 4533. Only the seven most important parameters are shown in the activity-score plot.

To evaluate the goodness of each approximation, Fig.

In this work we have applied the method of active subspaces to an integrated hydrological model of a small catchment with focus on subsurface flow. We used active subspaces to construct meta-models with two active subspaces rather than 32 uncertain parameters. The meta-model was used to constrain the stochastic sampling of the parameter space to five behavioral conditions. The active subspaces of the accepted full-model runs were used to compute the global sensitivity of four modeled observations to the parameters. The sensitivity analysis showed that not only hydraulic-conductivity values of the major layers but also their physical extent are important. However, depending on the location and type of observations, different sensitivities were found. This highlights the well-known fact that multiple dissimilar observations are needed to constrain uncertain variables of a catchment model. In the adaptive sampling we learned that certain combinations of unfavorable parameter values were clearly avoided. Most of the non-behavioral parameter combinations were not obvious beforehand but could be identified by applying the meta-model, which significantly improved the sampling efficiency.

The choice of meta-model used in this work (third-order polynomial of the first two active-subspace dimensions) was somewhat arbitrary. The number of different meta-models applied to hydrological problems is large

Choosing a low-order polynomial as a meta-model implies smoothness, and, hence, the meta-model does not exactly fit all model runs.

Even though the efficiency of the active-subspace sampler is higher than a rejection sampler without pre-selection, the rejection rate of

Overall, we draw the following conclusions from this work:

The method of active subspaces can be applied with little effort and good results to complex subsurface flow and transport problems. This holds not only when subsurface properties are uncertain but also for the geometries of geological units and boundary conditions.

The two-stage rejection sampling using a meta-model based on the active subspaces can drastically decrease the number of simulations needed to obtain a certain number of behavioral simulations. An additional positive aspect for the application by practitioners is the ease of visualization and intuitive understanding when using a one- or two-dimensional active subspace.

Using a quadratic rather than linear fit to estimate the gradients in the construction of the active subspace resulted in a much improved subspace decomposition. It is also a prerequisite to construct more than one subspace dimension.

The Supplement includes the MATLAB R2019a code Active Subspace Pilot to visualize the active subspaces and scores for all model runs. This code is written as a MATLAB app and also includes the data.

This Appendix contains additional plots showing the performance of the different gradient approximations for the four behavioral targets not presented in the main article. Each approximation is used to compute 1000 active subspaces based on a bootstrap resampling with 2500 samples from the original 4533 ensemble members. This holds for all observations apart from those for the flow across the top boundary, where the full ensemble is used to draw samples from. Only the seven most important parameters are shown in each activity-score plot.

Square root of activity score and corresponding NSE for flooding flow across the top boundary using different approximations to compute the gradients. In this case, the bootstrap sample is drawn from the full sample population of 10 000.

Square root of activity score and corresponding NSE for flow in stream A (see Fig.

Square root of activity score and corresponding NSE for unwanted flow in stream B (see Fig.

Square root of activity score and corresponding NSE for the ratio of stream discharge to incoming recharge using different approximations to compute the gradients.

The supplement related to this article is available online at:

Simulations and code development were performed by DE. Both authors contributed to developing and writing the paper. OAC was responsible for acquisition of the funding.

The authors declare that they have no conflict of interest.

This work was supported by the Collaborative Research Center 1253 CAMPOS (Project 7: Stochastic Modeling Framework of Catchment-Scale Reactive Transport), funded by the German Research Foundation (DFG; grant agreement SFB 1253/1).

This research has been supported by the German Research Foundation (grant no. SFB 1253/1).This open-access publication was funded by the University of Tübingen.

This paper was edited by Alberto Guadagnini and reviewed by two anonymous referees.