Sensitivity Analysis and Quantification of the Role of Governing Transport Mechanisms and Parameters in a Gas Flow Model for Low-Permeability Porous Media

Recent models represent gas (methane) migration in low-permeability media as a weighted sum of various contributions, each associated with a given flow regime. These models typically embed numerous chemical/physical parameters that cannot be easily and unambiguously evaluated via experimental investigations. In this context, modern sensitivity analysis techniques enable us to diagnose the behavior of a given model through the quantification of the importance and role of model input uncertainties with respect to a target model output. Here, we rely on two global sensitivity analysis approaches and metrics (i.e., variance-based Sobol’ indices and moment-based AMA indices) to assess the behavior of a recent interpretive model that conceptualizes gas migration as the sum of a surface diffusion mechanism and two weighted bulk flow components. We quantitatively investigate the impact of (i) each uncertain model parameter and (ii) the type of their associated probability distribution on the evaluation of methane flow. We then derive the structure of an effective diffusion coefficient embedding all complex mechanisms of the model considered and allowing quantification of the relative contribution of each flow mechanism to the overall gas flow. Relative importance of parameters driving gas flow in low-permeability media is assessed. The influence of parameter probability distribution on gas flow statistics is appraised. A simple effective diffusion model embedding major methane flow mechanisms is derived. Relative importance of parameters driving gas flow in low-permeability media is assessed. The influence of parameter probability distribution on gas flow statistics is appraised. A simple effective diffusion model embedding major methane flow mechanisms is derived.


Introduction
Methane is recognized as a potential energy source to assist transition to a carbon-free energy landscape (Hughes 2013), considerable reserves of methane being associated with subsurface reservoirs worldwide (U.S. Energy Information Administration 2015). After its generation, this gas typically accumulates in reservoir regions subdued to low-permeability layers (i.e., caprocks) that prevent its upward migration (Dembicki-Jr. 2017). Due to the partial sealing efficiency of caprocks, some amount of gaseous phase hydrocarbons might cross such barriers and reservoir gas can then be released into the overburden to (eventually) reach the surface (Schlömer and Krooss 1997;Schloemer and Krooss 2004). In this context, appropriate modeling approaches to quantify gas migration in low-permeability geomaterials can assist the appraisal of the feasibility of a methane recovery project.
A variety of models depicting gas movement in low-permeability geomaterials have been proposed (Wu et al. 2016;Sun et al. 2017;Rani et al. 2018;Wang et al. 2019). These models typically estimate the mass flow rate of gas as the result of a combination of various gas transport mechanisms taking place across the porous system. Parameters associated with these models, describing chemical, mechanical, flow, and transport features governing feedbacks between gas and the host rock matrix are always affected by uncertainty. The conceptual model of Wu et al. (2016) depicts the mass flow rate of a gas across a low-permeability medium as the sum of three key processes: (i) a surface diffusion, and two weighted bulk diffusion components corresponding to (ii) slip flow and (iii) Knudsen diffusion. This model takes into account changes in the porous system caused by mechanical deformation and adsorption/desorption dynamics. The model embeds numerous parameters which are typically estimated through (direct or indirect) laboratory-scale experiments. Considering the set of complex mechanisms involved, these types of experiments are costly, time demanding, and their results are prone to uncertainty. The latter is also related to the intrinsic difficulties linked to replicating operational field conditions at the laboratory scale as well as to the challenges stemming from transferability of results to heterogeneous field scale settings (Pan et al. 2010;Yuan et al. 2014;Tan et al. 2018).
Due to our still incomplete knowledge of the critical mechanisms driving gas movement in low-permeability media (Singh and Myong 2018;Javadpour et al. 2021) and the complexities associated with the estimation of model parameters, model outputs should be carefully analyzed considering all possible (aleatoric and epistemic) sources of uncertainty. In this sense, sensitivity analysis approaches are important tools enabling us to (i) quantify uncertainty, (ii) enhance our understanding of the relationships between model inputs and outputs, and (iii) tackle the challenges of model-and data-driven design of experiments (Dell'Oca et al. 2017). Hence, sensitivity analysis techniques may be effectively used in the context of methane flow modeling efforts to (i) quantify and rank the contribution of our lack of knowledge on each model parameter to the uncertainty associated with model outputs; (ii) identify model input-output relationships; and (iii) enhance the quality of parameter estimation workflows, upon focusing efforts on parameters with the highest influence to target model outputs (Saltelli et al. 2010;Dell'Oca et al. 2020). In cases where parameters associated with a model have already been estimated (e.g., through model calibration), the main purpose of a Global Sensitivity Analysis (GSA) is to assist quantification of the uncertainty still remaining after model calibration, thus guiding additional efforts for its characterization (e.g., Dell'Oca et al. (2020) and references therein). The probability density function (pdf) related to each model parameter at this step might differ from the one employed before model calibration and some model parameters might be associated with a reduced uncertainty. In cases where processes are described through black-box models, GSA can be employed to quantify the influence that the variability of hyperparameters embedded in these models can have on their outcomes. We note that if uncertainty of some model parameters is further constrained, for example through (stochastic) inverse modeling (e.g., Ceresa et al. (2021)), results of the uncertainty quantification might also change. In this work we illustrate the methodological framework and the workflow required for GSA of a methane flow model and provide the elements to perform such an analysis for diverse scenarios. In order to assist this process, we provide a repository with scripts developed during this work (see declaration section).
In this work, we rely on GSA approaches to study the behavior of the aforementioned gas migration model targeting low-permeability media. While previous works focus on only a few selected model parameters (Song et al. 2016;Wu et al. 2017;Sun et al. 2017), a comprehensive diagnosis of the system behavior based on rigorous and modern GSA approaches taking into account the way all model parameters influence model output uncertainty is still missing. Here, we do so by implementing two GSA techniques, respectively based on the evaluation of (i) the classical (variance-based) Sobol' indices (Saltelli and Sobol' 1995) and (ii) the recent moment-based GSA metrics proposed by Dell'Oca et al. (2017). We recall that GSA approaches relying on Sobol' indices are widely used to quantify the relative expected reduction of variance of the target model output due to the knowledge of (or conditioning on) a given parameter. These have been employed in several applications including diagnosis of models related to, e.g., flood risk assessment (Koks et al. 2015), overpressure risk assessment in sedimentary basins (Colombo et al. 2017), and energy storage (Xiao et al. 2021). A critical limitation of variance-based GSA methodologies is that the uncertainty of the output is considered to be completely characterized by its variance. Such an assumption can lead to an incomplete characterization of the system behavior. The moment-based GSA approach introduced by Dell' Oca et al. (2017) is designed to enhance our capability to evidence model behavior upon including the quantification of model parameter uncertainty on the (statistical) moments of the pdf of a model output of interest. As such, this comprehensive approach yields information enabling us to characterize various aspects of uncertainty, without being limited solely to the concept of variance. The ensuing indices (termed AMA indices, after the initials of the authors (Dell'Oca et al. 2017)) have been effectively employed in a variety of contexts, including geophysical analyses related to gravimetric responses due to pumping tests (Maina et al. 2021), biochemical degradation of compounds such as glyphosate in soils (la Cecilia et al. 2020), and groundwater flow, including its feedbacks with evapotranspiration (Bianchi Janetti et al. 2019;Maina and Siirila-Woodburn 2020).
This work is organized as follows: Sect. 2.1 briefly illustrates the complete model we consider to describe methane flow in low-permeability media. The main theoretical elements of the GSA approaches employed are described in Sect. 2.2. Key results of the GSA are presented in Sect. 3, where we also derive and discuss novel formulations describing an effective diffusive behavior and encapsulating all physical-chemical mechanism included in the full methane flow model described in Sect. 2.1. Finally, conclusions are drawn in Sect. 4.

Gas Flow in Low-Permeability Media
Models adopted to quantify gas migration in low-permeability media can be classified according to their complexity, in terms of, e.g., conceptualization and mathematical rendering of the embedded processes, as well as number of their characteristic parameters. Among existing models associated with a high degree of complexity and including multiple transport processes jointly contributing to the total gas migration across the system (Mehmani et al. 2013;Wu et al. 2015aWu et al. , 2016Wu et al. , 2017Sun et al. 2017;Zhang et al. 2018;Javadpour et al. 2021), here we consider the model of Wu et al. (2016). The selected model allows considering mechanical deformation as well as relevant features associated with real gases such as variations in the gas viscosity ( ), and the effects of the compressibility ( C g ) and gas deviation (Z) factors caused by pressure and temperature changes.
The model introduced by Wu et al. (2016) rests on a conceptual picture according to which the total mass flow rate of gas per unit of area (J) is rendered through the sum of (i) a surface diffusion ( J s ) and two weighted bulk diffusion components, corresponding to (ii) slip flow ( J v ), and (iii) Knudsen diffusion ( J k ), i.e., The surface diffusion component is given by Wu et al. (2015b) where p is (gas) pore pressure and p l represents the strength of the driving force through the system, corresponding to the spatial gradient of gas pore pressure. The (dimensionless) coefficient ms is intended to take into account the possibility of applying the model (originally developed for capillary tubes) to a complex pore space and is defined in Eq. (17) of the Appendix where it is shown that ms depends on porosity ( ), tortuosity ( ), pore size (r) (i.e., pore radius), and gas coverage on the geomaterial ( ). The term D s in Eq. (2) is the surface diffusion coefficient, which is expressed (as shown in Eq. (25)) in terms of gas temperature (T), isosteric adsorption heat of the geomaterial ( ΔH ), a parameter ( ) related to the blockage/migration ratio of the adsorbed molecules, and . Finally, C sc , defined in Eq. (28), is the adsorbed concentration, which in turn depends on and on the gas molecule diameter ( d m ).
The model proposed by Wu et al. (2016) allows representing the mechanical deformation of the pore space (in terms of variation of permeability and porosity with pressure) through power-law relationships and making use of the classical Kozeny-Carman equation. Here, we rest on their original model formulation, which naturally leads to Equations (19) and (20), clearly evidencing that both r and evolve with p as a function of a reference pore radius ( r o ) and reference porosity ( o ), respectively.
The weight coefficients of the slip flow ( w v ) and Knudsen diffusion ( w k ) components in Eq. (1) are given by (Wu et al. 2016) Here, K n is the (dimensionless) Knudsen number defined as with where M and R are the gas molar mass and universal constant, respectively. Note that K n relates the mean free path of the gas molecules ( ) to a representative length of the system (Civan 2010), here taken as the pore diameter.

3
Here, mb is intended to take into account the possibility of applying the slip flow formulation (7) to a complex pore space (see Eq. (23)) and is the rarefied effect coefficient for a real gas which, according to Karniadakis et al. (2005), is evaluated through Eq. (24). The Knudsen flow component is dominant in systems where K n > 10 (Ziarani and Aguilera 2012) and is evaluated as (Darabi et al. (2012); Liu et al. (2016)) Here, D f represents the fractal dimension of the pore surface and denotes the ratio between d m and r.
We conclude by noting that the model here described includes a total of 15 parameters, which are related to the richness of physical processes embedded therein (See also Sect. 3.3). All quantities here introduced are listed in Table 1 and in the list of symbols and nomenclature Section.

Global Sensitivity Analysis
We perform a rigorous sensitivity analysis of the model illustrated in Sect. 2.1 to diagnose its behavior with reference to the estimate of methane flow as driven by imperfect knowledge of the associated parameters. Here, we note that uncertainties associated with the selection of the interpretative model is not analyzed. GSA can also be tailored to consider   Karniadakis et al. (2005) quantification of uncertainty of model outputs in the presence of uncertain interpretive models. In this context, uncertainty of a target variable which might result from the use of a collection of interpretive (conceptual and mathematical) models could be assessed upon relying, for example, on the approach illustrated by Dell'Oca et al. (2020). Our analysis is intended to yield a robust quantification of the relative importance of uncertain model parameters to a model output of interest. As mentioned in the Introduction, we rely on two GSA approaches, corresponding to (i) the classical variance-based technique grounded on the evaluation of the well-known Sobol' indices (Saltelli and Sobol' 1995) and (ii) the moment-based GSA framework introduced by Dell'Oca et al. (2017). Model parameters are treated as statistically independent, as the amount of available information does not enable us to clearly identify cross-correlations among parameters and to quantify joint distributions. We consider three differing characterizations of pdf describing uncertainty of model parameters: (a) all parameters are represented through uniform pdfs, (b) all parameters are represented by truncated Gaussian pdfs, and (c) the reference pore radius is characterized by a (truncated) log-normal pdf, while all remaining parameters are associated with uniform pdfs. Case a is representative of an approach where information on the considered parameters is limited so that all parameter values within the identified range of variability are equally weighted in the analysis (other studies relying on the same assumption include, e.g., Ciriello et al. (2013) 2020)). Case b is implemented as an alternative uninformed case, making use of the widely adopted hypothesis that model parameters are normally distributed. Case c takes advantage of the findings of Naraghi et al. (2018) who provide some experimental evidence suggesting that the pdf of pore radii in shales can be interpreted through a log-normal model. Our choice of performing sensitivity analyses according to configurations associated with diverse pdfs characterizing uncertain model parameters enables us to analyze the influence of model parameter pdf (which is generally unknown a priori) on the results of the GSA and, ultimately, on gas flow forecasting.
Considering the computational cost associated with multiple model evaluations (corresponding to 10 −4 seconds per simulation on an Intel Xeon Gold 6148 CPU @ 2.4 GHz) required for these analyses and the corresponding cost for random sampling across the considered high dimensionality model parameter space, our analyses rest on 10 8 model evaluations. The latter has been deemed to constitute an acceptable trade-off between the need to obtain stable results and computational efforts (details not shown). The pressure gradient acting on the system is set as a given boundary condition (and equal to 0.1 MPa/m) in all test cases.

Variance-Based Sobol' Indices
Sobol' indices (Saltelli and Sobol' 1995) can assist the appraisal and quantification of the relative expected reduction of the variance of a target model output due to knowledge of (or conditioning on) a given model parameter, which would otherwise be subject to uncertainty. In this context, considering a model output y, which depends on N random parameters collected in vector = (x 1 , x 2 , … , x N ) and defined within the space corresponding to the support of the ith parameter, x i ), the principal Sobol' index S x i associated with a given model parameter x i is evaluated as Here, E[⋅] and V[⋅] represent expectation and variance operators, respectively; the notation y|x i denotes conditioning of y on x i . Note that S x i describes the relative contribution to V[y] due to variability of only x i . Joint contributions of x i with other model parameters included in to the variance of y are embedded in the total Sobol' indices (details not shown). We recall that relying on Sobol' indices to diagnose the relative importance of uncertain model parameters to model outputs is tantamount to identifying uncertainty with the concept of variance of a pdf. As such, while Sobol' indices are characterized by a conceptual simplicity and straightforward implementation and use, they provide only limited information about the way variations of model parameters can influence the complete pdf of model outputs.

Moment-Based AMA Indices
The recent moment-based GSA approach proposed by Dell'Oca et al. (2017, 2020)  Here, AMAM x i represents the indices associated with a model parameter x i and a given statistical moment M of the pdf of model output y (considering the first four statistical moments of y, M = E for the mean, M = V for the variance, M = for the skewness, and M = k for the kurtosis). The AMA indices are intended to quantify the expected change of each statistical moment of y due to our knowledge of x i . Large values of these indices indicate that variations of the associated parameter strongly affect the statistical moments of y.

GSA of Methane Flow Model
The 15 uncertain model parameters of model (1) are considered to vary across the support defined through the ranges of variability listed in Table 1. These ranges have been designed upon considering available literature references (values typically employed for the model parameters in low-permeability geomaterials). With reference to three of the model parameters i.e., L p o , 1 , and , only very limited information is available from the literature, to the best of our knowledge (Karniadakis et al. 2005). Thus, we take the values considered by Wu et al. (2016) and Karniadakis et al. (2005) as the centers of corresponding ranges of variability. We then consider their (uniform) distributions to be characterized by a given coefficient of variation (that we set as 30%), thus enabling us to imprint these parameters with a sufficiently broad range of variability, which is also consistent with the degree of variability documented for the remaining uncertain parameters (see Table 1). Finally, we allow the fractal dimension D f to vary within its theoretical bounds (i.e., 2 < D f < 3) (Coppens 1999;Coppens and Dammers 2006). Methane properties (such as viscosity, compressibility, and deviation factor) are estimated using miniREFPROP (Lemmon et al. 2018), a tool that incorporates equations of state for a variety of gas species. With reference to methane miniREFPROP relies on the equation of state proposed by Setzmann and Wagner (1991). Table 2 lists the moment-based GSA indices related to mean (AMAE x i ), variance (AMAV x i ), skewness (AMA x i ), and kurtosis (AMAk x i ) of J as well as the principal Sobol' indices ( S x i ) evaluated for methane flow rate values rendered by Eq. (1) for the case in which all model parameters are modeled as independent and identically distributed random variables, each characterized by a uniform pdf (Case a).
While the strength of the influence of the reference pore radius ( r o ) on the model output is not the same for the (first four) statistical moments, the AMA indices clearly suggest that conditioning on r o has (overall) the strongest impact on the first four statistical moments of methane flow. This is then followed by reference porosity, pore pressure, tortuosity, and temperature. While the remaining model uncertain parameters still exert some influence on the (first four) statistical moments of J (as evidenced by the non-zero values of AMA indices), the strength of their influence can be considered as marginal when compared to the above mentioned quantities, which are seen to be key in driving the main features of the pdf of methane flow. In the following, we denote as most influential parameters for metrics AMAM , respectively. Parameters identified as most influential by each of these metrics are reported in bold in Table 2. Values of the Sobol' principal indices are generally Table 2 Moment-based GSA indices AMAM x i and Sobol' principal indices S x i for all x i parameters included in Eq. (1) All model parameters are described by uniform pdfs (Case a). Values of each metric identifying the most influential parameters are reported in bold consistent with the results stemming from the moment-based GSA, even as and T are not identified as influential to the model output according to the Sobol' principal index. This result is consistent with the observation that conditional variance can be larger or smaller than its unconditional counterpart (see also Fig. 1b) in a way that its integral over Γ T vanishes. A similar effect associated with the principal Sobol' indices was identified by Dell'Oca et al. (2017) with reference to the Ishigami function, which is a widely used analytical benchmark in sensitivity studies. Figure 1 depicts the first four statistical moments of J conditioned on values of the five most influential uncertain parameters selected on the basis of Table 2. Uncertain parameters are normalized to span the unit interval, for ease of interpretation. Unconditional moments are also depicted as a reference. We note that, when considering conditioning on the model parameters which have been identified as non-influential according to the metrics employed, the difference between conditional and unconditional moments is negligible (details not shown).
As expected, conditioning on values of the reference pore radius ( r o ) yields the most marked effects to all of the statistical moments considered (see black dotted curves in Fig. 1). Mean and variance of methane flow generally increase with r o . A minimum mean methane flow value is attained for 2 < r o < 15 nm (corresponding to the range of normalized values comprised between 0 and 0.15 in Fig. 1). The dominant transport mechanism for r o < 15 nm is surface diffusion, the strength of its contribution decreasing with increasing r o . As r o increases, the strength of the contribution related to surface diffusion decreases faster than the corresponding increase of the slip flow contribution, thus resulting in a minimum value for the expected methane flow for values of the reference pore radius comprised in the aforementioned range. Otherwise, skewness and kurtosis (i) are affected by variations of the reference pore radius when the latter is smaller than 20 nm (corresponding to a normalized value of 0.18); and (ii) are generally constant for r o > 20  Table 2): a expected value, b variance, c skewness, and d kurtosis. The corresponding unconditional moments (i.e., SM Y ) are also depicted (gray bold horizontal lines). Intervals of variation of the uncertain model parameters are rescaled within the unit interval for graphical representation purposes. All model parameters are described by uniform pdfs (Case a) nm. Nevertheless, we note that these (statistical) moments are still remarkably different from their unconditional counterparts even for large r o values, thus evidencing the impact of acquired knowledge on r o on reducing the asymmetry (as rendered by the skewness) and the peakedness and tailedness (i.e., the probability associated with extreme values, as rendered by the kurtosis) of the methane flow pdf. Conditioning on pore pressure imprints variations to the statistics of the model output which are qualitatively similar to those associated with r o . Larger values of mean and variance of J are linked to larger values of p. This result descends from the linear relationship between pore pressure and slip flow (Eq. (7)), the latter being the dominant mechanism in systems formed by larger pores. Conditional skewness and kurtosis are constant (albeit different from their unconditional counterpart) across most of the variability range of p, sharp variations of these quantities being associated with conditioning on low values of p (i.e., corresponding to pore pressure values smaller than 10 MPa). Our findings about the influence of p on J are consistent with the results of Sun et al. (2017). These authors find that increasing values of pore pressure lead to an increase of apparent permeability (which is in turn linearly proportional to gas flow) for r o > 10 nm. Wu et al. (2016) document a similar behavior due to the dominance of the slip flow component (which is proportional to p; see Eq. (7)) in systems characterized by large pores.
While the impact of reference porosity and tortuosity is not analyzed in any of the available previous studies (Wu et al. 2015b(Wu et al. , 2016(Wu et al. , 2017Sun et al. 2017;Zhang et al. 2018), our results rank o and as the second and fourth most influential parameters in the evaluation of the pdf of J, respectively (see Table 2). The correction factors for bulk (Eq. (23)) and surface (Eq. (17)) diffusion flow increase linearly with reference porosity. Thus, increased values of o yield corresponding increases of the methane flow (and hence of its first two statistical moments) independent of the dominant transport mechanism. Conditional mean and variance of J decrease with increasing values of tortuosity. This is in line with the observation that all gas transport mechanisms are characterized by an inverse proportionality between J and through the correction factor which is related to these processes taking place within a porous domain. These elements are consistent with a physical picture according to which fluid flow rates across a porous geomaterial are expected to increase and decrease with increasing porosity and tortuosity, respectively. Unlike pore pressure and reference pore radius, conditioning on reference porosity and tortuosity yields a reduction of skewness and kurtosis of the pdf of J, whose conditional values remain constant independent of the value of o and/or .
Conditioning on temperature (T) affects the mean and variance of the methane flow pdf in a way which is qualitatively similar to the effect of tortuosity (albeit quantitatively to a lesser extent) due to the inverse proportionality between J and T. Otherwise, the overall shape of the pdf of J is not significantly influenced by the knowledge of T, values of conditional skewness and kurtosis practically coincide with their unconditional counterparts.
The results listed in Table 2 suggest that statistical moments of methane flow are virtually insensitive to the remaining parameters (i.e., 10 of the 15 model parameters). Therefore, setting any of these parameters at given values within the variability space considered in our analysis yields only minor changes in the prediction of J. In this context, our results suggest that methane flow can be assessed with an acceptable degree of reliability even in the presence of scarce information about several parameters embedded in Eq. (1) such as, e.g., the overburden pressure (i.e., p c ), the power-law exponents related to porosity (i.e., q) and pore radius (i.e., t), the fractal dimension of the pore surface (i.e., D f ), or the isosteric adsorption heat of the geomaterial (i.e., ΔH ). Further to this, our results suggest the opportunity to prioritize allocation of resources to robust characterization of (in descending order) reference pore radius, reference porosity, pore pressure, tortuosity, and temperature.

Impact of the Model Parameter pdfs on GSA Results
In this section, we analyze the impact of the choice of model parameter distribution on the pdf of J. As described in Sect. 2.2, we compare the GSA outcomes obtained with a uniform pdf for all model parameters (Case a) and illustrated in Sect. 3.1 against those computed when (i) all model parameters are characterized through (truncated) Gaussian pdfs (Case b) and (ii) r o is described by a (truncated) log-normal pdf while the remaining parameters are described as in Case a (Case c). To provide a consistent comparison, Gaussian and lognormal pdfs are defined to honor the same mean and variance of the scenario associated with Case a. Table 3 lists values of AMA and principal Sobol' indices for each of the parameters embedded in Eq. (1) for Case b. Results of Tables 3 and 2 are qualitatively similar, i.e., the GSA yields similar results considering a uniform or a (truncated) Gaussian pdf for all model parameters. Our results imbue us with confidence about the documented ranking of parameter importance, with reference pore radius, reference porosity, pore pressure, tortuosity, and temperature identified as the model parameters being the key drivers to the evaluation of the major features of the pdf of methane flow. Values of statistical moments of J conditioned on model parameters for Case b are very similar to those depicted in Fig. 1 for Case a (details not shown). Table 4 lists the AMA and the principal Sobol' indices associated with J for Case c. In this case, it is even more evident that the uncertainty of r o is strongly dominant on the evaluation of the pdf of methane flow. Additionally, the blockage/migration ratio of the Table 3 Moment-based GSA indices AMAM x i and Sobol' principal indices S x i for all x i parameters included in Eq. (1) All model parameters are described by truncated Gaussian distributions (Case b). Values of each metric identifying the most influential parameters are reported in bold adsorbed molecules ( ) gains importance with respect to previous cases, quantitatively impacting the pdf of J to an extent which is similar to what exhibited by temperature. This feature is attributed to the abundance of small pores in this scenario, which favors the dominance of the surface diffusion flow mechanism (linked to parameter ). Figure 2 depicts the first four statistical moments of methane flow conditioned on values of influential uncertain parameters for Case c (see Table 4). Unconditional moments are also shown as a reference. Overall, the results are qualitatively similar to those embedded in Fig. 1 for Case a. The unconditional mean and variance of J in Case c are reduced (to approximately one-fourth and one-sixth, respectively) with respect to the corresponding values for Case a. Otherwise, unconditional skewness and kurtosis increase by about 2.6 and 6 times, respectively. These behaviors are attributed to the larger frequency of small reference pore radius values considered in Case c with respect to Case a (and b). Low values of reference pore radius are associated with large values of surface diffusion and to small values of mean and variance of methane flow. Conditioning on r o and o imprints variations to the model output mean and variance across the entire range of variability of these parameters (Fig. 2). We further note that conditioning on r o strongly reduces skewness and kurtosis of the pdf of J, thus reducing the probability associated with extreme (large) values of J.
Conditioning on p induces variations in the (first four) statistical moments of the model output. Conditioning on larger values of this quantity yields the highest values of mean and variance of the model output. A minimum in the values of conditional variance, skewness, and kurtosis is observed in the interval 1 MPa< p <15 MPa. Finally, the blockage/migration ratio of adsorbed molecules displays (a small but noticeable) influence on the model output pdf. Mean and variance of J decrease with increasing values of . This behavior is expected, given the nature of , high values of this parameter being related to significant blockage of gas molecules on the geomaterial surface. Reference pore radius ( r o ) is described by a (truncated) log-normal distribution and the remaining model parameters are described by uniform distributions (Case c). Values of each metric identifying the most influential parameters are reported in bold

Scaling of Gas Flow Model and Identification of Dominant Flow Mechanisms
A pure diffusion modeling approach has been shown to represent with an acceptable degree of accuracy the movement of methane in low-permeability media (Lu et al. 2015). Such a model embeds all physics governing the system dynamics in a unique parameter (i.e., a diffusion coefficient D) and, under steady-state conditions, the mass flow rate of methane can be expressed as where C∕ l represents the spatial gradient of methane concentration (C), i.e., the driving force of the system. Considering an isothermal system under single-phase flow and introducing the density of methane, = pM∕RTZ , Eq. (11) can be written as We complete our set of results and discussion by noting that the model illustrated in Sect. 2.1 coincides with a pure diffusion model (Eq. (12)) under single-phase conditions, as we illustrate in the following.
Equation (1) can be written as  Table 4): a expected value, b variance, c skewness, and d kurtosis. The corresponding unconditional moments (i.e., SM Y ) are also depicted (gray bold horizontal lines). Intervals of variation of the uncertain model parameters are rescaled within the unit interval for graphical representation purposes. Note that r o is described by a truncated log-normal pdf and the remaining model parameters are described by uniform pdfs (Case c) Comparing Eqs. (12) and (13) (16), where the contribution of each of the processes described in Sect. 2.1 is clearly recognizable. Figure 3 depicts color maps quantifying the relative strength of the contribution of the three flow mechanism (slip flow in red, Knudsen diffusion in green, and surface diffusion in blue) to the overall diffusion coefficient defined by Eq. (15) considering various combinations of all uncertain parameters embedded in Eq. (1) for all scenarios investigated. Each sub-plot depicts the average value of the ratio D i ∕D as a function of two parameters (i.e., averaging is performed with respect to uncertain parameters with the exception of the two varying along the (normalized) axes of the subplots), selected among those which were classified as most influential to the system (see Sects. 3.1 and 3.2).
Our results indicate that the dominant flow mechanism in defining the overall diffusion coefficient (and consequently the methane flow) is slip flow (in red in Fig. 3) in all of the analyzed cases. An exception is observed for small values of the reference pore radius and/ or small pore pressure, where surface diffusion is dominant. The contribution of Knudsen diffusion mechanism is always negligible. This suggests that it is possible to simplify Eq.
(1) by neglecting the Knudsen diffusion mechanism in the evaluation of methane flow.
Further simplifications of the methane flux model illustrated in Sect. 2.1 can be considered when the dominance of a given flow mechanism can be clearly established. For example, Fig. 3 suggests that the identification of the dominant flow mechanism is affected by the pdf of the uncertain model parameters. If r o is represented by a Gaussian (or uniform) pdf, J is mainly dominated by slip flow or surface diffusion with a sharp transition zone between these two mechanisms. Otherwise, when r o is represented by a log-normal pdf both mechanisms (i.e., slip flow and surface diffusion) may play an important role in the estimation of methane migration independent of the value of the model parameters.
Finally, we evaluate the pdf of the overall diffusion coefficient (D) by making use of Eqs. (15) and (16) for all scenarios analyzed. Sample pdfs as well as corresponding (14) B v = w v mb r 2 pM 8 ZRT (1 + K n ) 1 + 4K n 1 + K n , Maximum Likelihood (ML) fits of log-normal distributions are depicted in Fig. 4 in logarithmic and natural scales. Positive skewness and large kurtosis are evident for all cases, these being larger for Case c, as illustrated in Sect. 3.2. These results reinforce the observation of higher frequencies of low J values in Case c with respect to the other settings investigated. Sample statistical moments (mean, variance, coefficient of variation, skewness, and kurtosis) of the pdf of D are listed in Table 5 together with the parameters of the ML-based log-normal models. The overall diffusion coefficient can vary across about four orders of magnitude (i.e., between 10 −9 and 10 −5 m 2 /s). As expected, the largest variance of D is associated with Case a, where all parameters of model (1) are characterized by uniform pdfs. Otherwise, the largest coefficient of variation of D is associated with Case c. Finally, we remark that the results embedded in Fig. 4 can be of practical assistance, as they allow for fast evaluations of the probability that methane flow in low-permeability media exceeds a given threshold value.

Conclusions
We perform a rigorous Global Sensitivity Analysis (GSA) to assess the impact of uncertain model parameters on the evaluation of methane flow (J) in low-permeability media, such as caprocks. We study three scenarios that consider differing characterizations of the probability density function (pdf) describing model uncertain parameters to assess the impact of this choice on the results of the analysis. In details, we consider settings according to which (i) all model parameters are represented through uniform pdfs, (ii) all model parameters are represented through (truncated) Gaussian pdfs, and (iii) the reference pore radius is characterized by a (truncated) log-normal pdf, all remaining parameters being associated with uniform pdfs. Our work leads to the following main conclusions: 1. The uncertainty of methane flow is governed by uncertainty in the reference pore radius, followed (in decreasing order of importance) by reference porosity, pore pressure, tortuosity, temperature, and (to a lesser extent) blockage migration ratio of adsorbed molecules. The remaining parameters of the investigated model (Sect. 2.1) being practically uninfluential. This result can assist future efforts to allocate resources during experimental activities aimed at characterizing methane flow in caprocks.  can be related to a simple pure diffusion model by introducing an overall diffusion coefficient (D). The latter represented by the contribution of three effective diffusion coefficients, each associated with a well-defined flow mechanism. The ensuing mathematical structure of D allows distinguishing the relative contribution of all flow mechanisms to the overall methane flow. The relationship we derive also enables one to estimate the pdf of D when the model parameters are uncertain. The latter is a useful tool which can assist the probabilistic evaluation of J even in the absence of the detailed amount of information which is typically required to characterize the full methane flow model. 3. The shape of the pdf employed to characterize uncertain model parameters affects the results of our GSA. Additionally, it has a marked effect in the definition of the dominating transport mechanisms of the model. With reference to the model parameter variability considered in this study, as evaluated on the basis of available information, our results suggest that the dominant transport mechanism is slip flow. Surface diffusion plays also an important role, especially for low values of reference pore radius and pore pressure, while Knudsen diffusion is negligible in all of the test cases analyzed.

Appendix: Additional mathematical details related to the description of the gas flow model introduced in Sect. 2.1
The correction factor ms is given by with where r ad is thickness of the adsorbed gas layer, d m is gas molecule diameter, p c is overburden pressure, p o is atmospheric pressure, and is evaluated through a Langmuir equilibrium isotherm as where p L is a Langmuir pressure evaluated with The correction factor mb is expressed as The value of (in Eq. (7)) is evaluated through where uncertain parameters 0 , 1 , and allow representing the variation of as a function of the Knudsen number ( K n ). Here, 0 represents the maximum value of for large values of K n , 1 governs the values of for small values of K n , and is a shape parameter. The surface diffusion coefficient ( D s ) is given by with The adsorbed concentration ( C sc ) is given by where N A is the Avogadro Constant (6.02×10 −23 /mol).
Funding The authors acknowledge financial support from Geolog srl.
Data Availability All data used in the paper will be retained by the authors for at least 5 years after publication and will be available to the readers upon request.
Code Availability Codes used in this paper are available in the following github repository https:// github. com/ rlsan dovalp/ Sensi tivity_ Analy sis.

Conflict of interest Not applicable
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.