Global sensitivity analysis using multi-resolution polynomial chaos expansion for coupled Stokes–Darcy flow problems

Determination of relevant model parameters is crucial for accurate mathematical modelling and efficient numerical simulation of a wide spectrum of applications in geosciences. The conventional method of choice is the global sensitivity analysis (GSA). Unfortunately, at least the classical Monte-Carlo based GSA requires a high number of model runs. Response surfaces based techniques, e.g. arbitrary Polynomial Chaos (aPC) expansion, can reduce computational effort, however, they suffer from the Gibbs phenomena and high hardware requirements for higher accuracy. We introduce GSA for arbitrary Multi-Resolution Polynomial Chaos (aMR-PC) which is a localized aPC based data-driven polynomial discretization. The aMR-PC allows to reduce the Gibbs phenomena by construction and to achieve higher accuracy by means of localization also for lower polynomial degrees. We apply these techniques to perform the sensitivity analysis for the Stokes–Darcy problem which describes fluid flow in coupled free-flow and porous-medium systems. We consider the Stokes equations in the free-flow region, Darcy’s law in the porous-medium domain and the classical interface conditions across the fluid–porous interface including the conservation of mass, the balance of normal forces and the Beavers–Joseph condition for the tangential velocity. This coupled problem formulation contains four uncertain parameters: the exact location of the interface, the permeability, the Beavers–Joseph slip coefficient and the uncertainty in the boundary conditions. We carry out the sensitivity analysis of the coupled model with respect to these parameters using the Sobol indices on the aMR-PC expansion and conduct the corresponding numerical simulations.


Introduction
Fluid flows through a coupled system consisting of a freeflow region and an adjacent porous-medium domain appear in many technical applications and environmental settings such as industrial filtration, wood drying, water-gas management in fuel cells, interaction of surface and subsurface flows [7,13,16,26,48].Solving such coupled flow problems at the microscale is computationally demanding, if not impossible, for applications as the detailed pore-scale geometry should be resolved.Therefore, macroscale model formulations containing two different flow models in the free-flow and porous-medium domains together with an appropriate set of coupling conditions on the fluid-porous interface are considered in practice [4,17,42,71].Such coupled flow models contain several effective coefficients whose values are typically uncertain.The quantification and correct choice of the relevant model parameters are important for physically consistent modelling and accurate numerical simulations of flow and transport processes in coupled systems for various applications in geosciences.
In this paper, we consider the most commonly used coupled flow model containing the Stokes equations in the free-flow domain and Darcy's law in the porous-medium region.The classical set of coupling conditions including the conservation of mass across the interface, the balance of normal forces and the Beavers-Joseph condition on the tangential velocity component is applied [8,18,41,42,50].Alternative coupling concepts for such flow problems exist, e.g.[4,21,32,40,49,51].Some of them include unknown parameters which still have to be determined in order to make these models computable [4,51,66].Also in case of the classical coupling conditions considered in this work, several effective parameters conceptualising the macroscale model are typically uncertain, even if the detailed pore-scale geometry is available.We consider three uncertain parameters in the model of interest in this work: the exact position of the fluidporous interface, the permeability of the porous medium and the Beavers-Joseph slip coefficient.Additionally, we introduce an uncertainty in the outflow boundary conditions to be able to model coupled problems with different flow directions to the interface.Several computational methods were developed in the past to address some of these uncertainties, e.g., [3,37].
There are several attempts in the literature to specify the uncertain parameters in the Stokes-Darcy model.However, to the best of the authors knowledge, there is no such systematic study for general flow problems.For periodic porous media, permeability can be computed in the classical way using the theory of homogenisation and solving additional cell problems [28].For arbitrary porous media, numerical upscaling techniques with different pore-scale resolved models can be applied in order to determine the intrinsic permeability, e.g.[58,67] and references therein.In this case, there are uncertainties in the choice of boundary conditions for the local flow problems in the representative elementary volume.Even less information exists on where to locate the sharp fluid-porous interface so that the pore-scale resolved solutions fit to the macroscale solutions.The interface is located on the top of the first row of solid inclusions in [8,39], other possibilities are considered in [14,57].The most challenging problem is the determination of the Beavers-Joseph slip coefficient.One can find the values of the slip coefficient for four materials in the original work [8].Several values for periodic isotropic porous media are provided in [46,57,70], but correct determination of this parameter still stays an open question for non-periodic and anisotropic porous structures.
The current paper aims to investigate sensitivity of the coupled Stokes-Darcy model to the exact position of the fluidporous interface, porous-medium morphology reflected in the intrinsic permeability, the Beavers-Joseph slip coefficient and the outflow boundary conditions.The knowledge about parameter relevance is necessary for model reduction, model calibration and parameter inference.
Importance of the model parameters and their influence on the overall modelling procedure is not trivial to understand even for the experienced modeller.On the one hand, uncertain parameters could contain unseen uncertainty, on the other hand, they could produce a joint effect that may become significant during the modelling procedure.Moreover, any manual identification of relevant parameters and their combinations is usually limited by the available resources.Sensitivity analysis provides identification of the parameters with the highest impact of uncertainty that could serve as the basis for model reduction.While the local sensitivity analysis studies the variation of the model response close to some selected basis point in the parameter space [63], the global sensitivity analysis (GSA) considers the averaged variation of the model response over the whole domain [27,61,63].In the present work, we use the Sobol and total sensitivity indices [27,60,61] to perform the GSA.
The GSA is commonly used for many applications in geosciences, e.g.[23,55].Nevertheless, the classical Monte-Carlo (MC) averaging based GSA [61,63] often requires high number of model runs, which becomes challenging with increasing computational costs for a single model run.Instead, the Polynomial Chaos Expansion (PCE) based surrogate models allow to avoid additional model evaluations after building the surrogate model [23,62].
Research on orthogonal polynomials goes back at least into the 19th century.The Ansatz of the homogeneous chaos was introduced by Wiener [60] and extended by Cameron and Martin [12].Since Ghanem and Spanos applied the PCE in context of partial differential equations [60], many researches have used and extended the concept of PCE in several uncertainty related areas [6,19,22,34,38,45,53,59,62,69].Of particular interest for this work is the data-driven Ansatz of the arbitrary Polynomial Chaos (aPC) expansion, which was introduced in [52].
The idea of the multi-resolution framework goes back to the early 90s [1,2,9].In the area of Uncertainty Quantification (UQ) multi-resolution framework was introduced in [43,68] and also widely used in research, in particular, in context of nonlinear model response [10,11,25,30,34,35,44,54,65].One of the drawbacks of the multi-resolution framework is the rapidly increasing number of stochastic subdomains for increasing number of uncertain parameters.Adaptivity strategies that reduce the number of stochastic subdomains w.r.t. the data are proposed, e.g., in [36,65].
In this work, we introduce GSA for arbitrary Multi-Resolution Polynomial Chaos (aMR-PC) [36], which combines the flexibility of the data-driven concept of the aPC [52] and the reliability of the multi-resolution/ multi-element based PCEs [11,19,34,35,43,65,68].More precisely, construction and performance of aMR-PC and multi-wavelet based adaptivity for aMR-PC were discussed in [36].In present work, we introduce the piecewise-polynomial-coefficients-based GSA for aMR-PC that allows to perform GSA without additional model-or surrogate model evaluations after building the aMR-PC-surrogate model such that the concept of GSA for plain PCE introduced in [62] becomes available on the multi-resolution framework.For the sake of comparison, we provide the application of the introduced GSA method on the Ishigami function, which is one of the standard benchmarks in context of GSA [62].We apply the proposed GSA method for the sensitivity analysis of the coupled Stokes-Darcy flow model, which contains several uncertain parameters, and lay the basis for future research on parameter identification and model reduction for coupled flow problems.
The paper is organised as follows.In Section 2, the coupled macroscale Stokes-Darcy problem formulation is described and the typical values and intervals of the uncertain model parameters are introduced.In Section 3, we provide a brief overview of the aMR-PC expansion, which is used as the surrogate model.Section 4 is devoted to the global sensitivity analysis employing the Sobol indices.Numerical simulation results on the sensitivity analysis for the Stokes-Darcy model are presented in Section 5. Finally, the discussion follows in Section 6.In Appendix, the Sobol indices for the aPC and aMR-PC expansions of the Ishigami function are provided.

Geometrical setting and assumptions
We consider steady-state non-compositional single-fluidphase flows at constant temperature both in the freeflow domain and in the porous-medium region.In noncompositional flow systems only one component (chemical species) is present in the fluid phase, therefore, flow equations for the phase as a whole are considered without modelling any transport of individual chemical components.The flow is assumed to be creeping at low Reynolds numbers.The fluid is incompressible and the viscosity is constant.The porous medium is considered to be homogeneous, isotropic and having constant porosity.The solid phase is rigid.
The coupled flow system contains the free-flow domain Ω ff ⊂ R 2 and the porous-medium region Ω pm ⊂ R 2 , which are separated by a sharp fluid-porous interface Γ = Ω ff ∩ Ω pm (Fig. 1).The interface is considered to be flat and simple, i.e. no thermodynamic properties can be stored in or transported along Γ .The precise location of the fluid-porous interface is uncertain.It can be placed on the top of the first row of solid inclusions, slightly above, between the first and the second rows of solid grains, or at a distance of maximum 5 pore diameters from the top of solid inclusions [14].

Flow problem formulation
Under assumptions given in Section 2.1, fluid flow in the free-flow domain can be described by the Stokes equations where v ff is the fluid velocity, p ff is the fluid pressure, T is the rate of strain tensor, I is the identity tensor, ρ is the fluid density and g is the gravitational acceleration.
The following boundary conditions on the external boundary Γ ff = ∂Ω ff \ Γ are considered where n ff is the unit outward normal vector from the domain Ω ff on its boundary, Γ ff = Γ D,ff ∪ Γ N ,ff and v ff is given.In this work, we consider the outflow boundary Γ N ,ff with the uncertain opening width h (Fig. 1) in order to model different flow directions to the fluid-porous interface Γ .Note that for small opening widths h the flow will have bigger angle to the interface (see Fig. 5 in [21]).Under the assumptions provided in Section 2.1, the Darcy flow equations can be applied to describe fluid flow through the porous medium The boundary conditions on Γ pm = ∂Ω pm \ Γ read where v pm is the Darcy velocity, p pm is the pressure, q is the source term, n pm is the unit outward normal vector from the domain Ω pm on its boundary, Γ pm = Γ D,pm ∪ Γ N ,pm , and p pm , v pm are given functions.The intrinsic permeability tensor K is symmetric, positive definite and bounded.Depending on their morphology, porous media can be classified as -isotropic, when K = kI with k > 0, -orthotropic, when K = diag{k 11 , . . ., k dd }, where d is the number of space dimensions, -anisotropic, with a full tensor K = k i j 1≤i, j≤d , where k i j = 0.
In this paper, we consider homogeneous isotropic porous media with K = kI, k > 0. Permeability of periodic porous structures can be computed using the theory of homogenisation.It can be estimated using the porosity-permeability relationships, e.g.Kozeny-Carman, or determined by means of numerical upscaling techniques, e.g.[28,58,67] and references therein.However, the permeability is uncertain in general.
In order to obtain a closed problem formulation for the Stokes-Darcy problem, Eqs.1-4 have to be completed by an appropriate set of coupling conditions on the fluid-porous interface Γ .

Interface conditions
Different sets of interface conditions are proposed in the literature to couple the Stokes Eq. 1 with the Darcy flow model Eq. 3 depending on the flow direction and the application of interest, e.g.[4,21,32,40,50].In this paper, we consider the classical set of interface conditions Eqs. 5-7 and analyse the sensitivity of the coupled model Eqs.1-7 to the exact position of the fluid-porous interface Γ and the choice of material parameters.
The conservation of mass across the fluid-porous interface is given by and the balance of normal forces is The Beavers-Joseph coupling condition [8,33] for the tangential velocity reads where √ K = √ kI, n = n ff = −n pm is the unit normal vector on the interface pointing outward from the free-flow domain Ω ff , τ is a unit vector tangential to the interface (Fig. 1) and α BJ > 0 is the Beavers-Joseph slip coefficient.
The Beavers-Joseph parameter characterises pore-space morphology near the fluid-porous interface including the microscale interface roughness.Usually, the slip coefficient is taken to be α BJ = 1, however, this is not optimal for many flow problems.It was shown that the model behaviour is very sensitive to the choice of this parameter [20,57].
Some researchers distinguish between the permeability of the porous bulk used in Eq. 3 and the permeability of the nearinterface region used in the Beavers-Joseph condition Eq. 7, e.g.[40].In this case, the slip coefficient α BJ can be scaled.We consider the same permeability value in Eqs. 3 and 7, and analyse model sensitivity to the Beavers-Joseph coefficient α BJ and possible simultaneous effects of permeability k and parameter α BJ .

Model parameters
In this paper, we analyse sensitivity of the coupled Stokes-Darcy model Eqs.5-7 with respect to the following model parameters that are uncertain: permeability k, exact location of the fluid-porous interface Γ (x 2 = γ ) and the Beavers-Joseph slip coefficient α BJ .In addition, we study different flow directions which are controlled by the opening h at the outflow boundary (Fig. 1).
We choose the intervals for the Beavers-Joseph slip coefficient α BJ , interface location γ and permeability k as well as their distributions based on the data from the literature.We take a larger interval α BJ ∈ (0, 10) than in the original paper [8] in order to take into account a wider range of porous materials.Based on the most commonly used values in the literature, we set the mode of the probability density function (PDF) of the Beavers-Joseph slip coefficient α BJ = 1 and use a scaled beta distribution of the appropriate shape.The exact location of the interface Γ is typically set on the top of the first row of solid inclusions (γ = 0), however, some authors move this location above or below at a distance of a few solid inclusions.We choose the mode of the corresponding PDF and the parameters' range accordingly, and use again an appropriate beta probability distribution.
The values of the intrinsic permeability are taken from the interval k ∈ (10 −8 , 10 −5 ) to reflect commonly considered porous materials.Here, we use the log-normal probability distribution [56] to generate the initial dataset and reduce it to the interval of interest.The outflow opening width h is varied in order to simulate different flow directions to the fluidporous interface and it is related to the considered geometry.The uniform distribution is chosen for h since there is no preference in flow direction.The intervals for these uncertain parameters are summarised in Table 1.The histograms and related PDF for the uncertain parameters are plotted in Fig. 2. Using these assumptions we generate the dataset for the construction of related aPC/aMR-PC (piecewise) polynomial bases.

Introduction to aMR-PC
In this section, we provide a brief overview of the aMR-PC expansion, which is used as the surrogate model in this work.For a comprehensive description of the method we refer the reader to [36].The core idea of aMR-PC is to combine the data-driven polynomial basis construction provided by aPC [52] and the multi-resolution Ansatz [1,2,9].The latter reduces by localization the undesired oscillations of the response surface and compensates the reduction of the highest polynomial degree of the expansion.

Arbitrary polynomial chaos expansion
In the following, we assume that Y (x, θ(•)) ∈ L 2 (Ω) is a random field, depending on the random variable θ on the probability space (Ω, A, P) with the sample space Ω, σalgebra A and probability measure P.
In the classical (arbitrary) polynomial chaos Ansatz the random field Y ≡ Y (x, θ) can be expressed by the following expansion where N o is the maximal polynomial degree and • denotes the expectation E[•] with respect to the probability measure P. The response coefficients c p (x) are deterministic and can be computed as well by numerical quadrature as by regression.In present work, we use an appropriate least squares regression.In particular, in the context of UQ the response coefficients c p usually denote deterministic functions c p (x).
For the sake of brevity, we omit the function related notation and write, e.g., c p instead of c p (•).The polynomials ϕ p (•) of degree p are orthonormal with respect to the probability distribution of θ , i.e.
for p, q = 0, . . ., N o .For the common probability distributions, the orthonormal polynomial bases {ϕ p } p∈N 0 are already known.E.g., the Hermite polynomials are used for the Gaussian distribution and the Legendre polynomials are applied for the uniform distribution.For more comprehensive description of the construction and properties of orthonormal polynomials we refer to [5,15,22,24,63,64].
The aPC Ansatz introduced in [52] allows to generate the orthonormal polynomial bases from arbitrary data samples using properties of the stochastic moments without any additional knowledge on the related probability distribution.The major requirement for the construction of the aPC expansion is the existence and boundness of the stochastic moments In this case monic polynomials of type form the orthogonal basis.Here, we obtain the coefficients p i by solving the system of linear equations Orthonormal basis {ϕ k } k∈N 0 can be obtained by dividing monic polynomials by related normalizing coefficient √ ϕ k , ϕ k that also can be obtained without use of quadrature as following: For a comprehensive introduction into aPC we refer to [52].
In Fig. 3, we present the aPC bases for the datasets introduced in Section 2.4 and visualised in Fig. 2.These plots demonstrate a strong relation between the underlying data and the shape of the related polynomial basis.

Extension to the multi-resolution framework
The main idea behind the multi-resolution/multi-element framework is the decomposition of the stochastic domain into subdomains.In the present work, we use the dyadic decomposition with respect to the (empirical) cumulative density function F θ (η) := P(θ ≤ η) for η ∈ R. We assume that F θ (•) is a continuous, strictly monotonically increasing function over the interval (a, b) for a < b, a, b ∈ R with F θ (a) = 0 and F θ (b) = 1, e.g.[36,43].These assumptions provide the existence of the unique inverse F −1 θ such that the quantile function can be used as a numerical approximation of the inverse F −1 θ .This allows us to decompose the sample space/dataset Ω into stochastic subdomains (SD) by using dyadic decomposition To obtain the piecewise polynomial basis we construct the aPC polynomial basis on each SD such that the orthonormality relation on the whole stochastic domain Ω is satisfied For practical implementation it could be more convenient to work with the piecewise polynomial basis ψ N r l, p satisfying the orthogonal relation This is, e.g., the case, if one applies the aPC-framework introduced in [52] on samples contained in the SD Ω N r ,l .In Fig. 4, we present the aMR-PC bases for the datasets introduced in Section 2.4 and visualised in Fig. 2. In Eqs. 10 and 11, the zero-order polynomials are given by functions that are constant on the related SD and zero elsewhere

Multivariate formulation
The extension to the multi-dimensional case requires some additional formalism.Let M be the stochastic dimension, representing the number of uncertain model parameters.For the sake of readability, we assume the same refinement level N r ∈ N 0 in each direction.To address each SD we use the concept of multi-index set.Let I N r = 0, . . ., 2 N r − 1 be a set of indices related to SDs Ω N r ,l , l = 0, . . ., 2 N r − 1, such that Ω = l∈I Nr Ω N r ,l .Then let I M N r be the set of multiindices denoted by l = (l 1 , . . ., l M ) with l i ∈ I N r for all i = 1, . . ., M. Applying this to multi-dimensional SDs on the refinement level N r , we obtain For the extension to the polynomials in M variables θ = (θ 1 , . . ., θ M ) we use the standard multi-index notation α Now we define the space of the multivariate piecewise polynomial functions P M,N r N o spanned by the orthonormal multivariate polynomials with maximal total polynomial degree N o constructed as else.
Here, we use the reverse-lexicographic ordering in multivariate polynomial degree α for the construction of the basis.The multivariate piecewise polynomial functions N r l,α are orthogonal, i.e.
due to orthonormality of ψ N r l, p in each variable, as in Eq. 10.Use of the piecewise polynomials ψ N r k, p satisfying the orthogonality relation Eq. 11 leads to The total number of basis functions on the refinement level N r after truncation on the total polynomial degree N o is Further, we denote the expansion of the second order random field Y (θ ) by In present work, we obtain the polynomial coefficients by an appropriate least-squares approach on Gaussian quadrature points [24, chap. 3.1] related to the polynomial degree N o +1.
(Gaussian) quadrature can be used as well.The mean and variance of the truncated aMR-PCexpansion of the random field Y can be computed as

Global sensitivity analysis
The variance is the common measure for the uncertainty.

Sobol indices
The core idea of the Sobol orthogonal decomposition [27,[60][61][62][63] is to write the square integrable random field Y (θ 1 , . . ., θ M ) as a sum of the form where N := {1, . . ., M}.The decomposition of Y can be obtained by use of the operators Here • d P(θ) −{ j} denotes the integration over all variables except θ j .
According to [27,61], P I Y has the property This property leads to the orthogonality of the summands in decomposition Eq. 15 in the following sense This decomposition implies also the decomposition of the variance where The Sobol sensitivity indices are defined as follows The total sensitivity indices S T i are defined as

Sobol indices of the PCE
Before we proceed with the new results in Section 4.3, we would like to recapitulate the notation and calculation of the Sobol indices for the (global) PCE.One of the advantages of the PCE is the possibility to calculate the Sobol indices without additional evaluations of the response surface using the polynomial coefficients only.We will follow the notation introduced by Sudret [62].Let α = (α 1 , . . ., α M ) be a tuple corresponding to the polynomial degree of a multivariate polynomial α such that the second order random field Y (θ ) has a PC expansion Here we should keep in mind that the mean is given by the zero degree coefficient and the variance is computed from the coefficients related to the polynomial degrees higher than zero Using this property together with the orthonormality of the basis polynomials it is possible to compute the Sobol indices directly from the expansion coefficients Y α .More precisely, let us define the set of α tuples such that only the indices I = (i 1 , . . ., i s ) ⊆ N are nonzero This means that A N o i corresponds to the multivariate polynomials depending only on the random variable θ i .Now we can write the Sobol decomposition for the PCE of Y : For I = (i 1 , . . ., i s ) ⊆ N the PC based Sobol indices of a random field Y are defined by The total sensitivity index is given by

Extension of Sobol indices to the multi-resolution formulation
Now we extend the concept of the Sobol indices to the multi-resolution formulation.The main difference from the classical setting is the fact that in the multi-resolution setting also the polynomial of degree zero contribute to the variance.Therefore, we define the restriction of the piecewise polynomial N r l,α on I ⊆ N : This definition allows to restrict the contribution of the zero degree coefficients related to i ∈ N \ I to the mean and avoid their contribution to the variance.Due to orthonormality of N r l,α we have Since in the aMR-PC expansion also the zero degree polynomials contribute to the variance, we have to redefine the set of polynomial degree tuples α corresponding to the index set I : According to the definition of N r l,α,I given in Eq. 18, we have for α Here # (N \ (I ∩ J )) denotes the count of elements in N \ (I ∩ J ) and Let us define the main effect of Using the orthonormal relation of the aMR-PC piecewise polynomials, we obtain 123 The interaction of the components I ⊆ N is given by Here we should keep in mind that P ∅ Y M R = μ M R .Using induction on the lengths of the set I , we can transform formulation Eq. 20 into l,α,I (θ ) for constants C I J = (−1) #(I \J ) , J I N .It could be easily shown by induction that property Eq.16 is also valid for P I Y M R that means Sketch of the proof: Assume that the statement is valid for all subsets of N of length ≤ l.Let I ⊂ N , # (I ) = l + 1, j ∈ I .We set J = I \ { j}.Then, we have Using the argumentation from [27,61,62], we obtain for I = J , I , J ⊆ N .Now we can write the Sobol decomposition for the aMR-PC expansion of the random field Y M R using an elegant recursive formulation of P I Y M R given in Eq. 20: In the next step, we consider the computation of the variance σ 2 I of P I Y M R : Using orthogonality Eq. 22 and recombining Eq. 20, we obtain Therefore, due to Eqs. 22 and 19 the variance σ 2 I of P I Y M R for index set I ⊆ N is given by Using similar arguments as in Eq.21, we obtain the following formulation for the variance index Also here we should keep in mind that the term related to J = ∅ is equal to μ 2 .Now, using Eq.23 or Eq. 24, we can compute the Sobol sensitivity index SY M R i and the total sensitivity index SY T ,M R i for the aMR-PC expansion by the following formulae

Numerical validation of Sobol indices
In this section, we validate the expression of the Sobol sensitivity indices SY M R i and the total sensitivity indices SY T ,M R i of the aMR-PC expansion provided in Section 4.3.To approach the main application scenario of the aMR-PC: random distributions and datasets for which the higher polynomial degrees of the polynomial expansions like aPC are not suitable by several reasons, e.g., vanishing stochastic moments, we use the lowest possible nonlinear polynomial degrees N o = 2, 3 and vary the resolution level N r .For the sake of completeness, we also provide the results computed on the aPC [52] (without multi-resolution) for several polynomial degrees generated on the same dataset.
To validate accuracy of computation of the Sobol indices of the aPC and aMR-PC expansions we use, similar to [62], the so-called Ishigami function Here, X i , i = 1, 2, 3 are the over the interval [−π, π] uniformly distributed random variables and a, b ∈ R are deterministic parameters.The Ishigami function is widely used in the literature to benchmark sensitivity methods and it could not be easily captured by polynomial representation due to its nature.The most relevant properties of the Ishigami function are nonlinearity, non-monotonicity and the known formulae to compute the variance and the Sobol sensitivity indices analytically For the numerical experiments we use the values a = 7 and b = 0.1.Due to our main interest on the data-driven sensitivity analysis we use the randomly generated dataset of 10 6 samples for the construction of the aPC and aMR-PC expansions, and evaluate the Ishigami function Eq. 25 on the related roots.In Table 2, we provide the analytically computed Sobol sensitivity indices S I according to Eq. 26 and the total sensitivity indices S T I of the Ishigami function for all possible combinations of the input parameters I .In Fig. 5, we visualise the comparison of accuracy of the total sensitivity indices of the aMR-PC and aPC surrogate models.Here, we consider the relationship between the accuracy of the approach and two key factors: the number of polynomial coefficients, which determines the problem size, and the number of training samples, which determines the quantity of PDE-model evaluations needed.Since the Ishigami function is smooth, the most accurate and less training samples demanding approach can usually be obtained by high polynomial degree N o and lowest possible number of MR refinements N r .Nevertheless, real-world applications of aMR-PC are usually related to low N o and increasing N r , due to the higher reliability.Therefore, to validate this type of approach we compare the results with aPC that is equivalent to aMR-PC with N r = 0.
For the comprehensive comparison we refer the reader to the tables in Appendix.In Table 4, we present the Sobol sensitivity indices SY I and the total sensitivity indices SY T I for the aPC expansion of the Ishigami function for N o = 3, 5, 7, 9, 11, 13 as well as their absolute errors in comparison to the reference solution shown in Table 2.The obtained accuracy is comparable to the results presented by Sudret [62].However, also for the comparatively high number of used random samples in the dataset the accuracy can slightly vary for each new realisation of the random dataset.2.
This numerical experiment demonstrates the convergence of the Sobol sensitivity indices SY M R i and the total sensitivity indices SY T ,M R i of the aMR-PC expansion of the Ishigami function for increasing N r .Accuracy of the Sobol indices computed using aMR-PC expansion, that could be achieved by increasing the number of refinements N r , is at least comparable with the accuracy of the approximation of the Sobol sensitivity indices SY i and SY T i provided by the aPC expansion of higher polynomial degree, which also shows the convergence to the analytical solution for increasing polynomial degree.Partially higher computational costs and complexity of the aMR-PC expansion is an acceptable price for use of the PCE based framework if a global PCE expansion is not suitable for higher polynomial degree.

Global sensitivity analysis for coupled Stokes-Darcy flow problems
In this section, we apply the methods proposed in Sections 3 and 4 to the model introduced in Section 2. In this setting, we are interested in the impact of uncertain model parameters on the velocity distribution v = (u, v).Therefore, we will proceed as follows.In Section 5.1, we discuss accuracy of the aMR-PC surrogate models of the tangential and normal components of velocity u and v, respectively.We consider the free-flow domain and the sharp fluid-porous interface Γ = [0, 2] × {γ }, i.e., H = 1, B = 1 and L = 2 (Fig. 1).The inflow boundary conditions in Eq. 2 on the top boundary and side boundaries are taken where the top boundary is , and the outflow boundary is Γ N ,ff = {2} × (γ , h].The intervals for the outflow opening width h and the exact position of the interface γ are given in Table 1.
We consider the no-flow boundary conditions for the porous-medium domain, i.e., in Eq. 4 we have where the no-flow boundary is The dynamic viscosity is μ = 10 −3 [Pa s] and the gravitational effects are neglected both in the free-flow and porous-medium models, i.e. g = 0 in Eqs. 1 and 3.The typical values and the considered intervals for the permeability tensor K = kI in Eqs. 3 and 7 as well as for the Beavers-Joseph slip coefficient α BJ in the interface condition Eq. 7 are provided in Table 1.
We discretize the coupled Stokes-Darcy model Eqs.1-7 using the second-order finite volume method on staggered grids that are conforming at the fluid-porous interface.We use the Cartesian mesh and consider 201 × 202 control volumes in the coupled domain Ω ff ∪ Ω pm .This means that we have to compute 40 602 piecewise polynomial expansions, one for each node in the spatial mesh.However, in this paper we omit use of sparsity or adaptivity concepts and work with the full-tensor expansion to avoid possible influence of the adaptivity related artefacts on the results.Finally, in Section 5.2, we apply the aMR-PC surrogate models to compute the Sobol indices which describe the global sensitivity of the surrogate model.

Uncertainty quantification
Mean (E) and variance (σ 2 = V) or standard deviation (σ ) are the moments which are commonly used for description of random fields.For the validation of the aMR-PC surrogate model we compare mean and standard deviation of the aMR-PC expansion of the tangential velocity u and the normal velocity v using Eq. 13, Eq. 14 with the reference mean and standard deviation obtained from the MC approach for M = 50 000 realisations.To ensure a realistic selection of (N r , N o ) we have chosen to focus on a polynomial degree N o = 2 for aMR-PC.Due to our experience, the lowest nonlinear polynomial degree N o = 2 is a very common compromise between reliability and accuracy if no additional information is available.Another consideration is the limitation imposed by the number of model runs.To address this, we have set the number of refinements N r to 2, which allows for a reasonable number of samples.Furthermore, assessment of a more accurate approximate solution will also require a more accurate reference solution.In Fig. 6, we visualize the computational effort in terms of number of training samples and number of polynomial coefficients that is necessary to address a certain polynomial degree N o for given N r = 0, 1, 2 in case of four uncertain model input parameters (M = 4) (Fig. 7).
We visualise the results in Fig. 8 providing mean, standard deviation and log-variance of velocity v = (u, v) computed using the coefficients of the aMR-PC expansion for polynomial degree N o = 2 and refinement level N r = 2.We present mean, standard deviation and log-variance of the MC reference solution for M = 50 000 realisations in Fig. 7. Visually, the results of the aMR-PC approach conform with the MC results.Also the log-variance of u and v used to visualise values of variances, which are close to zero, shows the same behaviour.Beyond of the visual examination, we compute the spatial averaged L 2 -error of mean and standard deviation in comparison with the reference solution provided by MC and provide the results in Table 3.These results show the error decay for increasing polynomial degree N o as well as refinement level N r , where the aMR-PC expansion allows to obtain the same or slightly higher accuracy.The presented aMR-PC concept could be also seen as a generalization of the aPC representation, where the aMR-PC with N r = 0 is equivalent to pure aPC expansion.The pure aPC expansion can only handle the nonlinearity of the underlying problem by increasing the degree of expansion, which is only appropriate for smooth problems.Unfortunately, increasing the aPC expansion degree results in a substantial increase in the number of expansion terms that could be hardly handled in the post-processing (e.g.mentioned above 40 602 polynomial expansions).In contrast, the aMR-PC could also handle the nonlinearity of the underlying problem via decomposition onto the stochastic subdomains N r , without the need for a high expansion degree in each subdomain.Moreover, due to the decoupled structure of aMR-PC with N r > 0, each stochastic subdomain can be post-processed independently.For example, the low polynomial degree N o = 2 requires only 15 realizations (cf.Section 3 for details) to be processed simultaneously.
Now that we have shown that the aMR-PC-surrogate model describes the Stokes-Darcy flow problem with appropriate accuracy, we can proceed with the GSA in Section 5.2.

Global sensitivity analysis using Sobol indices
In this section, we use the Sobol sensitivity and total sensitivity indices for the global sensitivity analysis to describe the impact of the single uncertain parameter or parameter combinations on the model output of the surrogate model.
We visualise the distribution of the total sensitivity indices in space computed with aMR-PC for N o = 2 and N r = 2 in Fig. 9, where the total Sobol indices of both velocity components u and v are presented for the uncertain Beavers-Joseph slip coefficient α BJ , the interface location γ , the permeability k and the outflow opening width h.Spatial distribution of the total Sobol indices in Fig. 9 indicates that the interface  location and outflow opening significantly influence both velocity components.However, Fig. 9 reveals the detailed information about the relationship between the interface location, outflow opening and their impact on the spatial distribution of velocities.It appears that both velocity components demonstrate highly nonlinear behaviour in response to changes in the interface location and outflow opening, particularly in the middle of the coupled domain indicated by the sharp fronts in Fig. 9 (c, d, g, h).Horizontal sharp fronts are related to the common interface between two domains.The exact location of the fluid-porous interface mainly influences the velocity in the free-flow domain (Fig. 9 (c, d)).Since the free-flow velocity is much higher than the porous-medium velocity, the change in the domain size due to the interface position influences the velocity in the free-flow region much stronger in comparison to the porous medium.Since we consider low permeable media, fluid flow through the porous medium is more sensitive to permeability variations than the free flow (Fig. 9 (e, f)).The outflow opening controls flow direction in the free-flow domain.Therefore, it influences mainly fluid flow near the right boundary (outflow) and in the porous medium, since the direction of the flow to the interface is varying and thus affecting the infiltration velocity (Fig. 9 (g, h)).Additionally, the strong correlations between the interface location and outflow opening are illustrated in Fig. 9 (c, d) and (g, h).Here we should keep in mind that the Sobol indices (and also the total sensitivity indices, which are derived from them) are by definition normalised by the variance.Therefore, the results in areas where the variance is close to zero are not necessarily meaningful from the modelling point of view, since the influence of the numerical error can be emphasised.
To obtain a concise representation we use the space averaged sensitivity indices.In Tables 6 and 7  We note that also in this case the space averaged sensitivity indices computed using the aPC (N r = 0) and aMR-PC (N r > 0) expansions demonstrate similar behaviour.Visualisation of the space averaged sensitivity indices for N r = 2, N o = 2 is presented in Fig. 10 for both velocity components.From Fig. 10(c) and (d) we observe that the coupled Stokes-Darcy model Eqs.1-7 is very sensitive to the exact location of the sharp fluid-porous interface γ , the outflow opening width h and the permeability k.The Beavers-Joseph slip coefficient α BJ , at least with respect to the proposed probability distribution on the considered interval (Fig. 2(a) and Table 1), does not play such an important role for the velocity field in the whole coupled domain.However, it is essential for the tangential velocity component in the nearinterface region as can be seen from Fig. 9(a).In addition to contribution of parameters γ , h and k on the overall flow behaviour, the domains of relevance for these parameters on both velocity components can be seen from Fig. 9(c)-(h).Note that, as mentioned in the beginning of this section, the Sobol indices are in general not meaningful in domains, where the variance is close to zero, due to normalisation by the variance.More precisely, here it can not be excluded that the variance of the model uncertainty is closer to zero than the variance of the computational uncertainty/numerical error that is almost always inherently included in results of numerical simulations, with the consequence that the computational uncertainty becomes dominating in the results of GSA.Therefore, the sensitivity indices should be considered together with the variance, e.g., zero variance regions in Fig. 8 and the related regions in Fig. 9.

Discussion and conclusions
In this paper, we proposed a new method to compute the Sobol sensitivity indices and performed with them the GSA on aMR-PC surrogate models.This provides the flexibility and reliability which are necessary for realistic and datadriven models.First, we verified the numerical accuracy of the method on the Ishigami function, which is one of the standard benchmark problems for GSA, and compared the computed Sobol indices not only with known analytical solutions but also with the Sobol indices computed from the related PCE coefficients.
The goal of this work was to improve the applicability of the GSA for real world applications.We considered the coupled Stokes-Darcy flow model with four uncertain parameters and provided an overview over the application of the aMR-PC based GSA workflow.We extended the Stokes-Darcy problem formulation introducing uncertain model parameters and boundary conditions and proposed the appropriate probability distribution/dataset for each input parameter.In the next step, we developed and validated the (piecewise) polynomial based surrogate models for the coupled Stokes-Darcy flow problem.Then, we proceeded with the computation of the sensitivity indices of the aMR-PC surrogate model and the interpretation of the obtained results.We also provided the sensitivity indices computed from the PCE coefficients for the sake of comparison.
Analysing the total sensitivity indices we observed that the aMR-PC surrogate of the Stokes-Darcy flow model with the classical interface conditions is very sensitive to the exact position of the fluid-porous interface, the outflow opening width and the permeability.The Beavers-Joseph slip coefficient does not significantly influence the overall flow behaviour, at least for the proposed probability distribution on the considered interval.However, this parameter impacts the tangential velocity component near the interface.Since we consider low permeable porous media in this work, the porous-medium velocity is very slow.The permeability from the considered interval k ∈ (10 −8 , 10 −5 ) has almost no influence on the free-flow velocity.Therefore, the contribution of the permeability to the variance is much lower (in space average) in comparison to those of the interface location and the outflow opening.The contribution of the crosseffects seems to be proportional to the single contributions of the related uncertain parameters such that no unexpected cross-effects are observed.For example, the influence of the Beavers-Joseph coefficient α BJ is very low in general, therefore, the influence of the related cross-effects is also almost negligible.

Appendix B
An important question is how accurate the aPC and aMR-PC based surrogate models, provided in Sec. 5, approximate the numerical solution of the coupled Stokes-Darcy problem.Partially this question was already addressed by considering of the L 2 -error in mean and variance presented in Table 3.In Fig. 11, we compare the predictions of the velocity component u made by evaluation of aPC/aMR-PC surrogate models on 10.000 (randomly sampled) input parameter sets with model output to assess the prediction performance of the used surrogate models.More precisely, Fig. 11 In all tests, aMR-PC achieves moderately but visible higher accuracy.Since aPC works mostly accurate but is also outperformed by aMR-PC already on similar number of training samples, this can be also interpreted as an indication for small but measurable presence of heterogeneity's in the model response.

Fig. 1
Fig. 1 Schematic geometrical setting for the coupled problem

Fig. 2
Fig. 2 Probability density functions and histograms of the uncertain parameters

Fig. 3
Fig.3aPC bases for probability distributions shown in Fig.2

Fig. 6
Fig. 6 Computational effort in terms of (a) number of training samples and (b) number of polynomial coefficients to obtain a polynomial degree N o = 2, . . .8 for N r = 0, 1, 2 in case M = 4

Fig. 7
Fig. 7 Mean E, standard deviation σ and log-variance log σ 2 of tangential (left) and normal (right) velocities computed with MC for M = 50 000

Fig. 9
Fig. 9 Total Sobol indices of velocity (u, v) for uncertain Beavers-Joseph parameter α BJ , interface location γ , permeability k and outflow opening h computed with aMR-PC for N r = 2 and N o = 2 the space averaged sensitivity indices of velocity components u and v are presented.We use the following uncertain parameter indices: I = 1 for the Beavers-Joseph slip coefficient α BJ , I = 2 for the interface location γ , I = 3 for the permeability k, I = 4 for the outflow opening width h.

Fig. 10 Table 4
Fig. 10 Space averaged Sobol and total sensitivity indices of velocity: (a), (b) mean Sobol indices; (c), (d) mean total sensitivity indices (a) shows the scatter plot for one arbitrary chosen point in space for aPC/ with N o = 8 and aMR-PC with N r = 2, N o = 2.In Fig. 11(b) we present the LogLog plots of Mean Squared Error (MSE) over all samples and spacial expansions vs. number of training samples/model evaluations for aPC on N o = 2, 4, 6, 8 and aMR-PC for N o = 2, N r = 1, 2.

Fig. 11
Fig. 11 Surrogate model accuracy comparison aPC vs. aMR-PC on 10.000 evaluation samples

Table 1
Typical values and intervals for uncertain parameters for the Stokes-Darcy model Eqs.1-7

Table 2
Analytical Sobol indices of the Ishigami function Eq

Table 3
Spatial averaged L 2 -error in mean E[•] and standard deviation σ [•] computed by aPC and aMR-PC in comparison with reference solutions provided by MC with 50 000 realisations aPC Error