Time-varying sensitivity analysis for mixing in chaotic flows: a comparison study

TL;DR

This study compares different global sensitivity analysis methods applied to time-varying mixing in chaotic flows, highlighting their efficiency and reliability in complex, high-dimensional models.

Contribution

It provides a comparative evaluation of Sobol, Morris, and activity scores methods for sensitivity analysis in chaotic flow mixing, especially in high-dimensional settings.

Findings

Morris method is significantly more computationally efficient than Sobol.

All methods show comparable sensitivities in simpler flow models.

Morris and activity scores provide consistent results in high-dimensional models.

Abstract

Engineered injection and extraction systems that create chaotic advection are promising procedures for enhancing mixing between two species. Mixing efficiencies vary considerably, so carefully selecting the design parameters, like pumping rates, well locations, or operation times, is crucial. While numerous studies investigate the conditions required to achieve chaotic flow, sensitivity analyses addressing its impact on mixing have rarely been performed. However, selecting a suitable sensitivity analysis method depends on the underlying system and is often restricted by the computational cost, especially when considering complex, high-dimensional models. Moreover, the most appropriate metric to quantify mixing (e.g., plume area, peak concentration) can also be system-specific. We perform a time-varying sensitivity analysis on the mixing enhancement of two chaotic flow fields with different complexities. The rotated potential mixing (RPM) flow is parametrized using two or four hyperparameters, while the quadrupole flow utilizes 16 hyperparameters. We compare three global sensitivity analysis methods: Sobol indices, Morris scores, and a modification of the activity scores. We evaluate the temporal evolution of the sensitivity of the design parameters, compare the performance of the three methods, and highlight their potential in analyzing parameter interactions. The analysis of the RPM flow shows comparable sensitivities for all methods. Additionally, our numerical experiments show that Morris is the cheapest method, needing at most four times fewer model evaluations than Sobol to reach convergence. This motivates us to only use the computationally cheaper but as reliable Morris and activity scores on the 16-dimensional model, yielding again consistent results.

Tables

Table 1. Table 1 : Hyperparameter ranges for the sensitivity analysis. We sample the hyperparameters for each experiment according to the uniform distribution from the respective intervals.

Non-randomized RPM flow, large intervals
$\Theta$	$\tau$
$[0,\pi ]$	$[0.001,1]$
Non-randomized RPM flow, small intervals
$\Theta$	$\tau$
$[0.4\pi ,0.6\pi ]$	$[0.4,0.6]$
Randomized RPM flow, large intervals
$\overline{\Theta }$	${\Theta }_r$	$\overline{\tau }$	${\tau }_r$
$[0,\pi ]$	$[0,0.2\pi ]$	$[0.1,1]$	$[0,0.2]$
Randomized RPM flow, small intervals
$\overline{\Theta }$	${\Theta }_r$	$\overline{\tau }$	${\tau }_r$
$[0.45\pi ,0.55\pi ]$	$[0,0.1\pi ]$	$[0.45,0.55]$	$[0,0.1]$
Quadruple flow
$log\left(k_1\right)$	$log\left(k_2\right)$	$log\left(k_3\right)$	$log\left(k_4\right)$	$q_{n,s,e,w}$	$r_{n,s,e,w}$	${\theta }_{n,s,e,w}$
$[-0.122,0.354]$	$[0.354,0.697]$	$[0.698,1.042]$	$[1.042,1.539]$	$[0.95,1.05]$	$[0,5]$	$[0,2\pi ]$

Table 2. (a) Summary of the model geometry and flow and transport input parameters from [ ziliotto_2025 ] .

Geometry of the model
Aquifer size	$201\times 201\times 10$ ${\text{m}}^3$
Discretization	$\Delta x=\Delta y=0.5$ m
Aquifer thickness $b$	$10$ m
Boundary conditions	west & east: constant head, north & south: no flow
Distance from the origin to the ideal wells $L$	$25$ m
Random hydraulic conductivity field
mean ${\mu }_{\log \left(K\right)}$	$0.7$
variance ${\sigma }_{\log \left(K\right)}$	$0.25$
correlation length $\lambda$	$10$ m
Flow model
Porosity $ϵ$	$0.25$
Storage coefficient $S$	${10}^{-5}$
EIE system
Total duration of the EIE	$75$ days
Number of EIE steps	$12$
Duration of EIE step $T$	$6.25$ days
Injection rate $Q$ of the wells	${\Lambda }^2=\frac{QT}{\pi ϵb{L}^2}$, ${\Lambda }^2$ from Table 1(b)
Transport model
Size of contamination	$11\times 11$ ${\text{m}}^2$
Initial concentration	$1$ $\frac{\text{kg}}{{\text{m}}^3}$
Isotropic hydrodynamic dispersion coefficient $D$	$4\cdot {10}^{-2}$ $\frac{{\text{m}}^2}{\text{day}}$

Table 3. (a) Summary of the model geometry and flow and transport input parameters from [ ziliotto_2025 ] .

Geometry of the model
Aquifer size	$201\times 201\times 10$ ${\text{m}}^3$
Discretization	$\Delta x=\Delta y=0.5$ m
Aquifer thickness $b$	$10$ m
Boundary conditions	west & east: constant head, north & south: no flow
Distance from the origin to the ideal wells $L$	$25$ m
Random hydraulic conductivity field
mean ${\mu }_{\log \left(K\right)}$	$0.7$
variance ${\sigma }_{\log \left(K\right)}$	$0.25$
correlation length $\lambda$	$10$ m
Flow model
Porosity $ϵ$	$0.25$
Storage coefficient $S$	${10}^{-5}$
EIE system
Total duration of the EIE	$75$ days
Number of EIE steps	$12$
Duration of EIE step $T$	$6.25$ days
Injection rate $Q$ of the wells	${\Lambda }^2=\frac{QT}{\pi ϵb{L}^2}$, ${\Lambda }^2$ from Table 1(b)
Transport model
Size of contamination	$11\times 11$ ${\text{m}}^2$
Initial concentration	$1$ $\frac{\text{kg}}{{\text{m}}^3}$
Isotropic hydrodynamic dispersion coefficient $D$	$4\cdot {10}^{-2}$ $\frac{{\text{m}}^2}{\text{day}}$

Table 4. (b) Details of the EIE sequence or the quadrupole flow. The first row gives the step number, the second row indicates the well that is operated at this step (W: west, E: east, S: south, N: north), and the last row shows the value of Λ 2 \Lambda^{2} that controls the injection rate Q Q . Λ 2 > 0 \Lambda^{2}>0 indicates injection [ mays_2012 ] .

Step	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$	$10$	$11$	$12$
Well	W	E	W	E	W	E	S	N	S	N	S	N
${\Lambda }^2$	$\frac{3.5}{\pi }$	$\frac{3.5}{\pi }$	$-\frac{1}{\pi }$	$-\frac{3}{\pi }$	$-\frac{1.6}{\pi }$	$-\frac{1.4}{\pi }$	$\frac{3.5}{\pi }$	$\frac{3.5}{\pi }$	$-\frac{1}{\pi }$	$-\frac{3}{\pi }$	$-\frac{1.6}{\pi }$	$-\frac{1.4}{\pi }$

Equations

$$\frac{dx}{dt}=\frac{\partial \psi }{\partial y},\frac{dy}{dt}=-\frac{\partial \psi }{\partial x},\left(\begin{array}{c}x(0)\\ y(0)\end{array}\right)=\left(\begin{array}{c}{x}_0\\ y_0\end{array}\right)$$

$$\psi (x,y)={\tan }^{-1}\left(\frac{2x}{1-x^2-y^2}\right).$$

$$\begin{array}{rl}{t}_{\text{adv}}(\theta ,\psi )& ={\csc }^2(\psi )\\ \times \{& \cot (\psi )\arctan \left[\frac{\sin (\theta )\cot (\psi )}{\sqrt{1+{\cos }^2(\theta ){\cot }^2(\psi )}}\right]\\ & +\sin (\theta )\sqrt{1+{\cos }^2(\theta ){\cot }^2(\psi )}\\ & -|\cot (\psi )|(\theta +\cos (\theta )\sin (\theta ))\}.\end{array}$$

$$T_{\text{res}}(\psi )=\left(t_{\text{adv}}\left(\frac{\pi }{2},\psi \right)-t_{\text{adv}}\left(-\frac{\pi }{2},\psi \right)\right).$$

$$\begin{array}{rl}t& {}_{\text{adv}}(\theta (t),{\psi }_0)+\frac{T_{\text{res}}({\psi }_0)}{2}=\\ & \left(t_{\text{adv}}({\theta }_0,{\psi }_0)+\frac{T_{\text{res}}({\psi }_0)}{2}-t\right)\mathrm{mod}{T}_{\text{res}}({\psi }_0).\end{array}$$

$$\begin{array}{rl}\frac{dx}{dt}=& \frac{4xy}{(y^2+x^2-1{)}^2+4x^2}\\ \frac{dy}{dt}=& \frac{-2(x^2-y^2+1)}{(x^2+y^2{)}^2+2x^2-2y^2+1}\\ \left(\begin{array}{c}x(0)\\ y(0)\end{array}\right)=& \left(\begin{array}{c}0\\ y_0\end{array}\right)\end{array}$$

$$x(t)=0,\frac{dy}{dt}=\frac{2}{y^2-1},y(0)=y_0.$$

$$t_{\text{adv}}(y)=\frac{y}{2}-\frac{y^3}{6},T_{\text{res}}=\frac{2}{3}.$$

$$\begin{array}{rl}\left(\begin{array}{c}x(t_k)\\ y(t_k)\end{array}\right)& =\prod _{j=1}^k{R}_{j-1}^{+}{\Psi }^{{\tau }_j}{R}_{j-1}^{-}\left(\begin{array}{c}{x}_0\\ y_0\end{array}\right)\\ & =R[{\Theta }_n]\prod _{j=1}^{k}R[-{\Theta }_j]{\Psi }^{{\tau }_j}\left(\begin{array}{c}{x}_0\\ y_0\end{array}\right).\end{array}$$

$$S\frac{\partial h}{\partial t}=\nabla \cdot (T\nabla h)+W,h(x,0)=h_0(x)$$

$$\frac{\partial c}{\partial t}=-\nabla \cdot (vc)+\nabla \cdot (D\nabla c),c(x,0)=c_0(x)$$

$$V$$

$$V_i$$

$$V_{i,j}$$

$$S_i:=\frac{V_i}{V},S_{i,j}:=\frac{V_{i,j}}{V}.$$

$$\sum _i{S}_i+\sum _i\sum _{j>i}{S}_{i,j}+...+S_{1,2,...,k}=1.$$

$$\begin{array}{rl}{S}_{Ti}& :=\frac{\mathbb{E}[\text{Var}(f(x)|\{x_j:j\ne i\})]}{V}\\ & =1-\frac{\text{Var}(\mathbb{E}[f(x)|\{x_j:j\ne i\}])}{V}.\end{array}$$

$$\begin{array}{rl}{S}_{Ti}& =\sum _{s=1}^n\sum _{\begin{array}{c}{j}_1<...<j_s,\\ \exists k:j_k=i\end{array}}{S}_{j_1,...,j_s}\end{array}$$

$$\begin{array}{rl}{V}_i& =\frac{1}{N}\sum _{j=1}^{N}f(A{)}_j(f(A_B^{(i)}{)}_j-f(B{)}_j)\\ V& =\frac{1}{2N}\sum _{j=1}^{2N}{\left(f(AB{)}_j-\frac{1}{2N}\sum _{k=1}^{2N}f(AB{)}_k\right)}^2.\end{array}$$

$$V_{i,j}=\frac{1}{N}\sum _{j=1}^{N}f(B_A^{(i)}{)}_{j}f(A_B^{(j)}{)}_j.$$

$$S_{Ti}=\frac{\frac{1}{2N}{\sum }_{j=1}^N(f(A{)}_j-f(A_B^{(i)}{)}_j{)}^2}{V}.$$

$$d_i(x)=\frac{f(x_1,...,x_{i-1},x_i+\Delta ,x_{i+1},...,x_n)-f(x)}{\Delta }.$$

$${\mu }_i=\frac{1}{N}\sum _{j=1}^N{d}_i^{(j)},{\sigma }_i=\sqrt{\frac{1}{N-1}\sum _{j=1}^N(d_i^{(j)}-{\mu }_i{)}^2}$$

$${\mu }_i^{\ast }=\frac{1}{N}\sum _{j=1}^N|d_i^{(j)}|.$$

$${\mu }_i^{2\ast }={\left(\frac{1}{N}\sum _{j=1}^N|d_i^{(j)}|\right)}^2.$$

$$\begin{array}{r}{\hat{\mu }}_i^{\ast }=\frac{100\cdot {\mu }_i^{\ast }}{{\sum }_j{\mu }_j^{\ast }},{\hat{\sigma }}_i=\frac{100\cdot {\sigma }_i}{{\sum }_j{\mu }_j^{\ast }},{\hat{\mu }}_i^{2\ast }=\frac{100\cdot {\mu }_i^{2\ast }}{{\sum }_j{\mu }_j^{2\ast }}.\end{array}$$

$$C={\int }_{[0,1{]}^n}\nabla f(x)\nabla f(x{)}^{T}dx=W\Lambda W^T.$$

$$\hat{C}=\frac{1}{N}\sum _{j=1}^N\nabla \hat{f}(x_j)\nabla \hat{f}(x_j{)}^T.$$

$${\lambda }_i=w_i^{T}C{w}_i={\int }_{[0,1{]}^n}(\nabla f(x{)}^T{w}_i{)}^{2}dx,$$

$${\alpha }_i(k)=\sum _{j=1}^k{\lambda }_j{w}_{ij}^2,k\le n.$$

$${\nu }_i=\frac{1}{N}\sum _{j=1}^N{\left(\frac{\partial f(x_j)}{\partial x_i}\right)}^2.$$

Full text

1 Introduction

Mixing processes in laminar flow regimes are often dominated by molecular diffusion and are hence slow [stone_2001, stremler_2004, rolle_2019]. This is a major limitation for applications where fast mixing is required, like for chemical reactions in which the reactants need to come into contact, [neupauer_2014, borgne_2014, bertran_2023]. In such regimes, enhanced mixing can still occur through chaotic advection. The stretching and folding caused by this process generate complex flow paths that increase the interfacial area between the reactive species [aref_1984, aref_1990, liu_2000, stremler_2004]. This enhanced interfacial area leads to sharp concentration gradients, accelerating mixing and hence chemical reaction rates [borgne_2014]. Applications of chaotic advection to enhance mixing can be found, for example, for in-situ groundwater remediation [mays_2012, zhang_2017, ma_2022, neupauer_2014, ziliotto_, bertran_2023, bagtzoglou_2007, zhang_2009], microfluidics [aref_2017, ward_2015], gene expression profiling [stremler_2004], the synthesis of biological molecules [teh_2008], wastewater treatment [bagtzoglou_2006], biogas production, and organic removal efficiency in Anaerobic hydrogen production systems [huang_2018]. One can implement chaotic advection using a passive design, where the geometry of the system generates the chaotic flow field [stremler_2004, song_2003, liu_2000], or actively using external perturbations such as magnetic, pressure, or acoustic disturbances [aref_2017, bagtzoglou_2006] for microstirrers, acoustic mixers, or flow pulsation [ward_2015].

In Engineered injection and extraction (EIE) systems, mixing efficiency depends on the properties of the flow field, like pumping rates or well locations. Mixing can, for example, be limited by Kolmogorov–Arnold–Moser (KAM) islands in the flow, where fluids are trapped around elliptic points in the periodic flow field [lester_2009, chate_1999]. For different configurations of the same EIE system, mixing efficiencies can vary considerably [gubanov_2009, feistner_2024].

To ensure the effectiveness of these EIE systems for real-world applications, it is essential to assess the uncertainty of the system with respect to its design parameters [ciriello_2013]. Nonetheless, a time-varying sensitivity analysis of mixing for EIE systems has been rarely performed [ziliotto_2025], and different sensitivity analysis methods have not been tested on these systems. To assess temporal changes, time-varying sensitivity analysis methods are mandatory and require validation [sarrazin_2016]. Time-varying refers to an analysis of the sensitivity at different times of the model instead of aggregating the overall sensitivity into one single value [pianosi_2016a]. When using a time-varying approach, it is essential to use a quantitative sensitivity metric, as it allows for comparison of the results among different time-steps. The choice of an appropriate method for the sensitivity analysis depends heavily on the computational complexity of the problem at hand, the dimensionality of the parameter span, and the properties of the model [pianosi_2016, iooss_2015, chiogna_2024]. A given problem might allow for different sensitivity analysis methods, which, in general, do not yield identical results [constantine_2017]. It is hence advised to apply multiple sensitivity metrics to support the general conclusion of the analysis [pianosi_2016]. However, especially for complex models with a high-dimensional input dimension, this analysis can be computationally expensive, and therefore it is conducted in only a few analyses [sarrazin_2016, dai_2024]. Assessing the computational costs of the sensitivity metrics, we concentrate on the number of Monte-Carlo samples, as the execution of the model is usually the most expensive part of the sensitivity analysis [pianosi_2016].

One of the most popular methods for sensitivity analysis is the Sobol index [sobol_2001, saltelli_2010, sobol_1993] used beyond others in [mazziotta_2024, chiogna_2024, ciriello_2013, li_2022, khorashadizadeh_2017]. It is a variance-based method that relates the variance of the output to the variance of the input [saltelli_2010], assuming that the second moment is sufficient to describe the uncertainty of the model [saltelli_2002a]. Approximating the Sobol indices can be done using Monte-Carlo sampling [saltelli_2010]. The computational effort grows with the number of uncertain parameters, requiring $N\cdot(2n+2)$ model evaluations, where $N$ is the number of samples, and $n$ is the input dimension. If the second-order indices are not needed, the number of evaluations can be reduced to $N\cdot(n+2)$. High-dimensional models usually require more samples to converge than low-dimensional ones [pianosi_2016]. Further, the work of [borgonovo_2006, borgonovo_2011, pianosi_2015, pianosi_2016, delloca_2017] shows that Sobol indices may represent output uncertainty incorrectly if the model output is highly skewed or multi-modal. In cases, however, where Sobol indices are applicable, the results can lead to a good understanding of the model uncertainties.

Further, we look at the Morris method [morris_1991, campolongo_2007], which is also used to perform sensitivity analysis in groundwater remediation in [ma_2022, ziliotto_2025]. This method quantifies sensitivity by the mean of the elementary effects for each input parameter on the model response [morris_1991, campolongo_2007]. Using Monte-Carlo sampling, the method generates sample trajectories of length $n+1$ by changing each input parameter one at a time. The number of model evaluations also grows linearly with the number of uncertain parameters $N\cdot(n+1)$ ($N$ number of trajectories). Compared to the Sobol indices, far fewer samples are needed to reach convergence as the works of [ma_2022, ziliotto_2025] both yield realistic results with only $10$ trajectories.

One more recent sensitivity analysis method is the activity score that can be computed from the results of the active subspace method (ASM) [constantine_2017, bittner_2020, parente_2019a]. The ASM is based on the derivatives of the model output with respect to the uncertain parameters and can lead to low computational costs if gradients are available. If not, we can use a finite difference approach, leading to $N\cdot(n+1)$ samples to create $N$ gradients. In addition to the activity scores, the ASM also computes the active direction of the model, which gives insights into the important linear combinations in the model, and we can directly construct a surrogate model [constantine_2013, constantine_2015]. However, following [constantine_2017], activity scores may cause different results than ANOVA-based methods, like Sobol indices, for highly non-linear systems.

In this work, we assess which sensitivity analysis methods are suitable for systems influenced by a chaotic flow field, considering the dimensionality of the uncertain model input, the computational cost, the significance of the results, and the metric used to quantify mixing. As we are not interested in the sensitivity of the chaotic flow with respect to an initial location, we do not face problems of divergence of sensitivity metrics like [lea_2002, chandramoorthy_2017]. Instead, we focus on the overall effect of chaotic advection to enhance mixing. We perform a time-varying sensitivity analysis on two chaotic flow fields with different input dimensionality

. The RPM flow [lester_2009, metcalfe_2010] utilizes two and four input parameters that control the rotation angles and the time between rotations of one source and one sink. While the quadrupole flow [mays_2012] represents a modification to the pulsating dipole in [jones_1988], adding periodic reorientations of the well locations to enhance mixing. Its $16$ input parameters describe pumping rates, well locations, and hydraulic conductivities. Both flow fields find applications for groundwater remediation [bertran_2023, cho_2019, neupauer_2014, ziliotto_2025]. We use total Sobol indices, Morris scores, and a modified version of the activity scores from ASM that reuses the Morris samples to reduce computational costs. As the Morris and the activity scores are both qualitative methods, we introduce a scaling such that the results get a quantitative interpretation in terms of percentage contribution to the overall variability. For the RPM flow, we quantify mixing through the fraction of the domain covered by solute after a specific time, while for the quadrupole flow, we consider the maximal solute concentration. We compare the results of the different sensitivity analysis methods on both models and show how they are used to assess sensitivity and parameter interactions. Finally, we analyze the convergence rates of the sensitivity metrics.

Our paper is structured as follows: In Section 2, we start by introducing the chaotic flow systems that we use for our sensitivity analysis, the RPM flow in Section 2.1.1 and the quadrupole flow in Section 2.1.2. We then introduce the sensitivity analysis methods in Section 2.2, Sobol indices in Section 2.2.1, Morris scores in Section 2.2.2, and activity scores in Section 2.2.3. We show the results of our numerical experiments in Section 3, starting with the two versions of the RPM flow in Sections 3.1.1 and 3.1.2 and then continuing with the quadrupole flow in Section 3.2. Finally, we conclude our work in Section 4.

2 Methodology

2.1 Chaotic flow fields

2.1.1 The RPM flow, a low dimensional problem

Model setup

We consider the two-dimensional rotated potential mixing (RPM) flow that was introduced in [lester_2009]. The flow is produced by using one source and one sink that lie opposite each other on the edge of a unit circle, which serves as the domain for the flow. Both wells operate simultaneously, while fluid that is pumped out of the domain at the sink is reinjected instantaneously at the source. After operating the source and the sink for time $\tau$, the location of the wells is rotated instantaneously by the angle $\Theta$ around the origin. The rotation makes the flow unsteady, which is needed to generate chaotic advection [ottino_1989, stremler_2004, aref_1990, zhang_2009]. Our quantity of interest is $\bar{M}$, the fraction of the domain covered by particles initially placed in a tiny region at the top of the domain. As we do not consider diffusion, $\bar{M}$ is in theory constant in time. Due to our discretization, we add numerical diffusion when evaluating $\bar{M}$, leading to a monotonically increasing $\bar{M}$ over time. This choice is motivated in [feistner_2024] where the computation of $\bar{M}$ on a non-diffusive RPM flow is a good indicator for the effectiveness of chaotic mixing in the same configuration of the RPM flow with diffusion. For $t\rightarrow\infty$, the particles occupy the entire chaotic region of the domain, leaving only the KAM islands particle-free. For large $t$, the value of $\bar{M}$ hence approximates the size of the chaotic region.

Uncertain parameters

By construction, the RPM flow is parametrized by two flow parameters, $\Theta$ and $\tau$. We give a graphical representation of the system in Fig. 1(a). In addition to the standard RPM flow, we investigate a setup in which we modify the RPM flow by introducing randomization. Therefore, for each rotation the angle $\Theta_{i}$ is sampled from the uniform distribution $\mathcal{U}(\bar{\Theta}-\frac{\Theta_{r}}{2},\bar{\Theta}+\frac{\Theta_{r}}{2})$, similar we sample $\tau_{i}$ from $\mathcal{U}(\bar{\tau}-\frac{\tau_{r}}{2},\bar{\tau}+\frac{\tau_{r}}{2})$. The flow is hence described using four parameters: the means and the deviations of the respective randomization intervals. We visualize the randomized version of the RPM flow in Fig. 1(b). To avoid confusion with the flow parameters $\Theta_{i}$ and $\tau_{i}$, we will refer to $\Theta$ and $\tau$ of the non-randomized system and to $\bar{\Theta}$, $\Theta_{r}$, $\bar{\tau}$ and $\tau_{r}$ of the randomized system as hyperparameters. During our numerical experiments, we consider both versions of the flow. We give the ranges of the hyperparameters in Table 1.

Flow equations

The Lagrangian form of the RPM flow with source at $(0,1)^{T}$ and the sink at $(0,-1)^{T}$ can be described using the Hamiltonian system [sposito_1998, durst_2022]

[TABLE]

with Lagrange stream function [lester_2009]

[TABLE]

[lester_2009] derive the analytical solution to (1) with stream function $\psi$ in (2). For the convenience of the author, we briefly describe the process in the following. We focus first on the positive side of the $x$-axis (i.e., $x>0$) and represent the particle position in terms of the angle to the origin $\theta\in(-\pi/2,\pi/2)$ and its streamline $\psi\in(0,\pi/2]$. This parametrization is convenient, as $\psi$ is constant in time. The advection time of a particle along its streamline, until it reaches $\theta=0$, was established in [lester_2009] as

[TABLE]

This formula yields a negative value for particles downstream of $\theta=0$. Using (3), we can compute the residence time of a particle inside the domain from injection to extraction as

[TABLE]

Combining (3) and (4), we compute the angle $\theta(t)$ of a particle $(\theta_{0},\psi_{0})$ after time $t$ as

[TABLE]

The modulo operator $\bmod$ is needed to account for the instantaneous reinjection of a particle. We cannot specify the solution of (5) analytically, hence we use a linear interpolation of the inverse $\theta(t_{\text{adv}},\psi)$ on a $100\times 100$ grid. Together, this yields the flow $\Psi_{+}^{t}$ for all particles with $x>0$. For $x<0$, we use the reflection symmetry of $\psi$ with respect to the $y$-axis, yielding the flow $\Psi_{-}^{t}\left(\begin{smallmatrix}x\\ y\end{smallmatrix}\right)=\left(\begin{smallmatrix}-1&0\\ 0&1\end{smallmatrix}\right)\Psi_{+}^{t}\left(\begin{smallmatrix}-x\\ y\end{smallmatrix}\right)$.

For $x=0$, the usage of $(\theta,\psi)$ coordinates is insufficient. Alternatively, we directly solve (1)

[TABLE]

using the definition of $\psi$ in (2). Given that $x(0)=0$ we find $\frac{dx}{dt}|_{t=0}=0$ and hence $x(t)=0$ for all $t\in\mathbb{R}$. This simplifies (6) to

[TABLE]

We solve (7) for $y$ by using the separation of variables, leading to

[TABLE]

We can use $t_{\text{adv}}$ from (8) in (5) where $y$ takes the place of $\theta$ to compute the flow $\Psi_{0}^{t}$ for all particles with $x=0$. Combining the three flows, we can define the flow $\Psi^{t}$ for the complete domain.

To solve the equations when the source and the sink are rotated after time $\tau_{i}$, we introduce the rotation operator $R[\Theta]:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ that rotates the particles by the angle $\Theta$. The operators $R_{k}^{+}$ and $R_{k}^{-}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ are defined accordingly by $R_{k}^{+}=R[\sum_{j=1}^{k}\Theta_{j}]$ and $R_{k}^{-}=R[-\sum_{j=1}^{k}\Theta_{j}]$ with $R_{0}^{+}=R_{0}^{-}=R[0]$. Using $R^{-}_{i}R^{+}_{i-1}=R[-\Theta_{i}]$, the position of a particle $\left(\begin{smallmatrix}x\\ y\end{smallmatrix}\right)$ after time $t_{k}=\sum_{j=1}^{k}\tau_{j}$ is given by

[TABLE]

Numerical model

We use the analytical solution of the RPM flow in (9) to perform Lagrangian particle tracking on $100\,000$ particles that start at the source. As two particles with the exact same initial location follow the same trajectory, we sample the initial particle locations according to $(0,1)^{T}+\epsilon$ with $\epsilon\sim(\mathcal{N}(0,10^{-5}))^{2}$ to ensure that the particles are initially concentrated in the same tiny region inside the domain without having the same initial coordinates. Due to this definition, particles may be initialized outside the flow domain. In this case, we resample the initialization until the particle ends up inside the unit circle. We evolve these particles until $t_{\text{max}}=20$. Using the resulting particle distribution for different times $t$, we approximate the fraction of the domain covered by particles $\bar{M}$ by dividing the domain into cells of size $0.02^{2}$ and counting the number of cells that contain at least one particle. We use a time step of $\tau$ for the particle tracking to reduce the computational complexity. Hence, we compute $\bar{M}$ only for times $t\in\{n\tau:n\in\mathbb{N}\}$. We use linear interpolation to get an approximation of $\bar{M}$ for all other times. The choice of $t_{\text{max}}=20$ was found to be sufficient to get close to the steady state of $\bar{M}$ for most configurations (see supplementary material, Figs. S10, S11 and S12), hence we do not expect large changes of the sensitivity metrics for $t>t_{\text{max}}$.

2.1.2 The quadrupole flow, a high dimensional problem

Model setup

As second model, we consider the quadrupole flow introduced in [mays_2012]. The flow is produced by four wells located in a diamond shape around the center of the domain that represents a two-dimensional confined aquifer (see Fig. 2). The wells are operated according to the extraction and injection protocol of [mays_2012], where the four wells inject and extract water one at a time at different pumping rates. In the ideal case, with constant hydraulic conductivity, the sequence is designed such that the solute, which is initially placed in the center of the domain, is not extracted. In our experiments, extractions appear due to our modifications to the well locations and the hydraulic conductivity field. If solute gets extracted, it leaves the domain and will not be reinjected. Due to this solute extraction, the fraction of the domain covered by solute $\bar{M}$ is an improper metric to quantify mixing in this setup. While $\bar{M}$ mostly increases with time, it decreases whenever the system experiences solute extraction (see supplementary material, Fig. S1). The results with this metric would hence be highly affected by the solute extraction. We hence use the maximal solute concentration at the end of each of the $12$ stress periods as our quantity of interest, which is also known to be a good indicator for mixing efficiency [chiogna_2017]. We highlight that there are large differences between the RPM flow and the quadrupole flow. While the RPM flow follows a more synthetic approach, using a non-diffusive system and applying assumptions on the reinjection protocol, the quadrupole is more physical. A detailed discussion of these differences can be found in [lester_2013a].

Uncertain parameters

Our analysis reproduces the work of [ziliotto_2025]. The flow is parametrized by $16$ parameters. The first $8$ parameters describe the position of the four wells. We determine the exact well location by perturbing the ideal well locations with polar coordinates $(r,\theta)$ sampled uniformly from their respective intervals. We show this process in Fig. 2. Also, the pumping rates of the four wells underlie some uncertainty. We introduce a multiplier $q$ for each well that introduces a $\pm 5$ % uncertainty to the pumping rates defined in [mays_2012]. Finally, we use an uncertain hydraulic conductivity field that specifies the distribution of four materials in the aquifer while each material covers $25$ % of the domain area (see Fig. 2). The corresponding hydraulic conductivity values $k_{1}$, $k_{2}$, $k_{3}$, and $k_{4}$ are sampled according to the log uniform distribution. The exact sampling intervals for the $16$ uncertain parameters are given in Table 1.

Groundwater flow and transport equations

We solve the groundwater flow and transport problem using an Eulerian framework. For the two-dimensional groundwater flow in an isotropic confined aquifer, we use the equation [fetter_2001, bear_1988]

[TABLE]

with hydraulic head $h$ [$L$], storage coefficient $S$ [-], transmissivity tensor $T=Kb$ [$\frac{L^{2}}{T}$] with hydraulic conductivity tensor $K$ [$\frac{L}{T}$] and thickness of the aquifer $b$ [$L$] and source and sink term $W$ [$\frac{L}{T}$]. The transport of a conservative solute due to the flow is given in [fetter_2001, bear_1988] by

[TABLE]

with solute concentration $c$ [$\frac{M}{L^{3}}$] and hydrodynamic dispersion tensor $D$ [$\frac{L^{2}}{T}$]. The velocity vector $v=-\frac{K}{n_{e}}\nabla h$ [$\frac{L}{T}$] is given by Darcy’s law and depends on the effective porosity $n_{e}$ [-] and the solution of (10) [fetter_2001].

Numerical model

We simulate the engineered injection and extraction (EIE) system using MODFLOW-2005 [arlenw.harbaugh_2017], to solve the groundwater flow and MT3D-USGS [tonkin_2016], to solve the transport of a conservative solute. For the simulation of the wells, we use the well package of [arlenw.harbaugh_2017]. The model consists of a square domain with extension $201$ m $\times$ $201$ m and a grid size $\Delta x=\Delta y=0.5$ m. For the boundaries, we apply a constant head on the east and west sides and no flow on the north and south sides. We generate the random hydraulic conductivity field using the field generator [chiang_2005], with $\mu_{\log(K)}=0.7$, $\sigma_{\log(K)}^{2}=0.25$, and correlation length $\lambda=10$ m. For the values in the lowest quartile, we assign $k_{1}$, in the second-lowest quartile $k_{2}$, in the second-highest quartile $k_{3}$, and in the highest quartile $k_{4}$. The total duration of the extraction and injection sequence is $t_{\text{max}}=75$ days, which is divided into $12$ stress periods of $6.25$ days. We save the concentration distribution of the solute after each stress period. The initial normalized concentration of the plume $c_{0}(x)$ occupies an $11\text{ m}\times 11\text{ m}$ square in the center of the domain, and we use a hydrodynamic dispersion tensor $D=4\cdot 10^{-2}$ m2/day. A more detailed summary of the parameters of the quadrupole system is given in the supplementary material Table S1 as well as in [ziliotto_2025].

2.2 Sensitivity analysis

Sensitivity analysis is a powerful tool to understand the influence of input parameters on the output of a model. It can be used to identify the most important parameters, to reduce the dimensionality of the problem, and to improve the model’s performance. In this work, we use three different methods for sensitivity analysis. We estimate how sensitive the mixing efficiency of two chaotic flow systems introduced in Section 2.1.1 and Section 2.1.2 depends on the flow parametrization. Given a parametrization $x=(x_{1},...,x_{n})\in\mathbb{R}^{n}$, we compute the mixing efficiency using $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$. For the RPM flow, $f(x)$ is given by the area covered by particles at time $t$, and for the quadrupole flow, we use the peak concentration at the end of each stress period.

2.2.1 Sobol indices

The Sobol index is a variance-based approach [sobol_2001, saltelli_2010, sobol_1993], derived on the ANOVA representation [archer_1997, jansen_1999, sobol_2001, sobol_2003, sobol_2005], of $f$. Details are given in the appendix Appendix A. For the model $f:\mathbb{R}^{n}\Rightarrow\mathbb{R}$ we compute the total variance $V$ and the conditional variances $V_{i}$ and $V_{i,j}$, that describe the output variance due to a single or pairs of input parameters, by

[TABLE]

[saltelli_2010, sobol_1993, homma_1996]. Using this, we can define the first- and second-order Sobol indices as

[TABLE]

Also the computation of higher-order interactions via $S_{i_{1},...,i_{s}}:=\frac{V_{i_{1},...,i_{s}}}{V}$ is possible [saltelli_2002a, homma_1996, sobol_1993, sobol_2001]. It holds that

[TABLE]

We also use the total Sobol index with respect to the parameter $x_{i}$ introduced in [homma_1996, saltelli_2002]

[TABLE]

It measures the contribution of the parameter $x_{i}$ considering its first-order effect and all interactions with other parameters [pianosi_2016]. Following [sobol_2001, archer_1997], the total Sobol index can also be written in terms of the classical Sobol indices as

[TABLE]

For the numerical estimation of the Sobol indices, one can use Monte-Carlo sampling [sobol_2001, saltelli_2010, saltelli_2002, sobol_1993]. In general, estimating Sobol indices of higher order is computationally expensive [saltelli_2002]. Therefore, we restrict our analysis to first-order, second-order, and total Sobol indices. Our numerical experiments use the Python package SALib [herman_2017, iwanaga_2022]. We use Saltelli Sampling and create the sample matrices $A$ and $B\in\mathbb{R}^{N\times n}$, with $N$ number of samples and $n$ input dimension of the model. The $i$’th column of the matrices contains input samples for $x_{i}$. We define $AB\in\mathbb{R}^{2N\times n}$ as the concatenation of the two matrices $A$ and $B$. Further, we use the matrices $A_{B}^{(i)}$ where all columns are taken from $A$, but the $i$’th column is replaced by the $i$’th column of $B$. In the same way, we also create $B_{A}^{(i)}$. The implementation of the first-order Sobol index follows [saltelli_2010]. The estimators for $V_{i}$ and $V$ for computing $S_{i}$ in (12) are given by

[TABLE]

For the second-order Sobol index, we use the estimator given in [saltelli_2002]

[TABLE]

Finally, the estimator for the total Sobol index is given in [saltelli_2010, jansen_1999] as

[TABLE]

The estimators in (15) and (17) are variations of a standard mean estimator for which the rate of convergence is in $O(N^{-1/2})$ [iooss_2015]. We hence expect a convergence rate of $\frac{1}{2}$.

2.2.2 Morris method

The second sensitivity analysis method under consideration is the Morris method, introduced in [morris_1991, campolongo_2007]. For this method, the input parameters $x$ are sampled as a set of trajectories from the unit hypercube $[0,1]^{n}$. From this distribution, they are later transformed into their actual distribution. Each element of the trajectory is bound to a discrete grid of $p$ levels ($p$ even number), allowing for $x_{i}\in\{0,\frac{1}{p-1},\frac{2}{p-1},...,1\}$. Starting from a random sample $x\in\mathbb{R}^{n}$, the trajectory is created by changing one parameter at a time until all $n$ parameters have been changed once. Thereby, the step size for this change is given by $\Delta=\pm\frac{p}{2(p-1)}$, to ensure a symmetric treatment of the input parameters [morris_1991], and avoid focusing on local behaviors [pianosi_2016]. We use the sign in front of $\Delta$ to guarantee that the new sample remains inside the unit hypercube $[0,1]^{n}$. A trajectory is hence given by $\{x,x+\Delta e_{j},x+\Delta e_{j}+\Delta e_{k},...\}$, with $e_{j}$ being the unit vector with a one entry at position $j$.

From a trajectory, we compute the elementary effect of each input factor $x_{i}$ by

[TABLE]

The values of these elementary effects can be interpreted as local measures of sensitivity at point $x$ in direction $x_{i}$. Estimating global sensitivity requires approximating the first and second moments of the elementary effects, which we approximate by creating $N$ trajectories [morris_1991, campolongo_2007]. The first and second moments are given by

[TABLE]

where $d^{(j)}\in\mathbb{R}^{n}$ is the vector of the elementary effects $d_{i}^{(j)}$ that correspond to the $j$-th trajectory. To avoid the problem of non-monotonic models where elementary effects may have different signs, [campolongo_2007] introduces a refined version by estimating the mean over the absolute value of the elementary effects

[TABLE]

Again, we see that the sensitivity metric is in the form of a mean estimator, and hence, we expect a convergence rate of $\frac{1}{2}$ with respect to the number of samples $N$.

Comparing (19) with the formula for approximating the conditional variances $V_{i}$ and $V_{i,j}$ in (15) and (16), we observe that the quantities have different units. Given that $f(x)$ has the unit $[F]$, then $V_{i}$ and $V_{i,j}$ are given in terms of $[F^{2}]$ while $\mu_{i}^{*}$ has unit $[F]$. The sensitivity ranking remains unaffected by this; however, if one compares the actual values of the two metrics, this can have a significant influence. Therefore, whenever we consider the actual values of the Morris scores we use

[TABLE]

The outputs of this method are $n$ values of $\mu_{i}^{*}$, $\mu_{i}^{2*}$ and $\sigma_{i}$, each corresponding to one input parameter. Unlike the Sobol indices, where the indices sum up to one in (13), the magnitude of $\mu_{i}^{*}$, $\mu_{i}^{2*}$ and $\sigma_{i}$ depend on the overall variability of the model. As we aim to compare the sensitivity scores for different times in our numerical experiment, we normalize the values of $\mu_{i}^{*}$, $\mu_{i}^{2*}$, and $\sigma_{i}$ such that the metric gets quantitative [pianosi_2016]

[TABLE]

The normalized values can be interpreted as a percentage contribution to the sensitivity. A large value of $\hat{\mu}_{i}^{*}$ indicates a high sensitivity of the model output with respect to the input parameter $x_{i}$. If additionally $\hat{\sigma}_{i}$ is small, i.e., $\frac{\hat{\sigma}_{i}}{\hat{\mu}_{i}^{*}}\leq 0.1$, we can assume that the model output depends almost linearly on $x_{i}$. Higher values of $\hat{\sigma}_{i}$ indicate the presence of interactions and/or non-linearities in the model. For $\frac{\hat{\sigma}_{i}}{\hat{\mu}_{i}^{*}}\in[0.1,1]$, we assume monotonous or almost monotonous influence and for $\frac{\hat{\sigma}_{i}}{\hat{\mu}_{i}^{*}}\geq 1$ highly non-monotonous non-linearities or interactions [garciasanchez_2014, merchan-rivera_2022, richieri_2024, pianosi_2016].

2.2.3 Active subspace method

Although, the active subspace method (ASM) is mostly used for its ability to identify the most important directions in the input parameter space of the model and reduce a model’s dimensionality [constantine_2013, constantine_2015, parente_2019, parente_2019a], it can also be used to compute a global sensitivity metric, the activity scores [constantine_2017, erdal_2019, erdal_2020, bittner_2020]. The variation of the model within its input parameter space $[0,1]^{n}$ is given by the covariance matrix [constantine_2015]

[TABLE]

As $C$ is symmetric, positive-semidefinite, there exists an eigenvalue decomposition with an orthogonal matrix of eigenvectors $W$ and a diagonal matrix of eigenvalues in decreasing order $\Lambda=\text{diag}(\lambda_{1},...,\lambda_{n})$. An estimation of $C$ in (22) can be computed using Monte-Carlo sampling [constantine_2017]. Given $N$ samples $x_{j}\in[0,1]^{n}$, we set

[TABLE]

To lower the cost when computing multiple sensitivity measures simultaneously, we use the elementary effects of the Morris method and set $\nabla\hat{f}_{i}(x)=d_{i}(x)$, where $d_{i}(x)$ is the elementary effect of the Morris method in (18). For each eigenpair $(\lambda_{i},w_{i})$ of $C$, it holds that

[TABLE]

[constantine_2017]. We can hence interpret $\lambda_{i}$ as the importance of the direction $w_{i}$ in the input parameter space on the model output.

Starting from the eigenvalue decomposition of the covariance matrix, the authors of [constantine_2017] define the activity score by

[TABLE]

The value of $k$ is typically chosen to be the location of the biggest gap between the eigenvalues in log scale, i.e., $\log(\lambda_{k})\gg\log(\lambda_{k+1})$ [constantine_2015]. In the context of surrogate modeling, the first $k$ dimensions, corresponding to $w_{1},...,w_{k}$, are the active dimensions that impact the model output most. The remaining $n-k$ dimensions are the inactive dimensions [constantine_2017]. In cases without a significant eigenvalue gap, $k=n$ is reasonable. Although $\alpha_{i}(n)$ also includes the directions of less importance, these directions are weighted by a smaller eigenvalue and contribute less to the sensitivity value. During our experiments, we estimate the sensitivity of the input parameters at different times. Since we find a different location for the gap depending on the time step considered, we use $k=n$ for all times to ensure a consistent approach. It is shown in [constantine_2017] that, in this case, the activity score is equivalent to the derivative-based global sensitivity measure

[TABLE]

For applications where the surrogate model is not of interest, this is advantageous, as the evaluation does not require the eigendecomposition of the covariance matrix, but only the sum of the squared derivatives in (25). It also highlights that for our implementation with $\frac{\partial f(x_{j})}{\partial x_{i}}=d_{i}(x_{j})$, the activity score can be seen as a modification of $\mu_{i}^{*}$ from the Morris method, using the mean of the square elementary effects in place of the mean of the absolute elementary effects. Analyzing the units of $\alpha_{i}$ based on (24), we obtain $[F^{2}]$ alike for the Sobol indices. The values are hence comparable without adding a transformation.

As for the Morris method, we normalize the activity scores as

[TABLE]

to make the metric quantitative [pianosi_2016], ensuring comparability of the results at different times. We compute the activity scores from the eigenvalues of the covariance matrix in (23). As $\hat{C}$ converges to $C$ with order $\frac{1}{2}$, we expect a convergence rate of $\frac{1}{2}$ for the activity scores $\hat{\alpha}$.

The signs of the eigenvector entries corresponding to the largest eigenvalues are also interesting in analyzing the relationships between input and output. If the sign is the same, a simultaneous increase of both hyperparameters results in the largest output change; if the sign is different, we need to decrease the value of one hyperparameter and increase the other [constantine_2017].

When constructing a surrogate model, we need the first $k$ eigenvectors $W_{k}\in\mathbb{R}^{n\times k}$ of the covariance matrix. We can then reuse the samples $(x_{i},f(x_{i}))$ from the approximation of $\hat{C}$ and compute the projection $y_{i}=W_{k}^{T}x_{i}$. The surrogate model is given by a regression surface $G:\mathbb{R}^{k}\rightarrow\mathbb{R}$ of the pairs $(y_{i},f(x_{i}))$ with $G(y_{i})\approx f(x_{i})$ [constantine_2015]. We use a polynomial fit of order $m$ to construct $G$ using the Python package Scikit-learn [pedregosa_2011].

3 Results and Discussion

During our numerical experiments, we analyze the sensitivity of the two flow systems introduced in Section 2.1. We compare the different sensitivity metrics by taking the dimensionality of the problem into account. For the comparison, we use the total Sobol index $S_{T}$, the Morris index $\hat{\mu}^{2*}$, and the activity score $\hat{\alpha}$, as all three quantities represent the total-order influence of one hyperparameter on the model output. After demonstrating that the results for the RPM flow are consistent over all three methods, we look at the higher-dimensional quadrupole flow. Here, we only consider the Morris method and ASM to reduce the computational cost. We motivate this choice by the convergence behavior of the metrics for the RPM flow, where, for reaching the same relative accuracy, Morris scores require significantly fewer model evaluations compared to Sobol indices.

3.1 Sensitivity analysis on low dimensional problems

In our first experiment, we analyze the sensitivity of the two versions of the RPM flow. We regulate the non-randomized versions by varying the two hyperparameters $\Theta$ and $\tau$ and the randomized version with the hyperparameters $\bar{\Theta}$, $\Theta_{r}$, $\bar{\tau}$ and $\tau_{r}$. For the sampling intervals, we choose one configuration that allows for a large range of values for the flow parameters and one configuration with a more restricted choice of flow parameters, to investigate the effects of the prior information on the parameters. We give the exact intervals in Table 1. For each configuration, we compute all three sensitivity metrics defined in Section 2.2. To do so, we create Saltelli samples for the first, second and total-order Sobol indices (12), (17) and Morris samples for the Morris scores (19) and the activity scores (24), with $N=4096$. Our convergence analysis demonstrates that this sample size is large enough to ensure that all methods converge.

We assess the result of our sensitivity analysis using two approaches. We start by ranking the parameters according to their sensitivity scores. The ranks are often used to present the sensitivity analysis results in the literature [pianosi_2016, sarrazin_2016, ziliotto_, richieri_2024, mazziotta_2024] and allow for an easier comparison between the different methods. To better understand the sensitivity, we also look at the exact values of the sensitivity metrics for each parameter. Doing so, we use a transformation on the Morris score to obtain $\mu^{2*}$ as presented in (20) and normalize the Morris scores and activity scores using (21) and (26) to allow for inter-comparison between the different time steps.

3.1.1 Sensitivity analysis for the non-randomized RPM flow

Starting with the sensitivity analysis for the non-randomized RPM flow, we show the rankings in Fig. 3. For both experiments, with large and small hyperparameter intervals, $\tau$ seems to be the most sensitive hyperparameter for small times. Hence, at the beginning of the experiment, the number of performed rotations is more important to the result than the rotation angle. At time $t>1.5$, there is a switch, and $\Theta$ gets more sensitive. Thereby, the switch happens simultaneously for the Morris method and the ASM ($t\approx 2.5$ for the large interval experiment and $t\approx 4.75$ for the small interval experiment). Although the switch for the total Sobol index appears earlier ($t\approx 1.75$ for the large interval experiment and $t\approx 3.25$ for the small interval experiment), we conclude that the results over the three methods are very consistent.

To better understand the sensitivity, we extend our analysis and show the sensitivity metrics’ values in Fig. 4. As the curves for the total Sobol index $S_{T}$, the Morris index $\hat{\mu}^{2*}$, and the activity score $\hat{\alpha}$ show the same behavior, the analysis reinforces that the results for different sensitivity metrics reflect the same underlying sensitivity for both experiments. However, considering the exact values instead of the ranking reveals additional information. For the experiment with the large hyperparameter intervals, we find a minimum in the sensitivity of $\tau$ and a maximum in the sensitivity of $\Theta$ at $t\approx 5$. We notice the variable line for the sensitivity scores at small times for the small hyperparameter intervals experiment in Fig. 4(b). This behaviour can be explained by the fact that the overall output variability is very small at the beginning of the experiment, and hence although $\tau$ reaches higher sensitivity scores, none of the hyperparameters is dominant. Due to the normalization of the sensitivity metrics, we cannot yield zero sensitivity for all hyperparameters. We show the output variability in the supplementary material, Fig. S2b.

Using the three different methods, we can also analyze hyperparameter interactions. For the Sobol method, we compute the first- and second-order Sobol indices. As we consider a two-dimensional model where $A_{B}^{(1)}=B_{A}^{(2)}$ and $A_{B}^{(2)}=B_{A}^{(1)}$, this comes with no additional cost, and the results give a detailed impression about hyperparameter interactions. The first-order Sobol index shows the influence of one hyperparameter without any interactions, while the second-order index represents the interaction between exactly two hyperparameters. We show all first and second-order Sobol indices in Fig. 4. We find differences comparing the experiment with large and small hyperparameter intervals. For large hyperparameter intervals, there is a strong interaction between the two hyperparameters for $t<5$ that decreases over time. For the small hyperparameter intervals, after we exit the time frame where the output variability is too small to draw reasonable conclusions on sensitivity at $t\approx 4$, the interaction stays more or less constant throughout the experiment. Using the Morris method, we consider $\frac{{\hat{\sigma}}}{\hat{\mu}^{*}}$ to analyze the type of interactions inside the model. Figure 5 displays the relationship between $\hat{\mu}^{*}$ and $\hat{\sigma}$ over time for the two experiments. We see nonlinear interactions for the experiment with large hyperparameter intervals, while for the small hyperparameter intervals, the interactions are at the transition between nonlinear and almost monotonous. Finally, one can also observe relationships between hyperparameters using ASM by considering signs of the eigenvector entries corresponding to the largest eigenvalue [constantine_2017]. In the small hyperparameter intervals experiment, the first eigenvector has the same sign for both hyperparameters over the whole experiment, indicating a positive correlation. We cannot find any indication of a linear relationship for the large hyperparameter intervals, as most of the time, one of the eigenvector entries is close to zero. We give the eigenvectors in the supplementary material, Fig. S7.

An interpretation of the sensitivity analysis tells us that to enhance mixing in the RPM flow, controlling the rotation angle $\Theta$ is important, as long as we do not have to operate the system only for a short time. We highlight that an arbitrary rotation angle $\Theta$ might not be applicable for all applications. In groundwater remediation, for example, this may lead to a new pair of wells being drilled after each rotation of source and sink, which is practically infeasible, while a modification of $\tau$ could be implemented more easily. Interestingly, a more constrained hyperparameter interval does not allow for identifying a leading parameter controlling mixing at short times. For longer times, the size of the parameter interval does not influence the ranking of the parameters controlling mixing.

3.1.2 Sensitivity analysis for the randomized RPM flow

Looking at the sensitivity analysis of the randomized RPM flow, we start with the experiment with large hyperparameter intervals. The ranking in Fig. 6 shows the same trends for all three methods. After the large switch at $t\approx 1$, there is a second switch of the ranking between $\Theta_{r}$ and $\tau_{r}$. The timing for this switch is different over all three methods ($t\approx 16$ for the Sobol indices, $t\approx 12.25$ for Morris, and $t\approx 14$ for ASM). This inconsistency is due to the values of the sensitivity metrics, which are very close to each other as shown in Fig. 7. We observe that the sensitivities of $\Theta_{r}$ and $\tau_{r}$ converge towards similar values, which makes locating the exact time of the jump very difficult. The analysis of Fig. 7 also reveals that the metrics yield very similar values. This could only be achieved due to the transformation of the Morris scores in (20). The modification due to the transformation is not reflected in the sensitivity ranking, as the transformation is strictly monotonically increasing. Overall, we find that the deviations of the randomization $\Theta_{r}$ and $\tau_{r}$ are less sensitive than the means $\bar{\Theta}$ and $\bar{\tau}$. This result are due to the choice of the larger sampling intervals for the means than for the deviations. The sensitivity of the means is similar to the sensitivity of the respective hyperparameters in the non-randomized case, while the sensitivity of the deviations is nearly zero but increases with time. We explain these sensitivity increases of the deviation parameters as they enable particles to enter the KAM islands of the RPM flow. These KAM islands are non-mixing regions that arise in the non-randomized RPM flow; the randomization slowly deconstructs these islands such that the particles can cover the entire domain. As the effect of the mean values is much larger than the effect of the deviation parameters, it is important to characterize the mean of the random rotations in order to control mixing.

This time, analyzing interactions using Sobol indices comes with an additional cost of $N\cdot n$ samples to create $f(B_{A}^{(i)})|_{i=1}^{n}$. We only show the second-order Sobol index $S_{\bar{\Theta},\bar{\tau}}$ in Fig. 7, as all other second-order indices yield values very close to zero. There is, hence, only one significant interaction that appears between the two mean parameters. As for the non-randomized case, the interaction is strong for $t<5$, but decreases over time. Using the finding in (14) and $S_{T\bar{\tau}}\approx S_{\bar{\Theta},\bar{\tau}}$, we highlight that the total sensitivity of $\bar{\tau}$ is solely due to its interaction with $\bar{\Theta}$. As the total Sobol index is the sum of all lower-order Sobol indices involving the respective hyperparameter, we can conclude that we have nearly no higher-order interactions. For the Morris method, the analysis of Fig. 8(a) demonstrates that the interactions are again mostly nonlinear. Looking at the largest eigenvector of ASM, we find that most of the time, the two means and the two deviations are positively correlated. In contrast, we find a negative correlation between the deviations and the means. The eigenvector entries are given in the supplementary material, Fig. S8a.

The situation is similar for the small hyperparameter intervals. Although the ranking of the three methods in Fig. 9 seems less consistent, we find a pattern that all three methods follow upon close inspection. Most of the time, the deviation parameters $\Theta_{r}$ and $\tau_{r}$ are the least sensitive until the sensitivity of $\tau_{r}$ reaches rank one at the end of the experiment at $t\approx 18.5$. Most of the time, the mean values $\bar{\Theta}$ and $\bar{\tau}$ are the more sensitive parameters. Thereby, there is a time frame for $t$ between $4$ and $13$ where $\bar{\Theta}$ has rank one. For the remaining time, $\bar{\tau}$ is the most sensitive parameter. Due to the increase in the sensitivity of $\tau_{r}$ at the end of the experiment, the ranking of the mean values drops by one. In the randomized RPM flow, a more informed prior (i.e., a smaller range of the parameters) has a larger impact on the sensitivity analysis than the standard RPM flow. This means that randomly operating the system increases the importance of the duration of the strikes. Our claim that all three methods yield consistent results is further supported in Fig. 10, where the values of $S_{T}$, $\mu^{2*}$, and $\alpha$ follow a very similar path over the whole experiment.

We observe many switches in the ranking for small times in Fig. 9. This behavior is again due to the variable behavior of the sensitivity metrics at small times in Fig. 10, which is due to the low output variability of the model. We show the output variability in the supplementary material, Fig. S3b. Looking at later times in Fig. 10, we observe why the switches in the ranking again appear at slightly different times, as the sensitivity metrics for $\bar{\Theta}$, $\bar{\tau}$ and $\tau_{r}$ yield very similar values towards the end of the experiment.

Analyzing the parameter interactions with Sobol, we only find interactions between $\bar{\Theta}$ and $\bar{\tau}$. Therefore, we only show this index in Fig. 10. For $t<3$, this index reveals that the sensitivity of $\bar{\Theta}$ is mainly composed of the second-order interaction with $\bar{\tau}$. For $t>3$, the interaction decreases until the second-order index reaches a value close to zero, while the first-order index yields larger values. A comparison of the first-order and total Sobol indices indicates that the total sensitivity is mainly composed of the first-order influence of the parameters and less due to interactions. The analysis of $\frac{\hat{\sigma}_{i}}{\hat{\mu}_{i}^{*}}$ in Fig. 8(b) shows that the interactions between hyperparameters are between nonlinear and almost monotonous for $\bar{\Theta}$ and $\Theta_{r}$ and between almost monotonous and monotonous for $\bar{\tau}$ and $\tau_{r}$. Considering the first eigenvector of ASM presents the same results as for the large intervals experiment. The two means and the two deviations are positively correlated, while the means and deviations are negatively correlated. For the exact values of the eigenvector entries, we refer to the supplementary material, Fig. S8b.

From an application perspective, the results of the sensitivity analysis indicate that, depending on the time horizon available for the mixing process, different hyperparameters are relevant for generating optimal mixing. Therefore, when considering a short time horizon (e.g., $t_{\text{max}}=5$) the choice of $\Theta_{r}$ and $\tau_{r}$ does not affect the expected mixing enhancement and can therefore be disregarded. As mentioned for the non-randomized configurations, controlling $\bar{\Theta}$ and $\Theta_{r}$ might be infeasible for some applications, where the locations of source and sink can not be chosen arbitrarily; therefore, the modification of $\bar{\tau}$ and $\tau_{r}$ seems much more reasonable. Given that at least one of the two parameters is classified as sensitive throughout the whole experimental time, we can still control our results using these two parameters.

We analyze the convergence to better understand the quality of the sensitivity metrics. We perform this analysis on the large intervals experiment of the randomized RPM flow. To do so, we compute the metrics for mini-batches of size $N\leq 512$ and compare the results to those with $N=4096$. A relative error estimate $\delta(N,4096)$ approximates the mean relative error between the mini-batch results with batch size $N$ and the result with $4096$ samples. Figure 11 shows the relative error estimate for different mini-batch sizes. All sensitivity metrics converge with order $\frac{1}{2}$, as expected. The same rate was also found in [constantine_2017] for the activity scores and total Sobol indices. Be aware that using the same value of $N$ leads to $N(n+2)$ samples for the Sobol indices and $N(n+1)$ samples for the Morris method and ASM. Although using more samples, the error of the Sobol indices is larger than the error of the two other methods. We find the highest accuracy for Morris, which needs at most $\frac{1}{4}$ of the samples to yield the same relative error as the total Sobol index. The results highlight that the Morris method is most efficient for this higher-dimensional problem, where the analysis is limited by the cost of computing the required Monte-Carlo samples. This finding is also supported in the literature [pianosi_2016, yang_2011]. Additionally, we can reuse the samples to compute the activity scores of the ASM, giving us more insights into the model without additional computational costs. For a more detailed analysis, we show the convergence also for both randomized RPM flow experiments, including also the first- and second-order Sobol indices, in the supplementary material, Figs. S5 and S6.

As described in Section 2.2.3, ASM enables us to construct a surrogate model. We showcase this process using a polynomial fit of order $5$ and input dimensions $k\in\{1,2,3,4\}$. The surrogate is constructed solely using the Morris samples. We use the Saltelli samples to compute the mean absolute error of the surrogate in Fig. 12. The results show that two input dimensions are sufficient to approximate the model output well. We notice that the mean absolute error of the surrogate is maximal for times with a high output variability in Fig. S3. One can increase the accuracy further by fine-tuning the polynomial order of the surrogate.

3.2 Sensitivity analysis on high dimensional problems

In our next experiment, we analyze the sensitivity of the quadrupole model from Section 2.1.2. This model is significantly higher-dimensional than the previous flow systems, having $16$ hyperparameters. Due to the higher computational cost of the Monte-Carlo sampling, we restrict the analysis and only compute the Morris sensitivity and ASM activity scores. Our results on the RPM flow show that these sensitivity metrics yield consistent results with the Sobol indices while requiring fewer samples to converge. For our analysis, the choice of $N=10$ is sufficient for the screening of the sensitive parameters. A detailed analysis of the results of the sensitivity analysis with Morris is given in [ziliotto_2025]; we therefore concentrate more on the comparison of the Morris and the ASM method. Further, due to the complexity of the model and the high dimensional input, a surrogate model as presented in Section 2.2.3 cannot be applied without additional sampling.

We start our analysis by showing the results for the sensitivity analysis of Morris after the last stress period using $N=10$ trajectories in Fig. 13. Despite the low number of samples, we can draw some interesting conclusions from these results. While it is generally enough for practical applications related to model calibration and system design to perform a sensitivity screening, i.e., identification of the group of most sensitive parameters, we can observe that the ranking is not fully converged yet, and small variations could be expected. Despite the exact values can vary significantly within the respective confidence intervals, we should consider the interaction among the parameters and hence an overlap of the confidence intervals of two parameters, does not necessarily mean that the two parameters could switch their ranking. Conclusions based on the exact values of the sensitivity index on the contrary need to be considered with special care. Such analysis is often not needed in practical applications focusing on inverse modelling and optimization, however it may be valuable in order to get insights about how dominant the sensitive hyperparameters are compared to the less sensitive ones. In the work of [sarrazin_2016], the authors use the size of the $95$ % confidence interval as a convergence criterion for the sensitivity analysis. They state that the method is converged if the largest confidence interval has a size below $5$. Knowing that Morris converges with a rate of $\frac{1}{2}$, we expect to get all confidence intervals below size five with $N=48$ sample trajectories. So, the convergence for the sensitivity indices comes at a significant computational cost, that may not be justified by the purpose of the study.

We show the complete results of the sensitivity ranking of Morris and ASM in Fig. 14(a). The cells with no values belong to wells activated at a later stress period, like the north well that starts pumping in stress period seven. The ranking of the two methods is very similar. We observe the same trends, like the high sensitivity of the conductivity values $k_{1}$, $k_{2}$, $k_{3}$, and $k_{4}$, followed by the parameters related to the east well location $r_{e}$ and $\theta_{e}$. We only find minor ranking differences between the two methods. The maximal deviation is three for the hyperparameter $\theta_{e}$ and stress period $7$. Other than that, we have seven deviations of two and $41$ of one. The remaining $101$ instances match exactly.

Next, we look at the exact values of the sensitivity metrics in Fig. 14(b). As discussed previously, the values require more samples to reach convergence, and hence the results need to be observed more carefully. Additionally, we highlight that the two methods are unlikely to obtain the same values, as they use the same data but different formulas to compute the sensitivity metrics. Looking at Fig. 14(b), we find that both methods (after using the transformation in the Morris scores in (20)) yield very similar sensitivity scores. Overall, the highest sensitivity is concentrated on the four conductivity values, followed by the hyperparameters corresponding to the east well and the discharge of the west well. The remaining hyperparameters are non-sensitive.

The comparison between Morris and ASM is useful to confirm the results presented by [ziliotto_2025]. In fact, it is good practice to compare different approaches in sensitivity analysis, since there can be a difference in the conclusion. From a physical point of view, we conclude that in the quadrupole experiment, the values of the hydraulic conductivities are the most important parameters, followed by the operational parameters of the wells activated at the beginning of the process. While the efficiency depends on both the design parameters of the system (configuration of the pumping wells) and the physical parameters of the aquifer (hydraulic conductivity), we can only vary the design parameters, e.g., pumping rates and well locations. The control of the hydraulic conductivity on the system is expected and was also identified in previous studies [bertran_2023]. On the contrary, the dependence on the design parameters is very relevant, although they are not the most sensitive ones, since this information can be used in the optimization of the design of the system.

4 Conclusion and Outlook

In our work, we use three different sensitivity analysis methods, Sobol indices, Morris scores, and activity scores of the active subspace method (ASM), to analyze the mixing enhancement of two chaotic flow fields depending on their design parameters. We use a time-varying approach to capture the temporal evolution of the sensitivity. To do so, we introduce a scaling to normalize the sensitivity metrics over time, making all sensitivity metrics quantitative. For the comparison of the values of the sensitivity metrics, we consider their units and introduce a transformation to make the units of the Morris scores the same as those of the other two metrics. Our analysis supports the decision-making on the best method for sensitivity analysis, depending on the input dimension, the computational budget, and the goals of the analysis. We show that the sensitivities are largely consistent over the three methods while having different computational costs and giving different insights into interactions between hyperparameters. With Sobol indices, we can directly compute second-order interactions at the cost of $N\cdot n$ additional samples. The Morris method can indicate the existence of interactions by looking at the ratio $\frac{\hat{\sigma}}{\hat{\mu}^{*}}$, and the first eigenvector of ASM gives insights into the important linear directions of the input hyperparameter space. Although the second-order Sobol indices come with increased computational costs, the effort might be reasonable, especially for low-dimensional applications, as the results give a precise and easy-to-understand estimation of second-order interactions. In contrast, when using high-dimensional or computationally expensive methods, we demonstrate the strength of the Morris method and, subsequently, of the ASM activity scores. Using our low-dimensional experiment, we demonstrate that Morris yields the same relative error as the total Sobol index with only $\frac{1}{4}$ of the model evaluations. Additionally, as we reuse the Morris samples to compute the ASM activity scores, we can analyze the model’s sensitivity with two approaches at the cost of one approach. Beyond sensitivity analysis, the ASM method offers the possibility to develop a surrogate model that approximates the original model’s behavior at a low cost. These surrogates can be used for a broad range of tasks, such as uncertainty quantification, optimization, or real-time decision-making, causing ASM to be a good starting tool for deeper analysis.

The results of our sensitivity analysis show that mixing in chaotic systems relies sensitively on their flow parametrization. The sensitivity rankings can be used to optimize the most influential design parameters of a chaotic flow field to enhance mixing. Additionally, when implementing the flow field, the more influential hyperparameters must be implemented precisely, as low aberration can lead to significant differences in the model response. The time-varying analysis reveals that the sensitivity ranking changes over time. For applications that generate effective chaotic mixing, this means that, depending on the duration of the process, different hyperparameters need special consideration. For the randomized RPM flow, for example, the sensitivity of the deviation parameters $\Theta_{r}$ and $\tau_{r}$ increases over time. The longer the mixing lasts, the more important the deviation parameters are to the final mixing enhancement, while for short times, they play a minor role. Additionally, we find that the sensitivity also depends on the hyperparameter intervals. To generate meaningful results, we recommend using hyperparameter intervals of comparable sizes. We emphasize that the mixing enhancement can also depend on physical parameters that can not be controlled to enhance mixing. If a physical parameter is highly sensitive, however, it may be beneficial to estimate its value through measurements before optimizing the design parameters to enhance mixing.

Especially for the quadrupole flow, the computational cost of the flow and transport model limits our analysis in terms of model realizations. In general, the sensitivity values require more samples for convergence than a ranking or screening of the hyperparameters [sarrazin_2016]. The authors of [campolongo_1997] show that $10$ trajectories are insufficient for their $35$-dimensional model. Other authors use more trajectories like [richieri_2024] with $500$ trajectories on a $17$-dimensional model. We reinforce the hypothesis that the optimal number of trajectories is model-specific and hard to predict. Regarding the ASM activity scores, we extend the original method to enable the usage of Morris samples. Consequently, the theorems presented in [constantine_2017] do not hold anymore. However, the connection between our activity scores and the Morris methods motivates that we yield valid results. If $m=n$, like in our experiments, the Morris and ASM activity scores are based on the same data and use a similar formula, which means that the results of these two methods will always be correlated.

For future work, we envision comparing our adapted version of the activity scores to the original activity scores. As a robust sensitivity analysis often requires multiple methods, a good agreement of our adapted activity scores and the original version would support our claim that we can get this sensitivity method “for free” when using Morris. Further, our analysis only features one method that can explicitly indicate second-order interaction. Including another method to analyze these effects and comparing the results to our second-order Sobol indices would be interesting. The approach in [campolongo_1999, cropp_2002] is an extension to Morris that requires a total of $N\cdot(\frac{n^{2}+n}{2}+1)$ model evaluations to compute these interactions (including $\hat{\mu}^{*}$ and $\hat{\sigma}$). As this number grows quadratically with the input dimension $n$, it becomes infeasible for most high-dimensional models. We would also like to investigate the influence of randomization on the mixing enhancement of the RPM flow. Noticing that our non-randomized RPM flow is equivalent to our randomized RPM flow with $\Theta_{r}=0$ and $\tau_{r}=0$, we found the sensitivity of the mean parameters $\bar{\Theta}$ and $\bar{\tau}$ to be comparable to $\Theta$ and $\tau$ in the non-randomized case. At late times, the sensitivity of the deviation parameters $\Theta_{r}$ and $\tau_{r}$ increases. The higher sensitivity of the deviation parameters indicates that the randomization influences the mixing enhancement and could be an important factor in optimizing chaotic flow fields to enhance mixing. The mixing efficiency depends on design parameters, like well locations and pumping rates. Hydrological parameters, like the hydraulic conductivity, are fixed and are, in general, subject to a large uncertainty. By applying risk-aware optimization, it may be possible to optimize the design parameters by considering the uncertainty of the hydrological parameters.

Statements & Declarations

\bmhead

Funding

This research is a result of the ChaosAD project (Chaotic ADvection and Mixing Enhancement in Porous Media: The Quest for Experimental Evidence), which is supported by the Deutsche Forschungsgemeinschaft (DFG) under the grant numbers CH 981/8-1 and SI 2853/2-1. Financial support for Francesca Ziliotto was provided by TUM International Graduate School of Science and Engineering (IGSSE), within the INE2 project (Innovative Engineering Injection Extraction systems for in-situ groundwater remediation. From model- and laboratory-based evidence to stakeholder involvement). The authors gratefully acknowledge the scientific support and HPC resources provided by the Erlangen National High Performance Computing Center (NHR@FAU) of the Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU). The hardware is funded by the German Research Foundation (DFG).

\bmhead

Competing Interests

The authors have no relevant financial or non-financial interests to disclose.

\bmhead

Author Contributions

Conceptualization was performed by Carla Feistner and Gabriele Chiogna. Carla Feistner also worked on the formal analysis and methodology. Software, investigation, and visualization of the RPM flow were performed by Carla Feistner. For the Quadrupole flow, the implementation of the flow and transport model and the sensitivity analysis with Morris was done by Francesca Ziliotto, and the analysis with the active subspace method was added by Carla Feistner. The supervision was done by Gabriele Chiogna and Barbara Wohlmuth. The original draft was written by Carla Feistner, and all authors contributed to the review and editing of the manuscript. All authors read and approved the final manuscript. The project administration and funding acquisition were done by Gabriele Chiogna. We thank Prof. David Mays for the fruitful discussion on the manuscript.

\bmhead

Supplementary information

We created a PDF file that contains supplementary figures and tables that support the results presented in this manuscript. The file contains a detailed table describing the parametrization of the quadrupole flow, a figure motivating the choice of our quantity of interest to describe mixing, an analysis of the output variability of both chaotic flow fields, all convergence figures for the different experiments of the randomized RPM flow, and the eigenvectors of the covariance matrices for all experiments.