Poro-Mechanical Coupling for Flow Diagnostics

Accounting for poro-mechanical effects in full-field reservoir simulation studies and uncertainty quantification workflows using complex reservoir models is challenging, mainly because of the high computational cost. We hence introduce an alternative approach that couples hydrodynamics through existing flow diagnostics simulations with poro-mechanics to screen the impact of coupled poro-mechanical processes on reservoir performance without significantly increasing computational overheads. In flow diagnostics, time-of-flight distributions and influence regions can be used to characterise the flow field in the reservoir, which depends on the distribution of petrophysical properties that are altered due to production-induced changes in pore pressure and effective stress. These extended flow diagnostics calculations hence enable us to quickly screen how the dynamics in the reservoirs (e.g. reservoir connectivity, displacement efficiency, and well allocation factors) are affected by the complex interactions between poro-mechanics and hydrodynamics. Our poro-mechanically informed flow diagnostics account for steady-state and single-phase flow conditions based on the poro-elastic theory and assume that the reservoir does not contain fractures. Fluid flow and rock deformation calculations are coupled sequentially. The equations are discretised using a finite-volume method with two-point flux-approximation and the virtual element method, respectively. The solution of the coupled system considers stress-dependent permeabilities. Due to the steady-state nature of the calculations and the effective proposed coupling strategy, these calculations remain computationally efficient while providing first-order approximations of the interplay between poro-mechanics and hydrodynamics, as we demonstrate through a series of case studies. The extended flow diagnostic approach hence provides an efficient complement to traditional reservoir simulation and uncertainty quantification workflows and enable us to assess a broader range of reservoir uncertainties.


Introduction
Coupled poro-mechanical studies are important for the sustainable management of geoenergy reservoirs, containing hydrocarbons, groundwater, or geothermal heat, or providing storage volume for CO 2 , natural gas, or hydrogen storage operations. Reliable, i.e. quantitative coupled poro-mechanical studies with meaningful uncertainty bounds, are needed when the reservoir is expected to deform due to production-or injection-induced changes of the pore pressure and present-day stress field, which may induce seismicity, subsidence, uplift, cap rock failure, fault reactivation, or changes in flow paths (Teufel and Rhett 1991;Dusseault et al. 1996;Fredrich et al. 1996Fredrich et al. , 2000Teatini et al. 2011;Rutqvist 2012;White et al. 2014;Göbel 2015;Feng et al. 2016;Gaitea et al. 2016;Goebel et al. 2017;Kim et al. 2018;Feng et al. 2019;Yuan et al. 2019). Fundamentally, the alteration in flow paths arises when production-induced pressure changes modify the pore structure of the pore space within the reservoir rock, and consequently alter porosity, permeability, capillary pressure, and relative permeability (Rossen and Kumar 1994;Haghi et al. 2018Haghi et al. , 2019Haghi et al. , 2021. Despite the development of different coupling strategies to link poro-mechanical simulations with reservoir simulations (e.g. Park 1983;Armero and Simo 1992;Lewis and Sukirman 1993;Settari and Mourits 1998;Settari and Walters 2001;Minkoff et al. 2003;Kim et al. 2011aKim et al. , 2011bKim et al. , 2011cDoster and Nordbotten 2015;Feng et al. 2016), improvements in computing technologies and numerical algorithms (e.g. Jeannin, et al. 2005Jeannin, et al. , 2007Pettersen and Kristiansen 2009;Pettersen 2012;Mustapha et al. 2016), the significant computing times needed to simulate coupled poro-mechanical processes in full-field reservoir models remains a challenge. The high computational cost of coupling of poro-mechanical simulations with full-field reservoir simulations (e.g. Black-Oil, compositional, or thermal models) renders their application in modern uncertainty quatification workflows, which require hundreds to thousands of simulations, time-consuming and often impractical. As a result, the analysis of how stress-dependent petrophysical parameters can impact reservoir performance remains limited and the uncertainty when forecasting fluid flow in stress-sensitive reservoirs can be high.
For the acceleration of reservoir management workflows, Rasmussen and Lie (2014), Lie et al. (2015) and Møyner et al. (2015) have proposed the use of grid-based flow diagnostics as a computationally efficient complement to traditional reservoir simulations for quantifying the impact of reservoir heterogeneities on fluid flow. Flow diagnostics approximate the reservoir hydrodynamics by computing the time-of-flight distribution in the reservoir and steadystate distributions of conservative tracers that are released at the injectors and, by inverting the flow field, producers. These quantities enable the identification of regions with fast and slow fluid flow, which allow us to estimate influence regions for injectors and producers, drained and swept reservoir pore volumes, connected pore volumes, breakthrough times at individual wells, and well allocation factors. Results from flow diagnostics can therefore be used to select a small subset of reservoir models based on their dynamic performance from a large model ensemble without compromising on the ability to forecast future reservoir performance and estimate the impact of uncertainties (e.g. Caers and Scheidt 2011;Scheidt and Caers 2009;Park et al. 2013;Watson et al. 2021). Our proposed poro-mechanical integration with flow diagnostics frameworks hence provides us with a computationally efficient estimation of how poro-mechanics might impact reservoir performance, which allows us to compare, contrast, and rank different stress-sensitive reservoir models and select individual scenarios for further detailed studies that use computationally demanding full-physics models that couple geomechanical and hydrodynamical simulations. It is important to note, however, that flow diagnostics do not replace full-physics simulations; flow diagnostics simplify the reservoir physics and are based on steady-state solutions, as will be discussed below. Hence flow diagnostics cannot capture how complex coupled processes impact the transient evolution of a reservoir system and subsequent flow processes but only provide information about the endmember states of the reservoir behaviour.
This paper builds upon our recent work (Gutierrez-Sosa et al. 2020) in which flow diagnostics methods (Rasmussen and Lie 2014;Lie et al. 2015;Møyner et al. 2015;Krogstad et al. 2016) are linked with poro-mechanical simulations for single-porosity reservoirs. In the study discussed here, we have extended existing flow diagnostics by first solving a coupled poro-mechanical and hydrodynamical problem before executing the flow diagnostics calculations, and the impact of poro-mechanically altered petrophysical properties on reservoir dynamics is investigated through the computation of the time-of-flight and influence regions. We further have extended and demonstrated the benefits of the analysis of poro-mechanical effects on fluid flow through the extended flow diagnostics simulations. In our proposed framework, the poro-mechanics are derived from the poro-elastic theory and considers stress-dependent-permeabilities. In Gutierrez-Sosa et al. (2022) we have extended the work presented here to fractured reservoirs that can be approximated using dual-continua formulations. The entire framework has been implemented in the opensource MATLAB Reservoir Simulation Toolbox MRST (Lie et al. 2012;Krogstad et al. 2015;Lie 2019). The structure of this paper is as follows. In the first part we introduce the governing equations, coupling strategy, and underlying assumptions. In the second part we present illustrative examples to demonstrate the benefits of including poro-mechanics in the flow diagnostics framework.

Model Formulation
We use the linear poroelastic theory (Biot 1941(Biot , 1955(Biot , 1957 within the macroscopic framework of Coussy (1995Coussy ( , 2004. We limit our work to steady-state single-phase fluid flow and hence the coupling between poro-mechanics and fluid flow manifests itself through the concept of effective stress and a strain dependent porosity and permeability. We briefly introduce the fundamental equation and refer to Gutierrez-Sosa et al. (2020) for further details.

Linear Poro-Elastic Theory
Using linear poro-elastic theory, we assume an isothermal, isotropic, and perfectly elastic material, i.e. we assume a linear and reversible mechanical behaviour with small and quasi-static deformation. The resulting linear interaction between stress and strain and the representation of the linearised strain tensor as a symmetric part of the displacement gradient leads to the governing equation of solid deformation expressed in terms of the displacement field which is given by where G = E(2(1 + )) −1 is the shear modulus, = E (1 + ) −1 (1 − 2 ) −1 is the Lamé constant with Young's Modulus E and the Poisson's ratio , is the displacement field, is the Biot coefficient of the material, p is the pore pressure, is the density of the solid material and b are the body forces.

Fluid-Flow Theory and Poro-Mechanics
Under the assumption of steady-state conditions, the equations governing isothermal single-phase fluid flow in a deformable porous medium (i.e. the solid velocity is nonzero, s ≠ 0 ) are given by Jaeger et al. (2006) as Darcy's law describes fluid flow relative to the solid phase flow as where is the density, is Darcy's velocity, s is the solid velocity, is the porosity, q is the source/sink term, is the intrinsic permeability tensor of the rock, is the dynamic fluid viscosity, g is the gravity vector, z is the vertical coordinate and the subscripts f and s refer to the fluid and solid, respectively.
Under the assumption of steady-state conditions and incompressible fluid flow by combining Eqs. 2 and 4 and substituting ∇ ⋅ s = 0 from Eq. 3 yields the traditional definition of the incompressible fluid pressure equation Equations 1 and 5 represent the constitutive model for steady-state fluid flow and solid deformation.

Stress-Dependent Rock Properties
We consider a stress-dependent relation between porosity and permeability. We represent the variation in porosity as a function of a normalised change of the volumetric strain (Tortike and Farouq, 1993;Li et al., 2006) and model the stress-dependent permeability using the modified Kozeny-Carman equation (Kozeny 1927;Carman 1956;MacMinn et al. 2016). Assuming compaction to be positive, the relationships between volumetric strain, porosity and permeability are given by where Δ v is the change in volumetric strain (calculated from the initial stress and pore pressure conditions of the simulation to the current deformation condition), ini is the initial porosity, ini represents the initial permeability tensor, and ( ) and ( ) are the ini updated stress-dependent porosity and permeability, respectively. For simplicity, we neglect any induced anisotropy in permeability. Equation 7 is the nonlinear coupling term that links rock deformation (Eq. 1) to changes in fluid flow (Eq. 5). Note that Eqs. 6 and 7 can readily be replaced by other constitutive relationships that evaluate how stress changes influence porosity and permeability, including correlations obtained from laboratory experiments. The numerical efficiency of the flow diagnostics calculations therefore allows us to not only evaluate how geological uncertainties impact reservoir flow but also how uncertainties in the poro-mechanical data impact subsequent reservoir dynamics.

Grid-Based Flow Diagnostics
Grid-based flow diagnostics (Rasmussen and Lie 2014;Lie et al. 2015;Møyner et al. 2015) approximates dynamic behaviour of the reservoir model in a fraction of the time needed by full-physics simulations by using efficient finite volume discretisation (Natvig et al. 2006(Natvig et al. , 2007Natvig and Lie 2008;Aarnes et al. 2007;Eikemo et al. 2006Eikemo et al. , 2009. Flow diagnostics can compute the reservoir partition, well-allocations, and swept and drained reservoir volumes based on the distributions of a conservative tracer concentration c and the time-offlight (or TOF). The time-of-flight corresponds to the time a non-reactive particle takes to travel from an injector to a certain point in the reservoir or from this point to a producer at a given velocity field in the reservoir. In flow diagnostics simulations the following equations are solved The solution of the concentration field (Eq. 8) is computed for each injector-producer pair (note that the flow field is reversed when evaluating the distribution of c for the producers). Equation 9 computes the distribution of in the reservoir and the volume-averaged values describe the flow within each well pair region. However, a drawback is that becomes inaccurate in regions where there are areas of fast and slow flow (i.e. small and large values of ) or where flow associated with different well pairs merges into a single grid cell. Hence Lie (2019) has proposed Eq. 10 as a slight modification of Eq. 9 to improve the computation of between each injector-producer pair, resulting in a more accurate computation of the field and breakthrough time for each influence region. , therefore, shows the impact of geological heterogeneity upon flow behaviour, enabling us, for example, to quickly identify reservoir regions with poor sweep efficiency or at risk of early breakthrough (Thiele et al. 2007(Thiele et al. , 2010. The threshold of the individual concentration fields of the volumetric regions can be approximate through the identification of the grid cells that have not met an established conditional of maximum value max (Eq. 11) in the reservoir model. This relationship is given by where PV is the pore volume, and f is the finite-volume discretised flux, i.e. the volumetric flow rate in a grid cell as defined by the permeability field of the reservoir and the given inflow boundaries, source terms, wells, or combinations of these. The subindex i refers to the i th grid cell. The advancement of the injected fluid (fluid displacement front) across the reservoir can therefore be tracked trough max . Assuming that the present flow field remains constant (steady-state condition), the evolution of the fluid displacement front can be estimated as a function of an equivalent actual time t e through the relationship of the total injected pore volume PV inj constrained by max with respect to the total injected flow rate ∑ n j q inj j of all j th injection wells in the model. This time t e is given by The dynamic Lorenz coefficient L c and F − Φ curves, as defined by Shook & Mitchell (2009) in the context of streamlines and Shahvali et al. (2012) in the context of finite volumes, are the common metrics that flow diagnostics use to characterise the dynamic heterogeneity of the flow field and displacement process. F , Φ , and L c depend on and the pore volume influenced by the displacement process. In this paper we refer to F , Φ and L c based on their finite-volume discretisation as where subscripts j and i refer to the cell index sorted according to ascending total time-offlight (fast to slow) and N is the total number of grid cells. F and Φ values are normalized so that both are within the unit interval. Note that L c = 0 indicates ideal piston-like displacement and L c = 1 corresponds to infinitely heterogeneous displacement.
Since the poro-mechanical change in petrophysical properties affect the velocity field (Eq. 1, 5, and 7) in the reservoir, the spatial distributions (x) and c(x) (Eqs. 8 and 9) will also be influenced and hence all subsequent flow diagnostics estimates will vary in response to any poro-mechanical changes.

Numerical Formulation
Our steady-state poro-mechanical coupling module has been implemented as part of the open-source MATLAB Reservoir Simulation Toolbox MRST (Lie et al. 2012;Krogstad et al. 2015;Lie 2019). The poro-mechanical problem described in Eqs. 1 and 5 is discretised using the virtual element method (VEM) and the two-point flux-approximation (TPFA) scheme, respectively, in MRST. Our implementation leverages the incompressible flow and mechanics solver modules that come as standard within MRST. As presented by Gutierrez-Sosa et al. (2020), the discretised coupled system of equations is given by where c and c are control point values of the linear displacement and pressure field variable vectors, is the transmissibility matrix, p is the vector of fluid body forces and fluxes, is the stiffness tensor, is the interaction matrix, u is the vector of body forces and surface traction. The coefficient matrices and contain mechanical parameters of the rock, the transmissibility matrix contains the stress-dependent permeability (Eq. 7). The coupled system of equations becomes a nonlinear system because the coefficient matrix is dependent on the permeability unknowns ( ) . The solution of the nonlinear system of equations (Eq. 16) is achieved by iterating sequentially between each MRST module through a new interface module which includes stress-dependent rock properties correlations that act as a coupling term. Numerically, this approach is represented by.
where superscript r refers the iteration level. We use a fixed-point iteration which starts with an initialisation stage where the flow and mechanics models are set to their initial conditions. The coupling loop starts with the computation of the fluid pressure, assuming that the initial conditions of the displacement fields c,r are not changing. With the solution of the pressure c,r+1 * , the mechanics solver is updated to account for the effect of the fluid pressure on the deformation of the rock. The solution for c,r+1 * is used to obtained v c,r+1 * , and ultimately the stress-dependent permeability c,r+1 * is calculated to recompute the transmissibility matrix ( c,r ) for the next calculation of the fluid pressure c,r+1 . This procedure is repeated until ‖ c,r+1 * − c,r+1 ‖ ≤ , i.e. convergence has been reached. Once the solution has converged, the final updated stress-dependent permeability and its corresponding pressure field is used to compute the velocity fields from which flow diagnostics can be calculated. We note that we chose pressure rather than displacement as the convergence criterion because we need to obtain an accurate solution of the pressure field for the subsequent flow diagnostics calculations, which depend on the accurate estimate for the fluid pressure (Eq. 5). It therefore could be possible that the displacement field still contains a significant residual error if the tolerance is too large. For further details about the demonstration and validation of our proposed poro-mechanics scheme we refer interested readers to the work of Gutierrez-Sosa et al. (2020).

Simple Case -Box Model
The integration of poro-mechanics with flow diagnostics is first illustrated by studying an idealised 5-spot box model with dimension of 300 m × 300 m × 50 m, consisting of 21 × 21 × 20 grid cells. All wells are bottom hole pressure (BHP) constrained. The model has homogeneous petrophysical properties with poro-elastic behaviour, and the system is confined and closed. For simplicity and without loss of generality, there are no acting forces or loads. The effective stress is only affected by the production-injection process. Note that the proposed framework is not limited to these boundary conditions and more complex boundary conditions with diverse displacements and applied loads are possible. The reservoir domain is divided into three regions (Fig. 1). The central region of the model (diagonal) represents a mechanically soft material, and the off-diagonal regions are mechanically stiffer. The properties of the model, boundary conditions, and input parameters are given in Fig. 1 and Table 1. This simple example does not represent a geologically realistic scenario and only intends to illustrate and compare the effect of productioninduced disturbance of stress-dependent petrophysical properties on reservoir dynamics when wells are allocated in regions of different mechanical stiffness, i.e. two injectors are placed in the stiff regions and the remaining injectors and producer in the soft region.
Before starting production-injection operations, the permeability of the model is initialised to the stress state obtained from the established mechanical boundary conditions, initial reservoir pressure and input properties and parameters. Then, this "mechanically initialised permeability" is used for the hydrodynamical simulation when poro-mechanics are neglected, i.e. no computation of coupled poro-mechanical solution and no production-induced permeability alteration (referred to as w/o poro-mechanics). For the simulation when poro-mechanics are considered, i.e. the computation of coupled poro-mechanical and hydrodynamical solutions with resulting permeability alteration (referred to as w/ poro-mechanics), the mechanically initialised permeability represents the start of the simulation before the stress state is altered by the production-injection operation. In this manner, having defined a common permeability is intended to fairly compare and quantify the resulting production-induced changes when accounting for poro-mechanics against the simulation case that neglects poro-mechanics.
When considering poro-mechanics, permeability decreases by up to 25% in the softer region of the reservoir, specifically in the vicinity of producer where the pressure drop is the largest (Fig. 2). Consequently, the well influx and are affected by the permeability change. The overall flux in the model is reduced, with the most significant reduction occurring between the injectors and the producer located in the stiff zone, i.e. wells INJ1, INJ4, and PROD1. Furthermore, the change in flow field impacts and the stationary producer-and injector concentration which alters the reservoir volumetric partitioning and well allocation factors. The result of these changes is the modification of Fig. 1 Setup of the 3D poro-elastic problem. Model with homogenous and isotropic petrophysical properties, representing a system with lateral confinement and vertical constrained movement at the bottom and top (a) with two Young's modulus regions that represent the materials stiffness and divide the system into three mechanical regions, two stiff regions and one soft region (b) reservoir connectivity and displacement efficiency of each injector (Fig. 4). We refer to produced concentration c p as the ratio between the individual concentration flow rate q p contribution of each injector i to the total produced concentration flow rate ∑ n i q p i of a given producer p which is given by  The contribution of each injector to the produced concentration primarily depends on the reservoir connectivity between well pairs. The produced concentration, therefore, reveals the evolution of pore volume connectivity over time, i.e. influence regions between well pairs as a function of time max (Eqs. 10 and 11), assuming that the present flow field remains constant. Note that when considering poro-mechanics, the produced concentration of each well pair (i.e. INJ2-PROD1 and INJ3-PROD1 versus INJ1-PROD1 and INJ4-PROD1) is different compared to the case where poro-mechanics are neglected.
In addition, we analyse the impact of including poro-mechanics at the well level in three different ways (Figs. 3 and 4) by measuring (i) the cumulative flux from the bottom to the top of the perforated layers, (ii) the individual flux per perforated layer, and (iii) by comparing production and injection profiles. Figure 3 shows that the well flow rates are  Fig. 4f are caused by the strong anisotropy in the flow field, the relatively coarse grid that is not aligned with the main flow direction, and the two-point flux approximation, all of which influence the accuracy at which Eqs. 8 and 9 are solved and the fluxes at different max values are evaluated reduced in terms of inflow per reservoir layer and cumulative flow rates (up to 12%). Note that the reduction of flow rate per layer is uniform due to the absence of gravitational and capillary forces and because the system experiences the same deformation in all layers, which implies that the permeability changes in each layer occur in the same proportion. The reduction in inflow rates (Fig. 4) causes an increase in the time it takes to produce the reservoir (Fig. 5). Note that same values of max in Figs. 4 and 5 refer to two different time scales (Eq. 12) and injected pore volume because changes due to the alteration in permeability and flow rates caused by poro-mechanics. The equivalent actual time t e and the injected pore volume are hence depicted as secondary axes on Figs. 4 and 5.

Effect of Different BHP Constraints
To further illustrate poro-mechanical effects under different flow conditions, we run additional simulations using the same box model described in Fig. 1 and Table 1 but consider three cases with different BHP constraints defined by values of 1, 10 and 100 bar above and below the initial average pressure p i for the injectors and the producer, respectively. This example aims to illustrate under which production conditions poro-mechanical effects have a noticeable impact on reservoir performance.
Under the prescribed operational conditions, only case 3 (BHP constraints are p i ± 100 bar, respectively) shows a considerable impact of poro-mechanical effects on the well allocation factors (Fig. 6). We observe significant changes which include different arrival times of each injector concentration at the producer and a substantial variation in the flux allocation of the different well pairs. This variation in the production profile of the allocation fluxes indicates a significant alteration in the preferential flow paths and reservoir connectivity. In contrast, poro-mechanical effects on reservoir performance appear to be negligible in the other two cases. This example therefore illustrates how the extended flow diagnostics framework could be used to screen the impact of poro-mechanical effects on reservoir dynamics and select individual scenarios, i.e. case 3, for further detailed fullphysics simulations.
We observe that the integration of flow diagnostics with poro-mechanical simulations enables us to compute the essential reservoir dynamics much faster than the full-physics  Fig. 5a are numerical simulations carried out with a commercial simulator. We compared our simulations with a fully coupled commercial simulator and observed a difference of approximately 180 (approximately 1 min and 3 h for our proposed implementation and the commercial simulator, respectively) to reach the steady-state solution but note that we did not aim for a rigorous comparison of CPU times. Hence, our approach is ideal for a preliminary screening of possible reservoir behaviours before commencing more detailed and CPU intensive fullphysics poro-mechanical simulations that better capture the transient evolution of the reservoir system due to the complex, nonlinear, and time-dependent interactions of the coupled hydrodynamical and poro-mechanical processes.

General Overview of the SPE10 Models
The 10 th SPE Comparative Solution Project, also known as SPE 10 model (Christie and Blunt 2001), considers a complex and heterogeneous reservoir geology represented by a regular Cartesian grid with box geometry (Fig. 7). The model dimensions are 366 m × 671 m × 52 m discretised into 60 × 220 × 45 grid cells, it represents a part of the Brent sequence, consisting of the Tarbert formation (first 20 layers) and the Upper Ness formation (remaining 25 layers). The Tarbert formation model represents a prograding near-shore environment, consisting of stacking patterns of sand bodies of different characteristics (Fig. 7) and the Upper Ness formation model represents fluvial channels bounded by shales (Fig. 7). As in Gutierrez-Sosa, et al. (2020), we rescale the original petrophysical properties with the purpose of inducing more significant changes in the poro-mechanical simulations, as suggested by Rutqvist and Stephansson (2003). We preserve the original porosity-permeability relationships for the four facies (fine sand, sand, coarse sand, and shale). Mechanical properties such as Young's modulus, Biot's coefficient and rock density are assigned to each facies in the model using the data from Graham (1997) for consolidated sandstones and Molina et al. (2017) for shales. Thus, the model is mechanically heterogeneous. The porosity and permeability ranges and the Young's modulus distribution across the facies for the Tarbert and Upper Ness formations are depicted in Fig. 7 and Table 2.  6 Production profile for: case 1 with BHP constraints of p i ± 1 bar (a), case 2 with BHP constraints of p i ± 10 bar (b), and case 3 with BHP constraints of p i ± 100 bar (c) for the injector-producer pairs when neglecting and accounting for poro-mechanics. The secondary axes in the upper part of each plot represent the equivalent actual time t e and the cumulative injected volume, assuming the present flow field in the reservoir remains constant. Note that fluid flow is symmetric when neglecting poro-mechanics and hence the corresponding curves overlap We study Tarbert and Upper Ness formations separately due to their different characteristics. However, both formations are examined under the same operational consideration, assuming a 5-spot injection pattern where all wells are BHP constrained. Each model is subject to production-induced poro-mechanical effects. There is no lateral displacement on any of the vertical sides, the bottom of the model cannot move vertically, and the top is free to displace in all directions. In a similar way to the box model case, prior to commencing production-injection operations the permeability of Tarbert and Ness formations is subjected to the initial stress state according to the defined mechanical boundary conditions, initial reservoir pressure, input properties, and parameters listed in Tables 2 and 3. The hydrodynamical simulation when poro-mechanics are neglected utilises the mechanically initialised permeability, which remains unaltered during the whole simulation. Thus, the permeability before stating production operations for the simulation case that consider poro-mechanics is identical to the simulation case that neglects poro-mechanics.

Case Study 1: Tarbert Formation
A considerable reduction of the reservoir permeability, as well as a noticeable change in permeability distribution, can be observed when considering poro-mechanics in the simulations (Fig. 8). The impact of the rock stiffness on the permeability change is most   Fig. 7), leading to a permeability reduction of up to 300 mD. Furthermore, changes in flow paths between well pairs (see well pair INJ3-PROD1 in Figs. 9 and 10) and well productivity can be observed. Well inflow rates in each layer are significantly altered, leading to a 50% reduction in overall inflow at the producer (Fig. 11). The reduction of reservoir productivity and injectivity results in a delay in the breakthrough time ( Fig. 12) and subsequently slower recovery, indicating that the complex interplay of poro-mechanics and reservoir dynamics warrants further analysis using full-physics simulations.

Case Study 2: Upper Ness Formation
In the Upper Ness formation, poro-mechanical effects only cause slight changes resulting in a subtle overall reduction of the reservoir permeability (Fig. 13). However, this small reduction influences the change in reservoir connectivity (Figs. 14 and 15), causing a 32% reduction in inflow rate and significant modification of the inflow profile (Fig. 16). The reduction of reservoir productivity and injectivity results in longer oil recovery times and a slight delay in breakthrough times (Fig. 17). The substantial change in reservoir productivity suggests that even if the heterogeneity of mechanical properties cause only minor changes in permeability, the impact on reservoir performance can be significant and, as for the Tarbert formation, this behaviour should be studied in further detail using the appropriate full-physics simulations.

Discussion of Case Studies
We compare the impact of poro-mechanics on fluid flow for Tarbert and Upper Ness formations using the F − Φ curves, dynamic Lorenz coefficient L c , and the permeability quotient to define the change in the spatial distribution of permeability using the grid cell-based quotient operator (Henrici 1964) in its percentage form Q x = x(u)∕x ini − 1 × 100 , which quantifies the percentage of change between the poro-mechanically-altered quantity x(u) to the original quantity x ini . The poro-mechanically induced changes in the Upper Ness and Tarbert formations models cause L c to increase, which indicates an increase in the dynamic heterogeneity and decrease in sweep efficiency (Fig. 18). As discussed above, the increased dynamic heterogeneity, which changes flow paths and connected pore volumes, also alters the production profiles (Figs. 12 and 17). The change in flow paths (Figs. 9 and 14) are a direct consequence of the spatial changes in the permeability, which lead to permeability quotients ±  Fig. 17a are numerical 8% for the Upper Ness formation and − 100% and + 30% for the Tarbert formation (Fig. 18). Although the permeability quotient for the Upper Ness formation is small compared to that of the Tarbert, the Lorenz coefficient increases by 23% for the Upper Ness but only 6% for the Tarbert formation, this counter-intuitive behaviour of the permeability reduction and changes in L c indicates that the impact of poro-mechanics on the reservoir dynamics is more severe for the Upper Ness formation where the heterogeneity in permeability is more randomly distributed in comparison to the Tarbert formation with its extensive and continuous sand bodies, which results in a more uniformly distributed permeability. Hence, the relatively small permeability change in the Upper Ness formation causes more pronounced alterations in flow paths and displacement efficiency. These results clearly indicate how important the inclusion of poro-mechanics can be for reliable reservoir performance forecasts because spatial changes in permeability can affect flow paths, connected pore volumes, and ultimately the recovery efficiency in unexpected and counter-intuitive ways.
Although the obtained results can be simulated using traditional coupled reservoir simulators, our proposed methodology provides a computationally efficient first screening to assess the likely importance of poro-mechanics at a significantly reduced computing time prior to more detailed reservoir simulation studies. For this particular example, the entire computations were carried out on a standard desktop PC for the Tarbert and Ness formations; it took only 16 and 23 min to simulate the whole workflow, respectively. Such fast computations allow us to analyse a much wider range of geological scenarios, parameter combinations, and well patterns, and hence enable us to screen and explore a broader range of uncertainties prior to selecting models and scenarios for more detailed full-physics simulations using appropriate clustering and ranking techniques (Caers and Scheidt 2011;Scheidt and Caers 2009;Park et al. 2013;Watson et al. 2021). Performing similar screening simulations would likely take days using coupled full-physics simulations. Hence the proposed extended flow-diagnostics framework is a natural pre-processing that helps to accelerate, but does not replace, coupled process modelling and modern uncertainty quantification workflows for poro-mechanical studies. Fig. 18 Comparison of F − Φ curves for the Tarbert formation and Upper Ness formation (a), and the permeability quotient for the Tarbert formation (b) and Upper Ness formation (c) to compare the changes in permeability due to production-induced poro-mechanical changes

Conclusions
We demonstrate an extended flow-diagnostics framework that integrates poro-mechanics with traditional flow diagnostic calculations that approximate the reservoir dynamics in unfractured reservoirs. The proposed poro-mechanically informed flow diagnostics framework establishes a sequential poro-elastic coupling (Eqs. 17 and 18) between the steady-state fluid flow and the rock deformation problem (Eqs. 1 and 5) considering stressdependent permeabilities (Eq. 7) that act as a coupling term. By solving the poro-mechanical problem first and then executing the computationally efficient flow diagnostics simulations, the impact of the change in the reservoir flow field due to the poro-mechanically altered petrophysical properties is investigated. The spatial distribution of the time-of-flight and the stationary producer-or injector concentrations (Eqs. 8 and 9) is utilised to characterise reservoir dynamics, sweep efficiency, and well inflow rates. This framework was implemented in the open-source MATLAB Reservoir Simulation Toolbox MRST.
Using two case studies, we demonstrated the capability of our extended flow diagnostics framework to quickly screen how the complex interactions between geomechanical and hydrodynamical processes could alter petrophysical properties and hence affect subsequent predictions of the reservoir dynamics. The extended flow-diagnostics framework can therefore be used to complement and accelerate modern coupled process simulation studies that aim to assess and quantify the impact of a broader range of geomechanical and hydrodynamical uncertainties as well as engineering parameters (e.g. well placements) on the reservoir performance. By using the extended flow diagnostics framework, it is now possible to quickly screen, compare and contrast, and rank different reservoir models so as to identify a smaller number of representative reservoir models that need to be taken forward for more detailed, but also more time-consuming full-physics simulation studies that properly represent the transient coupled processes, while still honouring the full range of uncertainties inherent to the reservoir.