The effects of stress perturbations on friction are crucial for understanding earthquake triggering. Previous experimental studies have primarily been conducted at room temperature, where fault gouge materials typically exhibit velocity-strengthening and frictionally stable behaviour. In this study, we investigate how variations in effective normal stress () influence fault (in-)stability by performing -perturbation experiments on simulated carbonate fault gouges under fluid-drained, hydrothermal conditions. Our results indicate that in the velocity-neutral or -weakening regime, perturbing can reinforce frictional instability, leading to accelerated slow slips or enhanced stick-slip events. This effect is particularly pronounced when the excitation period () approaches or exceeds the characteristic recurrence period () associated with pre-perturbation instabilities. Stress drops of the resulting events can have larger amplitudes than expected from quasi-steady-state. When cyclic perturbations are imposed, slip events tend to synchronize at specific phases when is close to - notably between 0.5 pi and 1 pi in radian, corresponding to maximum destressing rate and minimum , respectively. Additionally, short-period () perturbations can induce significant shear stress reduction (or fault weakening), with magnitudes comparable to the stress drops from stick slips, yet they are surprisingly associated with acoustically quiet slow slip, suggesting a stabilizing effect. These findings underscore the critical role of perturbation period in controlling fault response. In the context of induced seismicity, our results imply that cyclic or monotonic fluid injections should be carefully designed, considering both perturbation amplitude and period. Properly turned cyclic injections could potentially mitigate seismic risk by promoting quiet, slow slip over seismic fault slip.

Understanding fault slip nucleation within the reservoir interval and its propagation beyond the reservoir is essential. Analytical and numerical studies have shown that, depending on the type of operation (injection/depletion), fault slip can nucleate at external or inner corners along the displaced fault system, driven by positive peak shear stresses. In the case of depletion, slip patches gradually start at the inner corners and grow towards the inner part of the reservoir, merging with further depletion. Conversely, injection or increased pore pressure leads to slip patches at external corners, potentially propagating beyond the reservoir into the overburden and underburden. We conducted triaxial experiments on small-scale (mm scale) cylindrical samples containing an entirely displaced vertical fault to investigate fault reactivation and slip nucleation in such settings. Two types of stress paths, monotonic and cyclic, were applied to examine the effects of stress patterns on slip nucleation. For this purpose, we utilized strain gauges to measure differential compaction along the displaced fault directly on the small-scale samples. Direct measurements with a strain gauge network adjacent to the displaced fault system during the monotonic test revealed that differential compaction intensifies from the top of the sample towards the internal corner at the center of the fault where different layers are juxtaposed vertically, indicating a variation in the stress field surrounding the fault plane. Furthermore, results from the cyclic test showed that the differential compaction increases with an increasing number of cycles. Our direct measurements near the displaced fault plane confirm/match the anomalies and peaks in stress observed in previous numerical and analytical studies.

Abstract Pore pressure fluctuation in subsurface reservoirs and its resulting mechanical response can cause fault reactivation. Numerical simulation of such induced seismicity is important to develop reliable seismic hazard and risk assessments. However, modeling of fault reactivation is quite challenging, especially in the case of displaced faults, i.e., faults with non-zero offset. In this paper, we perform a systematic benchmarking study to validate two recently developed numerical methods for fault slip simulation. Reference solutions are based on a semi-analytical approach that makes use of inclusion theory and Cauchy-type singular integral equations. The two numerical methods both use finite volume discretizations, but they employ different approaches to represent faults. One of them uses a conformal discrete fault model (DFM) while the other employs an embedded (non-conformal) fault model. The semi-analytical test cases cover a vertical frictionless fault, and inclined displaced faults with constant friction and slip-weakening friction. It was found that both numerical methods accurately represent pre-slip stress fields caused by pore pressure changes. Moreover, they also successfully cope with a vertical frictionless fault. However, for the case with an inclined displaced fault with a constant friction coefficient, the embedded method can not converge for the post-slip phase, whereas the DFM successfully coped with both constant and slip-weakening friction coefficients. In its current implementation, the DFM is therefore the model of choice when accurate simulation of local faulted systems is required.

We present expressions to compute the inverse of a Cauchy-type singular integral equation representing the relation between a double-peaked Coulomb stress in a fault or fracture and the resulting slip gradient in two distinct collinear slip patches. In particular we consider a situation where the patches are close enough to account for the influence of the slip gradient in one patch on the slip-induced shear stress in the other patch and vice versa. This situation can occur during depletion-induced or injection-induced fault slip in subsurface reservoirs for, e.g., natural gas production, hydrogen or CO2 storage, or geothermal operations. The theory for a single slip patch is well-developed but the situation is less clear for a configuration with two patches although the monographs of Muskhelishvili (1953) and Weertman (1996) provide earlier results. We show that the general inverse solution for the coupled two-patch problem requires six auxiliary conditions to ensure six physical requirements: boundedness of the slip gradient at the four end points of the slip patches and vanishing of the integrals of the slip gradient over the patches. Mathematically, the presence of two additional conditions, as compared to earlier formulations, corresponds to two undetermined coefficients in the general solution of the governing integral equation. Numerical simulation confirms that at least one of these is always non-zero in the coupled situation. For a coupled double-patch case with a symmetric pre-slip Coulomb stress pattern, the general inverse solution requires three auxiliary conditions. Moreover the conditions for the asymmetric case may be reduced to a set of four again, but these are different from the sets of four obtained earlier by Muskhelishvili (1953) and Weertman (1996). We illustrate the theory with a numerical example in which the evaluation of the Cauchy integrals is performed with a modified version of augmented Gauss–Chebyshev quadrature that relies on analytical inversion.

We critically review the derivation of closed-form analytical expressions for elastic displacements, strains, and stresses inside a subsurface reservoir undergoing pore pressure changes using inclusion theory. Although developed decades ago, inclusion theory has been used recently by various authors to obtain fast estimates of depletion-induced and injection-induced fault stresses in relation to induced seismicity. We therefore briefly address the current geomechanical relevance of this method, and provide a numerical example to demonstrate its use to compute induced fault stresses. However, the main goal of our paper is to correct some erroneous assumptions that were made in earlier publications. While the final expressions for the poroelastic stresses in these publications were correct, their derivation contained conceptual mistakes due to the mathematical subtleties that arise because of singularities in the Green's functions. The aim of our paper is therefore to present the correct derivation of expressions for the strains and stresses inside an inclusion and to clarify some of the results of the aforementioned studies. Furthermore, we present two conditions that the strain field must satisfy, which can be used to verify the analytical expressions. Previous authors derived mathematical expressions to aid in understanding the nature of earthquakes as a result of the production of fluids from the subsurface or the injection of fluids into it. These expressions represent relationships between stresses and pore pressures in subsurface reservoirs (i.e., sealed bodies of porous rock containing fluids). Although the expressions are correct, the reasoning behind their derivation contained flaws related to mathematical subtleties in the underlying theory. In the present paper we explain those errors and present the correct derivation. Moreover we present two additional results that can be of help in the derivation of similar expressions for different configurations or purposes, and we give an example to demonstrate the use of these expressions in computing stresses caused by fluid production from a faulted reservoir. We clarify and correct expressions for poroelastic strains and stresses obtained with inclusion theory in previous studies Green's functions for strains and stresses can be integrated to obtain analytical expressions outside, but not inside poroelastic reservoirs We provide conditions to verify analytical solutions for poroelastic strains based on physical and mathematical requirements

Recent laboratory and field studies suggest that temporal variations in injection patterns (e.g., cyclic injection) might trigger less seismicity than constant monotonic injection. This study presents results from uniaxial compressive experiments performed on Red Felser sandstone samples providing new information on the effect of stress pattern and rate on seismicity evolution. Red Felser sandstone samples were subjected to three stress patterns: cyclic recursive, cyclic progressive (CP), and monotonic stress. Three different stress rates (displacement controlled) were also applied: low, medium, and high rates of 10 −4 mm/s, 5 × 10 −4 mm/s, and 5 × 10 −3 mm/s, respectively. Acoustic emission (AE) waveforms were recorded throughout the experiments using 11 AE transducers placed around the sample. Microseismicity analysis shows that (i) Cyclic stress patterns and especially cyclic progressive ones are characterized by a high number of AE events and lower maximum AE amplitude, (ii) among the three different stress patterns, the largest b-value (slope of the log frequency-magnitude distribution) resulted from the cyclic progressive (CP) stress pattern, (iii) by reducing the stress rate, the maximum AE energy and final mechanical strength both decrease significantly. In addition, stress rate remarkably affects the detailed AE signature of the events classified by the distribution of events in the average frequency (AF)—rise angle (RA) space. High stress rates increase the number of events with low AF and high RA signatures. Considering all elements of the AE analysis, it can be concluded that applying cyclic stress patterns in combination with low-stress rates may potentially lead to a more favourable induced seismicity effect in subsurface-related injection operations.

Pore pressure changes due to fluid injection or withdrawal alter the rock stresses, which may potentially induce seismic events. Activities that have been associated with induced seismicity include geothermal energy production, subsurface gas storage, and natural gas production. Physics-based models are required to gain insight in the processes that lead to induced seismicity. When computing induced stresses with such models, a commonly made simplification is the assumption of uniform pressure changes across the entire reservoir. In reality, pressure gradients arise due to fluid production or injection. Here, we assess the effect of non-uniform pressure fields under steady-state flow conditions on induced stresses. We employ (semi-)analytical techniques to compute the corresponding pressure field and fault stresses. We are particularly interested in reservoirs with displaced faults (i.e., cases with nonzero fault offset), as shear stresses tend to concentrate at the reservoir corners along the faults. The stress profile along the fault becomes asymmetric under steady-state flow. The effect of fluid flow on fault stresses is larger in case of injection than in case of depletion. Injection with up-dip flow results in increased zones of fault slip near the bottom of the reservoir, while injection with down-dip flow results in increased slip near the top of the reservoir. The significance of steady-state flow on induced stresses can be estimated from the ratio of the average pressure change along the fault and the applied pressure gradient. The effect of steady-state flow is most relevant at the start of production or injection and diminishes with time. Thus, the effect of steady-state flow is only expected to be relevant when initially critically stressed faults are present. For non-critically stressed faults, the assumption of uniform depletion or injection is expected to be reasonable. An order-of-magnitude estimate of the effect of steady-state flow across displaced faults in the Groningen natural gas reservoir shows that the effect on fault stresses is probably negligible. A similar estimate of the effect in typical low-enthalpy geothermal doublets indicates that steady-state flow may possibly play a small role, in particular close to the injector. Nevertheless, site-specific assessments are necessary to quantify the effect in greater detail.

ABSTRACT Quantification of the poromechanical response of subsurface formations due to human-induced pore pressure fluctuations is critical for the performance and stability assessment of many geo-energy systems. In particular, natural faults in the subsurface introduce the hazard of induced seismicity. Numerical modeling of fault reactivation is challenging, while the specific details of induced stresses and fault slip in reservoirs with displaced (i.e. non-zero offset) faults may cause additional challenges depending on the type of numerical formulation employed. To facilitate the systematic development and testing of numerical tools for the simulation of induced seismicity in faulted reservoirs we developed a set of semi-analytical test problems of increasing complexity, based on inclusion theory and Cauchy singular integral equations. With these we investigate the accuracy of two recently developed Finite Volume (FV) schemes with collocated and staggered arrangements of unknowns. One of them employs a conformal discrete fault model (DFM) which can guarantee sufficient accuracy at the cost of adaptive mesh refinement but may suffer from modelling and computational challenges when addressing large-scale realistic geological configurations. The second one employs an embedded (or non-conformal) discrete fault model (EDFM) which avoids the need for excessive mesh refinement, but of which the accuracy and the range of applicability are still to be investigated. We found that both numerical schemes accurately represent the pre-slip Coulomb stresses, but show different degrees of accuracy in representing the resulting depletion-induced fault slip. The semi-analytical benchmark data are available via DOI 10.4121/22240309. INTRODUCTION The kernel of this paper is formed by a series of semi-analytical poro-mechanical test problems of increasing complexity with the aim to systematically compare the capacities of two poro-mechanical finite-volume-based simulation codes: one developed by Novikov et al. (2022b) which employs a discrete fault model (DFM). It forms part of a comprehensive porous media simulation package, the Delft Advanced Reservoir Terra Simulator (DARTS) and will be referred to with that acronym. The second code, developed by Shokrollahzadeh Behbahani et al. (2022), is based on a smoothed version of the embedded discrete fault model (sEFVM) and will be referred to with that last acronym. Both codes are being developed as part of the DeepNL Science4Steer project (NWO, 2017), and Appendices A and B give a brief overview of their characteristic features. The grids used in this study are presented in Appendix C.

We consider steady-state single-phase confined flow through a subsurface porous layer containing a displaced, fully conductive fault causing a sudden jump in the flow path, and we employ (semi-)analytical techniques to compute the corresponding pressures and fault stresses. In particular, we obtain a new solution for the pressure field with the aid of conformal mapping and a Schwarz–Christoffel transformation. Moreover, we use an existing technique to compute the poro-elastic stress field with the aid of inclusion theory. The additional resistance to fluid flow provided by a displaced fault, relative to the resistance in a layer without a fault, is a function of dip angle, fault throw divided by reservoir height, and reservoir width divided by reservoir height. Fluid flow has a larger effect on fault stresses in case of injection than in case of depletion, where injection with up-dip flow results in increased zones of fault slip near the bottom of the reservoir. Opposedly, injection with down-dip flow results in increased slip near the top of the reservoir. An order-of-magnitude estimate of the effect of steady-state flow across displaced faults in the Groningen natural gas reservoir shows that the effect on fault stresses is probably negligible. A similar estimate of the effect in low-enthalpy geothermal doublets indicates that steady-state flow may possibly play a small role, in particular close to the injector, but site-specific assessments will be necessary to quantify the effect.

We address aseismic fault slip and the onset of seismicity resulting from depletion-induced or injection-induced stresses in reservoirs with pre-existing vertical or inclined faults. Building on classic results, we discuss semi-analytical modelling techniques for fault slip including dislocation theory, Cauchy-type singular integral equations and the use of Chebyshev polynomials for their solution and an eigenvalue-based stability analysis. We consider slip patch development during depletion for faults with zero, constant static and slip-weakening friction, and our results confirm earlier findings based on numerical simulation, in particular the aseismic growth of two slip patches that may subsequently merge and/or become unstable resulting in nucleation of seismic slip. New findings include improved approximate expressions for the induced seismic moment per unit strike length and a description of the effect of coupling between the slip patches which affects both forward simulation and eigenvalue computation for high values of the ratio of fault throw to reservoir height. Our implementation based on analytical inversion and semi-analytical integration with Chebyshev polynomials is more efficient and more robust than our numerical integration approach. It is not yet well suited for Monte Carlo simulation, which typically requires sub-second simulation times, but with some further development that option seems to be within reach. Moreover, our results offer a possibility for embedded fault modelling in large-scale numerical simulation tools.

We present a scalable collocated Finite Volume Method (FVM) to simulate induced seismicity as a result of pore pressure changes. A discrete system is obtained based on a fully-implicit fully-coupled description of flow, elastic deformation, and contact mechanics at fault surfaces on a flexible unstructured mesh. The cell-centered collocated scheme leads to a convenient integration of the different physical equations, as the unknowns share the same discrete locations on the mesh. Additionally, a generic multi-point flux approximation is formulated to treat heterogeneity, anisotropy, and cross-derivative terms for both flow and mechanics equations. The resulting system, though flexible and accurate, can lead to excessive computational costs for field-relevant applications. To resolve this limitation, a scalable processing algorithm is developed and presented. Several proof-of-concept numerical tests, including benchmark studies with analytical solutions, are investigated. It is found that the presented method is indeed accurate and efficient; and provides a promising framework for accurate and efficient simulation of induced seismicity in various geoscientific applications.

A smoothed embedded finite-volume modeling (sEFVM) method is presented for faulted and fractured heterogeneous poroelastic media. The method casts a fully coupled strategy to treat the coupling between fault slip mechanics, deformation mechanics, and fluid flow equations. This ensures the stability and consistency of the simulation results, especially, as the fault slip is implicitly found through an iterative prediction-correction procedure. The computational grid is generated independently for embedded faults and rock matrix. The efficiency is further enhanced by extending the finite-volume discrete space by introducing only one degree of freedom per fault element. The embedded approach can lead to an oscillatory stress field at the fault, which damages the robustness of the implicit slip detection strategy. To resolve this challenge, a smoothed embedded strategy is devised, in which the stress and slip profiles are post processed within the iterative loops by fitting the best curve based on a least-square error criterion. The sEFVM provides locally conservative mass flux and stress fields, on staggered grid. Its performance is further investigated for several proof-of-the-concept test cases, including a multiple fault system in a heterogeneous domain. Results indicate that the method develops a promising approach for field-scale relevant simulation of induced seismicity.

We explore and develop a Proper Orthogonal Decomposition (POD)-based deflation method for the solution of ill-conditioned linear systems, appearing in simulations of two-phase flow through highly heterogeneous porous media. We accelerate the convergence of a Preconditioned Conjugate Gradient (PCG) method achieving speed-ups of factors up to five. The up-front extra computational cost of the proposed method depends on the number of deflation vectors. The POD-based deflation method is tested for a particular problem and linear solver; nevertheless, it can be applied to various transient problems, and combined with multiple solvers, e.g., Krylov subspace and multigrid methods.

The Gutenberg-Richter law describes the frequency-magnitude distribution of seismic events where its slope, the 'b-value', is commonly used to describe the relative occurrence of large and small events. Statistically significant b-value variations have been measured in laboratory experiments, mines, and various tectonic regimes (Wiemer & Wyss, 2002). An inversely proportional dependency of the b-value on the differential stress has been observed across different scales (Amitrano, 2003; Schorlemmer et al., 2005). Layland-Bachmann et al. (2012) have shown that this could explain the observed pattern of induced seismicity spatial-temporal b-value variations in Enhanced Geothermal Systems. In our study, we look for a similar relation applied to the Groningen gas field in the Netherlands.It is well known that the poroelastic changes in differential stress during gas extraction are influenced by the offset of the reservoir layer across the fault. Recently, Jansen et al. (2019) and Lehner (2019) proposed an analytical solution for stress changes on offset faults due to reservoir depletion. In a parallel study, we extended this solution to include the development of aseismic slip under slip weakening and the derivation of the onset of seismic slip.We utilize this formulation to derive the onset of seismic slip on theoretical faults of variable fault offset, dip, and reservoir thickness. Subsequently, we map our theoretical faults onto the pre-existing faults in the Groningen gas field, deriving fault segment-specific depletion levels at which the segment would become seismically active. We then simulate reservoir depletion conditions over time and assign an event magnitude to fault segments that move past their seismic activation depletion. To assign a magnitude, we use the observation that b-values are inversely proportional to differential stress, which is governed by the pore pressure depletion. Hence, we assume a simple inverse linear relation with pore pressure depletion. Each event magnitude is then randomly drawn from the probability density function of the Gutenberg-Richter distribution with the b-value assigned.We aim to compare the obtained catalogue and its b-value distribution both in time and space to the observed event-size distribution of the Groningen gas field as derived by Muntendam-Bos and Güdük (EGU abstract 2021).

We develop a collocated Finite Volume Method (FVM) to study induced seismicity as a result of pore pressure fluctuations. A discrete system is obtained based on a fully-implicit coupled description of flow, elastic deformation, and contact mechanics at fault surfaces on a fully unstructured mesh. The cell-centered collocated scheme leads to convenient integration of the different physical equations, as the unknowns share the same discrete locations on the mesh. Additionally, a multi-point flux approximation is formulated in a general procedure to treat heterogeneity, anisotropy, and cross-derivative terms for both flow and mechanics equations. The resulting system, though flexible and accurate, can lead to excessive computational costs for field-relevant applications. To resolve this limitation, a scalable parallel solution algorithm is developed and presented. Several proof-of-concept numerical tests, including benchmark studies with analytical solutions, are investigated. It is found that the presented method is indeed accurate, stable and efficient; and as such promising for accurate and efficient simulation of induced seismicity.

Reactivation of pre-existing faults/fractures in the reservoir due to the deep injection is a key concern in designing and running geothermal and water/CO2 injection projects. Therefore, we investigate potential methods to manage injection-induced seismicity. Recent laboratory and field studies recommend that changes in injection pattern (e.g., cyclic injection) might trigger less seismicity than monotonic injection. This study presents results from uniaxial compressive laboratory experiments performed on high porosity Red Felser sandstone that provide new information about the effect of loading pattern and rate on injection‐driven seismicity. Red Felser sandstone samples with identical porosity and dimensions were subjected to three different loading patterns, including cyclic recursive (CR), cyclic progressive (CP), and monotonic loading. Besides, three different loading rates (displacement control) were applied for each loading pattern: low, medium, and high rates that are 10-4 mm/s, 5×10-4 mm/s, and 5×10-3 mm/s, respectively. Microseismicity analysis shows that (i) the maximum magnitude of seismic events and seismic radiated energy at failure decrease for lower loading rates and during the cyclic loading scenario, (ii) the b-value (magnitude-frequency distribution of events) increases on average 40% for a low-rate cyclic recursive loading in comparison with high-rate cyclic recursive and monotonic loading at different rates. The largest b-value resulted from a low-rate cyclic recursive (LCR) loading pattern. The b-value was estimated and compared using different methods, including a least-square regression on either an incremental frequency distribution or a cumulative frequency distribution, and with the maximum likelihood method (MLM) to provide a reliable b-value estimation. The analyses indicate that by considering the accurate magnitude of completeness, MLM, and, with a least-square regression, the incremental frequency distribution, both result in a reliable b-value. From a mechanical perspective, a low loading rate reduces the sample's final strength by 19%. Moreover, samples subjected to cyclic loading display more complex fracture patterns and more disintegration. In our laboratory study, a combination of low-rate loading and a recursive cyclic loading pattern resulted in reduced seismicity through decreasing the maximum seismicity magnitude and increasing the b-value.

The exploitation of subsurface hydrocarbon reservoirs is achieved through the control of production and injection wells (i.e., by prescribing time-varying pressures and flow rates) to create conditions that make the hydrocarbons trapped in the pores of the rock formation flow to the surface. The design of production strategies to exploit these reservoirs in the most efficient way requires an optimization framework that reflects the nature of the operational decisions and geological uncertainties involved. This paper introduces a new approach for production optimization in the context of closed-loop reservoir management (CLRM) by considering the impact of future measurements within the optimization framework. CLRM enables instrumented oil fields to be operated more efficiently through the systematic use of life-cycle production optimization and computer-assisted history matching. Recently, we have proposed a methodology to assess the value of information (VOI) of measurements in such a CLRM approach a-priori, i.e. during the field development planning phase, to improve the planned history matching component of CLRM. The reasoning behind the a-priori VOI analysis unveils an opportunity to also improve our approach to the production optimization problem by anticipating the fact that additional information (e.g., production measurements) will become available in the future. Here, we show how the more conventional optimization approach can be combined with VOI considerations to come up with a novel workflow, which we refer to as informed production optimization. We illustrate the concept with a simple water flooding problem in a two-dimensional five-spot reservoir and the results obtained confirm that this new approach can lead to significantly better decisions in some cases.

We consider fluid-induced seismicity and present closed-form expressions for the elastic displacements, strains, and stresses resulting from injection into or production from a reservoir with displaced faults. We apply classic inclusion theory to two-dimensional finite-width and infinite-width reservoir models. First, we simplify the fault model to the bare minimum while still maintaining its essential features: a vertical fault in a homogeneous reservoir of infinite width in an infinite domain. We confirm and sharpen findings from earlier numerical studies and furthermore conclude that the development of infinitely large elastic shear stresses in a displaced fault, at the internal and external reservoir/fault corners, implies that even small amounts of injection or production will result in some amount of slip or other nonelastic deformation. Another finding is that there is a marked difference between the shear stress patterns resulting from injection and production in a reservoir with a displaced fault. In both situations two slip patches emerge but at the start of injection some amount of slip occurs immediately in the overburden and underburden, whereas during production the slip may remain inside the reservoir region. Next we derive similar but more complicated expressions for displaced inclined (normal or reverse) faults and conclude that our findings for vertical faults also apply to inclined faults. Plain Language Summary Injection of waste water or CO2 in the deep subsurface, or production of natural gas from subsurface reservoirs, may produce earthquakes. Earlier studies have shown that these are especially likely to occur when the reservoir contains faults that have undergone earlier movements (displaced faults). We derive mathematical expressions that allow for an improved understanding of the stresses in these faults compared to earlier computer studies. We conclude, among other findings, that there is an essential difference between injection and production: for injection the fault movement is much more likely to propagate outside the reservoir than for production. Our theoretical insights do not have direct quantitative predictive value but are relevant for the interpretation of experimental and computer studies. They may also help to drastically speed up computer studies, for example, for hazard and risk assessments of injection into or production from deep reservoirs.

We present an efficient workflow that combines multiscale (MS) forward simulation and stochastic gradient computation - MS-StoSAG - for the optimization of well controls applied to waterflooding under geological uncertainty. A two-stage iterative Multiscale Finite Volume (i-MSFV), a mass conservative reservoir simulation strategy, is employed as the forward simulation strategy. MS methods provide the ability to accurately capture fine scale heterogeneities, and thus the fine-scale physics of the problem, while solving for the primary variables in a more computationally efficient coarse-scale simulation grid. In the workflow, the construction of the basis fuctions is performed at an offline stage and they are not reconstructed/updated throughout the optimization process. Instead, inaccuracies due to outdated basis functions are addressed by the i-MSFV smoothing stage. The Stochastic Simplex Approximate Gradient (StoSAG) method, a stochastic gradient technique is employed to compute the gradient of the objective function using forward simulation responses. Our experiments illustrate that i-MSFV simulations provide accurate forward simulation responses for the gradient computation, with the advantage of speeding up the workflow due to faster simulations. Speed-ups up to a factor of five on the forward simulation, the most computationally expensive step of the optimization workflow, were achieved for the examples considered in the paper. Additionally, we investigate the impact of MS parameters such as coarsening ratio and heterogeneity contrast on the optimization process. The combination of speed and accuracy of MS forward simulation with the flexibility of the StoSAG technique allows for a flexible and efficient optimization workflow suitable for large-scale problems.