Controls on the spatio-temporal patterns of induced seismicity in Groningen constrained by physics-based modelling with Ensemble-Smoother data assimilation

The induced seismicity in the Groningen gas ﬁeld, The Netherlands, presents contrasted spatio-temporal patterns between the central area and the south west area. Understanding the origin of this contrast requires a thorough assessment of two factors: (1) the stress development on the Groningen faults and (2) the frictional response of the faults to induced stresses. Both factors have large uncertainties that must be honoured and then reduced with the observational constraints. Ensembles of induced stress realizations are built by varying the Poisson’s ratio in a poro-elastic model incorporating the 3-D complexities of the geometries of the Groningen gas reservoir and its faults, and the historical pore pressure distribution. The a priori uncertainties in the frictional response are mapped by varying the parameters of a seismicity model based on rate-and-state friction. The uncertainties of each component of this complex physics-based model are honoured through an efﬁcient data assimilation algorithm. By assimilating the seismicity data with an Ensemble-Smoother, the prior uncertainties of each model parameter are effectively reduced, and the posterior seismicity rate predictions are consistent with the observations. Our integrated workﬂow allows us to disentangle the contributions of the main two factors controlling the induced seismicity at Groningen, induced stress development and fault frictional response. Posterior distributions of the model parameters of each modelling component are contrasted between the central and south west area at Groningen. We ﬁnd that, even after honouring the spatial heterogeneity in stress development across the Groningen gas ﬁeld, the spatial variability of the observed induced seismicity rate still requires spatial heterogeneity in the fault frictional response. This work is enabled by the unprecedented deployment of an Ensemble-Smoother combined with physics-based modelling over a complex case of reservoir induced seismicity.

Induced seismicity in the Groningen gas field 1283 The induced seismicity rate in Groningen is potentially controlled by the combination of two processes: (1) the development of the induced stress changes on faults during gas production and (2) the frictional response of each fault when subjected to these induced stresses. In particular, depending on the fault frictional response, strong and fast stress changes might potentially lead to only low seismicity rates. Recently, Candela et al. (2019) proposed that spatio-temporal variabilities of the observed Groningen seismicity rate could be entirely explained by heterogeneities in the fault frictional response. However, spatial heterogeneities in elastic properties [as a consequence of the variations in the reservoir porosity (NAM 2016) and indirectly revealed by subsidence inversion (Fokker & Van Thienen-Visser 2016;Smith et al. 2019)], which can directly control the development of induced stress changes, were disregarded in their modelling approach. The present contribution deploys a modelling strategy specifically designed to disentangle the relative importance between induced stress development and frictional response in controlling the induced seismicity at Groningen. To do so, a transparent modelling strategy is developed which allows to dissociate the contribution of each physical process.
When the objective is to disentangle the relative contribution of each physical process at work, intrinsically the modelling complexity increases as well as its computational demand and thus its run time on a standard PC. This calls for an efficient data assimilation procedure to reduce the number of models needed to appropriately constrain the posterior distribution of each model parameter for each physical process at work. Indeed, such complex and computationally heavier modelling strategy precludes the use of traditional brute-force or Markov Chain Monte Carlo (Foreman-Mackey et al. 2013) methods to screen the prior space of each model parameter. The present contribution develops an Ensemble-Smoother formulation with a single-step of data assimilation (Emerick & Reynolds 2013a,b) which allows to constrain the posterior distributions of each model parameter with a relatively small ensemble of models. Combining both (i) a transparent physics-based modelling strategy honouring all the available a priori knowledge and (ii) an efficient and robust data assimilation procedure, the causes of the spatiotemporal evolution of the seismicity at Groningen are effectively isolated and identified.

A T R A N S PA R E N T P H Y S I C S -B A S E D M O D E L L I N G S T R AT E G Y
The deployed modelling strategy is tailored to honour both: (i) the physics of the processes at work and (ii) all the pre-existing a priori knowledge. The starting modelling ingredient is the flow simulation of the entire Groningen gas field, computed by the field operator, incorporating all the prior knowledge in terms of geology and subsurface hydrogeology, and history-matched with the subsidence observations (Bierman et al. 2015;Fokker & Van Thienen-Visser 2016;Smith et al. 2019). The 3-D numerical flow model for the full field contains approximately 650 000 cells, with an average cell size of 725 × 725 m and cell thickness of 5-10 m.
From this discretized pore pressure evolution, the second modelling step consists in computing the induced stresses at the Groningen faults. We make use of the same semi-numerical approach as Candela et al. (2019), the so-called MACRIS method (Mechanical Analysis of Complex Reservoirs for Induced Seismicity, see van Wees et al. 2018van Wees et al. , 2019. MACRIS is a mesh-free approach in which there is no need to build a dedicated grid for the geomechanical analysis. The 3-D yearly pressure fields of the Groningen field are thus directly used as input, and each cell is considered as a compacting nucleus of strain (centre of compression ;Mindlin 1936;Geertsma 1973;Okada 1992). One key aspect of MACRIS is the use of the Barnes-Hut algorithm (Barnes & Hut 1986;van Wees et al. 2019) to rediscretize the initial reservoir grid for two purposes: (i) clustering the nuclei of strain close to the faults in order to increase the spatial resolution of stress and (ii) shortening the computation time. The induced displacement vector at a receiver point (hereafter also called observation point) caused by a compacting nucleus of strain (i.e. a point source) is defined as: where u i are displacement vector components (u 1 , u 2 , u 3 ), x 1 and x 2 are the x,y-coordinates of the receiver point relative to the point source located at (0,0,-c). x 3 and -c are z-coordinates of the receiver point and point source, respectively, relative to the surface. The term M 0 2π G represents the strength of the source, with the moment M 0 defined as: where P is the reservoir pressure change (with a negative value denoting depletion) in volume dV, represented by the point source, and the Biot coefficient α Biot is considered equal to unity. The strength of the point source M 0 2π G is thus dependent on the elastic properties (G stands for shear modulus) of the medium, here corresponding to an elastic homogeneous reservoir and burden. The shear modulus is related to the Young's modulus and Poisson's ratio by: where E and ν are the Young's modulus and Poisson's ratio, respectively. The first term on the right-hand side of eq. (1) represents the displacement field due to an inflation force in an infinite medium. For receiver positions far from the surface and close to the source, this term dominates. The second term of eq. (1) originates from a similar image source above the surface; the third and fourth terms also originate from the image position but have a different form. The contribution of each of these nuclei is integrated to compute the induced poro-elastic strain changes along each Groningen fault with a metre-scale spatial resolution. More precisely, the receiver points are chosen along the pillars of the reservoir mesh corresponding to the faults (Fig. 1). Along these fault pillars, the receiver points are regularly spaced every 5 m and extended along the fault dip at a distance of up to a few hundred metres above and below the reservoir. To convert induced poro-elastic strain changes ε into induced stress changes σ at the receiver points, the generalized Hooke's law is applied as: and with In order to account for the direct effect of the pore pressure on the effective normal stress at faults, the pore pressure at the faults is sampled from the fault compartments which have experienced the largest pressure changes. In a nutshell, MACRIS offers the unique opportunity to honour: (1) the heterogeneous spatio-temporal distribution of the pressure fields, (2) the fault geometry including its large-scale roughness and offset juxtaposing non-depleting rocks against the depleting reservoir, (3) the poro-elastic stress changes and (4) the direct pore pressure effect at faults. The restriction still present in MACRIS is that only one-way coupling is considered. We deem this acceptable for gas reservoirs, where the effect of compaction on the gas pressures in the pores is small.
Combining both the poro-elastic and direct pore pressure effects, the modified Coulomb stress function can be calculated as: Downloaded from https://academic.oup.com/gji/article/229/2/1282/6459727 by guest on 06 April 2022 where τ is the shear stress acting along the fault plane, σ n is the effective normal stress, μ the friction coefficient and α is a constitutive parameter (zero in this study). All the model parameters used in MACRIS up to the Coulomb stress calculation are displayed in Table 1. These model parameters are representative of the Groningen field (e.g. van Eijs 2015). The reader is referred to appendix A1 (and results of van Wees et al. 2019) to evaluate the excellent match between the MACRIS stresses and the solution obtained with a commercial finite element code (DIANA FEA, http://dianaf ea.com).
MACRIS gives access to the stress changes along the 3-D fault surfaces of the Groningen field. Fig. 1 displays the Coulomb stress changes along the 800 Groningen faults, corresponding to more than two million observation points, distributed over about 21 000 fault pillars. Because of the lack of constraint on the depth of the observed events, the seismicity catalogue can be seen as a 2-D field evolving over time. Consequently, for model-data comparison, the MACRIS output needs to be represented in 2-D through time as well. There are a number of ways in which the upscaling from 3 spatial dimensions to 2 spatial dimensions can be achieved. The 'strictly correct approach' consists in calculating first the seismicity rate at each receiver point along the length and width of the Groningen faults and then calculating the integral over each pillar to get the seismicity rate per pillar in 2-D. To gain in computational speed, we followed the approach of Candela et al. (2019) to represent each 3-D fault pillar by its observation point recording the maximum Coulomb stressing rate. Appendix A1 of Candela et al. (2019) demonstrates that this approximation captures the temporal evolution of the seismicity rate at each pillar location. This can be explained by the fact that for each pillar the integrated seismicity rate for all the observation points is expected to be dominated by the location where the seismicity rate is the highest (Kroll et al. 2017). The downside of this approach is that stress peaks at the reservoir edges can control the maximum Coulomb stressing rate. However, it is not possible to relax these stress peaks which are intrinsic to our elasticity assumption. Even for the 'strictly correct approach', these stress peaks are part of the integral.
From the 2-D upscaled induced stress changes on the fault, the seismicity rate can be computed incorporating the fault frictional response. Ignoring the fault frictional response, the traditional Coulomb failure model predicts that whenever the Coulomb stress reaches a threshold value, an earthquake is generated. This prediction is not in agreement with the observed seismicity, which generally shows a gradual decay following the onset of Coulomb stress decrease. Instead, the rate-and-state formalism reproduces the fact that the onset of frictional sliding is a non-instantaneous time-dependent process (as opposed to the instantaneity assumption of the Coulomb model), which introduces a time-dependent failure mechanism for the generation of earthquakes. Assuming a population of faults following a rate-and-state frictional behaviour, and where the time-to-failure of the nucleation spots along the faults is uniformly distributed, Dieterich (1994) derived the following seismicity rate model: and where R D is the seismicity rate, θ is a state variable, S is the modified Coulomb stress function defined in eq. (7). The constant r 0 is the steady-state background seismicity rate at the reference stressing rateṠ 0 . A is a dimensionless fault constitutive parameter. Segall & Lu (2015) reformulated this seismicity rate equation to eliminate the state variable θ. They defined a normalized seismicity rate, relative to the background rate, as: The differential equation for R, derived from eqs (8) and (9) (see Appendix A2 for the detailed derivation), is: where the Coulomb stressing rateṠ is defined as: In eq. (10), t a = Aσ n /Ṡ 0 is the characteristic time delay for the earthquake nucleation process. Note here that when solving the ordinary differential equation (eq. 10), and as implemented by Rubin & Ampuero (2007), the changes in effective normal stress during the induced stressing history are honoured. Indeed, the changes in effective normal stress are accounted for at three levels: (1) the Coulomb stress rateṠ, (2) the friction (defined as μ = τ/σ n ) and (3) t a .
In summary, the physics-based modelling strategy involved two key steps: (1) computing the development of induced stress changes along the Groningen faults during reservoir production and (2) computing the seismicity rate taking as input the induced stress changes and honouring the fault frictional response. Both the induced stress development and the fault frictional response can control the spatiotemporal evolution of the seismicity rate at Groningen. In the sequel of the manuscript, the induced stress development refers to the MACRIS calculation, and the fault frictional response is embedded in Dieterich's seismicity rate theory. Candela et al. (2019), instead of modelling the entire Groningen field, focused on two subareas: (i) the central area (C-area) where the seismicity rates are the highest and (ii) the south west area (SW-area) which has a much lower seismic activity (see Fig. 2). They concluded that the difference in seismicity rate between these two subareas can be explained by a difference in the frictional response of the faults. This conclusion was based on the assumption of spatially homogeneous elastic properties when computing the poro-elastic stress changes. However, it is well-known that the porosity of the Groningen reservoir is spatially heterogeneous. More specifically, Fig. 2 displayed the weighted averaged porosity of the Groningen reservoir flow model used as input for MACRIS; at the C-area location the Groningen reservoir is relatively more porous compared to the SW-area. This porosity contrast attests for the reservoir-scale heterogeneities in the uniaxial compaction coefficient C m (NAM 2016), which is linked to the reservoir Young's modulus and Poisson's ratio by: Inversions of subsidence measurements have also highlighted the spatial heterogeneities in the compaction coefficient of the Groningen field (Bierman et al. 2015;Fokker & Van Thienen-Visser 2016;Smith et al. 2019). In the present contribution, we test the hypothesis that the spatial heterogeneities in the Poisson's ratio can explain the difference in the seismicity history between the C-area and SWarea. Candela et al. (2019) assumed an identical Poisson's ratio of 0.2 for both areas of interest. However, using a lower Poisson's ratio for the C-area should lead to a much steeper stress path (Segall & Fitzgerald 1998;Hettema et al. 2000;Buijze 2020), and thus to relatively higher Coulomb stress rates. For uniaxial strain conditions (that is considering a laterally extensive depleting reservoir undergoing uniaxial compaction), the stress path is defined as: where σ h is the change in horizontal stress and α Biot = 1 as stated previously for the MACRIS calculation. The increase in horizontal stress and thus the increase in Coulomb stress scale linearly with the increase in the stress path caused by a lowering of the Poisson's ratio.
Assuming now that the frictional response of the Groningen faults is spatially uniform (that is, identical for both the C-area and the SWarea), these higher Coulomb stress rates might explain the earlier kick-off of the induced seismicity in the C-area relatively to the SWarea (see Fig. 2). Note here that since MACRIS assumes an elastic homogeneous reservoir and burden, varying the Young's modulus would lead to an identical stress path and thus identical Coulomb stress rate history.
With MACRIS the Poisson's ratio is involved at two levels in the calculation of the stress changes induced by the poro-elastic deformation of the reservoir (Fig. 3). First, the Poisson's ratio scales the strength of the compacting sources (see eqs 1-3), that is the compacting grid cells around the Groningen faults. Secondly, the same Poisson's ratio is used when the induced elastic strains along the Groningen faults are converted into induced elastic stresses by applying the generalized Hooke's law (see eqs 4-6). Varying the Poisson's ratio in MACRIS, one can generate multiple stress solutions and thus test our hypothesis. To test the spatial heterogeneity in the Poisson's ratio, the two areas of interest, C-area and SW-area, are treated separately. The Poisson's ratio is kept spatially uniform at the scale of each area, but an ensemble of stress solutions is computed by varying the Poisson's ratio.
Figs 4 and 5 display the spatio-temporal evolution of the Coulomb stress changes along the Groningen faults of respectively the SWarea and C-area, assuming an identical Poisson's ratio of 0.2. The C-area presents a large number of faults with high Coulomb stress changes kicking off early in the history. Instead, in the SW-area, a single fault strand (oriented NW-SE) concentrates almost all the high Coulomb stress changes, and it is only late in the history that multiple faults start to experience high Coulomb stress changes. The contrast of fault orientation between both areas explains part of the difference in the spatio-temporal evolution of the Coulomb stress changes. Most of the faults of the C-area are optimally oriented; in contrast, the SW-area contains only one optimally oriented NW-SE fault strand. Fig. 6    patch (that is, the distance between two pillars) is roughly constant for all the faults of the Groningen field. Overall the cumulative distributions of Coulomb stress rate of fault patches of the C-area are much narrower than those of the SW-area. Three key observations can be drawn. First, the C-area contains a higher fault density underlined by a higher total number of fault patches. Second, the Carea presents a higher number of fault patches experiencing average to high Coulomb stress rates (between 0 and 0.2 MPa yr -1 ). Thirdly, the SW-area displays a higher number of fault patches experiencing very high Coulomb stress rates (>0.35 MPa yr -1 ). To summarize, the C-area is characterized by a higher density of fault patches with an average to high Coulomb stress rate kicking off early in the history. Instead, the SW-area is characterized by a lower density of fault patches, with a relatively large number of them experiencing very high Coulomb stress rates at a late stage during the field history.
Focusing on the C-area and now assuming an uniform but relatively small Poisson's ratio of 0.05, Figs 7 and 8 indicate, respectively, the spatio-temporal evolution of the Coulomb stress changes and the cumulative distribution of the Coulomb stress rates. As expected, for each time step, both the magnitude of stress changes and the number of fault patches experiencing high Coulomb stress rates are larger than in the previous scenarios with a Poisson's ratio of 0.2.

E N S E M B L E -S M O O T H E R F O R S E I S M I C I T Y DATA
The objective is now to develop an efficient data assimilation scheme in order to find the optimum set of posterior model parameters that give the best agreement between the observed seismicity rate and the computed rate. More specifically, the objective is not only to find the best set of model parameters: (1) to explain the observed temporal evolution (as performed e.g. by Dempsey & Suckale 2017or Heimisson et al. 2020 and (2) but also to match the spatial pattern (as e.g. Bourne et al. 2018). Before describing the details of the data assimilation procedure, it is important to define the earthquake catalogue used as data. This earthquake catalogue has been compiled by the Royal Netherlands Meteorological Institute (KNMI). KNMI has monitored seismicity in the Netherlands since at least 1986. From the dedicated borehole network deployed since 1995 by KNMI and constantly upgraded over the years, the earthquake detection capability improved: the magnitude of completeness kept decreasing from 1.5 in 1995 to around 0.5 today. For sake of direct comparison of our results with the ones of Candela et al. (2019), we restrict our analysis to all events with M L ≥ 1.0. As mentioned previously, constrained by the large uncertainty attached to the depth of each observed event, the seismicity catalogue is considered as a 2-D field. In order to compare the computed 2-D seismicity rates to the observed ones, both uncertainties in model (e.g. pore pressure distribution, fault orientation, stress calculation) and data (e.g. earthquake location) should be accounted for. We applied a Gaussian smoothing to the seismicity rate R D fields to incorporate uncertainties. In order to compute the smoothed seismicity rate fields, R D , for each R D field, the following integration over space is performed: where Gaussian operator is defined as: As a result, the five model parameters to be optimized are: (1) the Poisson's ratio used in the computation of the stress development with MACRIS and (2) the three parameters [A, r 0 ,Ṡ 0 ] of the Dieterich's seismicity rate theory governing the fault frictional response, and the characteristic length scale (standard deviation) of the Gaussian smoothing σ s .
It is important to note here that MACRIS is probably the only modelling approach that resolves 3-D poro-elastic stress changes along multiple faults with a metre-scale spatial resolution while honouring the full details of the reservoir and fault geometries. Finiteelement numerical computations could achieve this only through a tremendous meshing effort and running time. Here, it takes 5 hr on a standard PC for MACRIS to compute the depletion-induced stress changes along the 76 faults of the C-area with a metre-scale spatial resolution. Even if such achievement is exceptional, this run time is still not adapted for a traditional data assimilation scheme such as brute-force grid search or the Markov Chain Monte Carlo algorithm. A much more efficient data assimilation procedure, with a low number of forward simulations, needs to be developed for a robust calibration of model parameters.
An ensemble-based approach is deployed to update the prior model ensembles with the use of the data. The Ensemble Kalman Filter has proven to be very effective in weather forecasting, but also in reservoir engineering approaches (Evensen 2009). Every time that new data are collected, the Ensemble Kalman Filter procedure is applied, and the model parameters and predictions are updated. A progressively developing forecast mean and bandwidth is created. Another alternative, often used for subsidence inversion (Fokker et al. 2012Baù et al. 2015), is to use an Ensemble-Smoother where the new posterior ensemble is constructed in a single step (van Leeuwen & Evensen 1996). This is the data assimilation procedure deployed here for the first time over a complex case with real seismicity data. Another alternative is to apply multiple steps of data assimilation, so called Ensemble-Smoother MDA (Emerick & Reynolds 2013a). In this last approach, the update of the model parameters is more 'progressive', and in some cases it allows to prevent the emergence of an ensemble collapse where the posterior ensemble converges towards one unique solution (Emerick & Reynolds 2013a).
One outstanding challenge is the discrete nature of the seismic events. Indeed, ensemble-based approaches, both Ensemble Kalman Filter and Ensemble-Smoother, are designed for assimilating continuous data. As developed by Tarrahi & Jafarpour (2012) and Tarrahi  2015) for the Ensemble Kalman Filter and for synthetic scenarios and data, a Gaussian smoothing is applied here to the real discrete seismicity data before applying the Ensemble-Smoother in the same fashion as defined in eqs (14) and (15). Thus two distinct Gaussian kernels are applied here. The first one σ s to smooth the model outcomes and honour model uncertainties, and as mentioned previously considered as a model parameter to be optimized during the data assimilation step. The second Gaussian kernel must be applied to the seismicity data where the standard deviation of the kernel is fixed and is intended to map the uncertainty in event locations. Note here that, in practice, this approach would allow to handle: (1) different uncertainties for each event location and (2) different uncertainties between x and y directions. In this study, a unique, isotropic Gaussian kernel is selected.
The Ensemble-Smoother approach consists of an inversion scheme, for which the goal is to maximize an objective function of the form (Menke 1989;Tarantola 2005): where m and G(m) are, respectively, the optimized (posterior) vector of model parameters [ν, A, r 0 ,Ṡ 0 , σ s ] and model predictions.
Following this approach, the objective function is integrated in an inversion scheme seeking the solution for the vector m of model parameters that optimize the match with data d and with prior information m 0 . The Ensemble-Smoother explicitly includes our prior knowledge m 0 and the model covariance matrix C m (explained in the sequel of this section) with the term (m − m 0 ) T C −1 m (m − m 0 ). The optimal 'least-square' solution for one particular realization, assuming a linear inverse problem, is given by: For an ensemble-based estimate with a non-linear problem we define G M 0 as the result of the non-linear forward model working on all the members N e of the ensemble m 0 . Then, G M 0 is the ensemble of prior event density predictions covering the time window of interest. For an ensemble-based estimate, the Ensemble-Smoother then gives as updated model parameter ensemble: Eq. (18) corresponds to an ensemble of seismic event density data realizations created adding to the vector data d different random noise vectors ε that lie within the data uncertainty range. Following Tarrahi & Jafarpour (2012), the error of the event density data is considered as proportional to the value of the event density data and the variance at the kth grid block of the event density data ρ is computed as: where σ max and σ min are, respectively, the minimum and maximum variances of the event density data, and N d is the total number of grid block of the event density data. The member j of the perturbed observed event density at the kth grid block, can then be written as: with ε Two ensembles of forward simulations thus need to be run, the prior ensemble with the prior ensemble of model parameters M 0 , and the posterior ensemble with the updated posterior ensemble of model parametersM. A simple sensitivity analysis revealed that once the size of the prior ensemble is as large as 100-200 members, M is stable. Thus the total number of forward simulations can be kept as low to 100-200 (prior + posterior) in order to obtain a robust model conditioning. This is an excellent achievement considering that with traditional estimation algorithms (e.g. brute-force grid search or the Markov Chain Monte Carlo algorithm), several thousands of simulations would be required for a proper coverage of the prior model parameter space. Fig. 9 illustrates the pre-processing procedure we described and to be applied to both the modelled predictions and observations before an Ensemble-Smoother update.

F U L L G RO N I N G E N F I E L D W I T H O N E S I N G L E S T R E S S S C E N A R I O
Before running prior ensembles of MACRIS stress realizations varying the Poisson's ratio for both the C-area and SW-area (covered in the next Section 5), the full Groningen field is considered with a single MACRIS stress realization computed with a fixed Poisson's ratio of 0.2. In other words, we seek to optimize, for the time window from 1 January 1993 to 31 December 2017 , solely 4 model parameters: the three parameters [A, r 0 ,Ṡ 0 ] of Dieterich's seismicity rate theory governing the fault frictional response, and the standard deviation of the Gaussian smoothing σ s . Posterior ranges of model parameters (Fig. 10) for the full Groningen field are similar to the ones derived by Candela et al. (2019) for the C-area (see Table 2). More specifically, for the full Groningen field and considering an average effective normal stress of 12.5 MPa (derived from MACRIS) the mean characteristic relaxation time of  seismicity ( t a = Aσ n /Ṡ 0 ) is 108 yr. In Candela et al. (2019), for the C-area t a was estimated to be 87 yr and for the SW-area to be 6700 yr. The similarity between the t a derived from the whole Groningen field and the t a derived for the C-area is not surprising: because most of the induced events are located in the C-area, during the data assimilation step, the C-area dominates the update-procedure of the spatially uniform sets of model parameters of the full Groningen field.
The temporal evolution of the mean posterior event rate for the full Groningen field closely follows the one of the data (see Fig. 11). This match is solely controlled by the update of the three parameters [A, r 0 ,Ṡ 0 ] modulating the fault frictional response. Indeed, the update of the Gaussian smoother σ s solely controls the spatial pattern. Consequently, since the calibration of the model parameters is dominated by the C-area, where most of the data-events are located, a significant residue is observed between the spatial patterns of event density of the mean posterior and the data (Fig. 12). The temporal event rate is nothing more than the spatial integration of the event density. Thus in order to match the spatial pattern of the data event density, the mean posterior model overestimates the event density at the locations (including the SW-area) surrounding the C-area, and it underestimates the event density in the C-area.
To conclude for the full Groningen exercise, the Ensemble-Smoother performs well, although the spatial pattern of event density of the posterior clearly diverges from that of the data. This mismatch must be attributed to our forward modelling strategy, which is certainly missing complexities when the model parameters of both the induced stress development and the fault frictional response are assumed spatially uniform.

D UA L C O N T R I B U T I O N O F S T R E S S D E V E L O P M E N T A N D F R I C T I O N A L R E S P O N S E
Varying the Poisson's ratio, two separate prior ensembles of MACRIS stress realizations are now computed for both the C-area and SW-area. For the sake of comparison with the results obtained by Candela et al. (2019), the calibration time window is now from 1 January 1993 to 31 December 2015 . We seek to optimize both the parameters controlling the stress development and the ones controlling the fault frictional response (that is, the Dieterich's model parameters).
The resulting posterior distributions of the frictional parameters [A, r 0 ,Ṡ 0 ] (including the Gaussian smoother σ s ) are significantly different for both areas of interest (Fig. 13). Considering an average effective normal stress of 10 and 12.5 MPa (both derived from MACRIS) for, respectively, the C-area and SW-area, the mean characteristic relaxation times of seismicity are 90 and 5470 yr, respectively.
The posterior distributions of the Poisson's ratio ( Fig. 14 and Table 2) are also significantly different for both areas of interest. The mean posterior of the SW-area is very close to the a priori constant value of 0.2 used in Candela et al. (2019). Instead, the inverted mean posterior Poisson's ratio of the C-area is very low, close to 0.015.
For both areas of interest, the posterior ensembles of event rates (Fig. 15) follow very closely the long-period fluctuations of the observations. In both cases, the range of variability of the posteriors is sufficiently large, demonstrating that the collapse of the posterior ensembles towards one unique solution has been prevented (Emerick & Reynolds 2013a). This is particularly interesting considering that a rather drastic update procedure of both prior ensembles were imposed since only one step of data assimilation has been performed. Following a more progressive update procedure applying multiple steps of data assimilation is apparently not warranted in our case. The larger variance of the posterior ensemble event rate of the SW-area is a direct consequence of the relatively broader posterior distributions of its model parameters.
Figs 16 and 17 display the posterior event densities for the Carea and SW-area, respectively. Visual inspection reveals a persistent mismatch between model and observations. However, compared to the full Groningen inversion (with a constant Poisson's ratio of 0.2), the spatial patterns of the posterior event densities have improved.  This improvement can be visualized for the C-area where now the posterior event density does not systematically underestimate the observed event density over the entire surface area anymore. Fig. 18 demonstrates that the difference in the posterior model parameters between the C-area and SW-area has a significant impact on the modelled cumulative number of events. Using the mean posterior model parameters of the SW-area for the C-area, the cumulative number of events is significantly lower than both (i) the observations at the C-area and (ii) the modelled cumulative number of events for the C-area using the mean posterior model parameters of the C-area.

R E S O L U T I O N O F T H E P O I S S O N ' S R AT I O
Can our inversion procedure effectively constrain the Poisson's ratio? Can we demonstrate that a seismicity history cannot be equally explained by two different Poisson's ratios even after tuning the Dieterich's model parameters? Our objective here is to assess if our inversion procedure can discriminate two models computed with two different Poisson's ratios. A corollary is: can we estimate why our results indicate that a model with a lower Poisson's ratio can better explain the C-area seismicity rate history? We address these Figure 11. Comparison of the predicted seismicity histories (yearly event rate) with the data (dark) for the full Groningen. Cyan thin lines: prior ensemble with mean (thick blue line). Magenta thin lines: posterior ensemble with mean (thick red line).
questions for the C-area by comparing the quality-of-fit of the data by the optimum model obtained in this study and the one derived for a Poisson's ratio of 0.2 by Candela et al. (2019). We use two metrics to assess the quality-of-fit. Our first metric is the R M S E, defined as the root-mean-square of the residuals between data and the optimum model: where N d is the number of observations. Our second metric is the χ 2 method, which judges if the average of the squared residuals is of the order of the sum of the covariances of the data and model. The normalized χ 2 N d reads: where C d and C G(m) are, respectively, the covariances of data and model. Lower values of RMSE indicate a better fit. A χ 2 N d close to unity means that the model matches the data at a level consistent with the error covariance of the data. For the C-area, both metrics indicate a better fit of the model derived here (R M S E = 5.4, χ 2 N d = 0.8) than of the one derived in Candela et al. (2019) (R M S E = 7.9, χ 2 N d = 2.8). Our inversion procedure can thus effectively discriminate predictions with two different Poisson's ratios, and using a lower Poisson's ratio for the C-area leads to a better match of the observations.
To further examine the resolution of the Poisson's ratio and to gain confidence in the results obtained by our inversion procedure, we consider an 'ideal' case of laterally extensive depleting reservoir undergoing uniaxial compaction. Appendix A3 demonstrates that a seismicity rate history obtained with one stressing history (that is generated by one Poisson's ratio) cannot be reproduced by a different Poisson's ratio adjusting the Dieterich's parameters. Both the performance metrics (R M S E and χ 2 N d ) of the real case and our analytical exercise (see Appendix A3), reveal that the Poisson's ratio (controlling the stress development during the Groningen depletion history) and the Dieterich's model parameters (controlling the frictional response of the Groningen faults) can be well constrained.

'Apparent' background stressing and seismicity rates
Our assimilation procedure of observed events started in 1993 when a significant earthquake activity started to be recorded. Consequently, Dieterich's seismicity rate eq. (4) was solved assuming an initial condition at steady state, that is R (0) = 1 in 1993. However, the start of the human-induced perturbation of the Groningen field goes back to 1968, and hence it is most likely that the background activity was not at steady state in 1993. The background stressing rate and seismicity rate inferred for 1993 should not be interpreted as real steady state background values (before the start of human-induced perturbations by gas production), but should be understood as 'apparent background values', which are actually the reference values at the initial time of the analysis, here 1993. Other alternative strategies could have been deployed to start the analysis at the onset of gas production in 1968 and to honour the fact that the Groningen faults were probably far from failure at this time. One such approach is to add an activation threshold as suggested by Zhai et al. (2019) to model Oklahoma induced seismicity, or more recently by Heimisson et al. (2020) to handle the Groningen dormant faults. Heimisson et al. (2020) reformulated the seismicity rate theory of Dieterich (1994) to take into account that seismic sources can be initially far from instability. Their proposed solution is the following (eq. 1 of Heimisson et al. 2020): where t b is the time when the Coulomb stress threshold S c was reached, and σ 0 is the effective normal stress at t = 0, when field production started. Setting t = t − t b and noting that S(t) = S (t ) − S c is the Coulomb stress change after time t b , eq. (24) has the same form as the 'integral form' of the original seismicity rate equation of Dieterich (1994), but integrated starting at time t b : Appendix A4 gives the detailed derivation of eq. (25) starting from eq. (8).
Future work should leverage our fine-scale Coulomb stress model and compare our current approach starting the analysis in 1993 with the one of Heimisson et al. (2020) starting the analysis at the onset of gas extraction. However, the objective of this paragraph is to show the validity of our approach. Solving the original differential eq. (8) starting in 1993, we assumed that the Groningen faults were less than one stress drop away from failure, that is, critically stressed in 1993. This assumption is valid if the Coulomb stress threshold S c was reached before 1993 at the C-area and SW-area; that is, if t b ≤ 1993. The average pressure-drop from the onset of gas extraction to 1993 is roughly twice the average pressure-drop from 1993 to 2016, that is, at least 15 MPa (NAM 2021). Following the analytical solution (A15) for a laterally extensive depleting reservoir undergoing uniaxial compaction, the effective normal stress at the start of gas extraction σ 0 was ∼5 MPa, that is ∼5 MPa lower than the average value of ∼10 MPa in 1993. Since A is assumed to be a time-independent fault constitutive parameter, one can use the values derived in our analysis at 1993 (see Table 2) to calculate the Coulomb stress threshold S c , which is expected to be in the Downloaded from https://academic.oup.com/gji/article/229/2/1282/6459727 by guest on 06 April 2022 Figure 12. Comparison of the spatial pattern of the models and data event densities (events per km 2 ) for the full Groningen. Top panels: model, middle: data, bottom panels: residue. The residual map is the difference between data and model. range 1 − 10 Aσ 0 (eq. A9 of Heimisson et al. 2020), that is 1-17.5 MPa for both the C-area and SW-area. On the other hand, since the average pressure-drop from the onset of gas extraction to 1993 is roughly twice the average pressure-drop from 1993 to 2016, to first order, the induced Coulomb stress changes from the onset of gas extraction to 1993 are also expected to be twice those derived from 1993 to 2016 (Figs 5 and 7), that is ∼20 MPa. The induced Coulomb stress changes caused by pressure depletion up to 1993 of ∼20 MPa are thus higher than the range of expected Coulomb stress thresholds S c = 1 − 17.5 MPa. Therefore, the threshold in Coulomb stress was probably reached before 1993, and our assumption that the Groningen faults were within one stress drop away from failure in 1993 is plausible. This assumption is also supported by the fact that at least ∼10 events of a magnitude higher than 1.0 were recorded before 1993 at Groningen (Dost et al. 2018), even if the spatial coverage of the Groningen seismic network was still suboptimal then.
We thus demonstrated that in 1993 the Coulomb stress threshold was very likely already reached at the C-area and SW-area. Furthermore, performing our analysis separately for the C-area and SW-area we intend to honour the difference in background seismicity rate r 0 at 1993 between the C-area and SW-area. The original differential eq. (8) are thus solved starting at 1993, with initial conditions r 0 and σ 0 taken at that time and separately for both the C-area and SW-area. Equivalently, the integral form (25) can be integrated starting at 1993 (the start time is encapsulated in the definition of t = 0). In our analysis, we make use of the Coulomb stress rate when solving the ordinary differential eq. (10). Thus, our analysis involves only Coulomb stress changes after 1993, relative to the stress values at 1993. In particular, as long as the Coulomb stress threshold is reached before 1993 and that the background seismicity rate can be considered spatially uniform over the C-area and SWarea, our modelling (and our inferred values of Dieterich's model parameters and Poisson's ratio) is not affected by whether or not there is a threshold reached at t b .
It is important to mention again that in our approach, the modelled seismicity rate accounts for the changes in effective normal stress during the depletion history. Instead, in eq. (24) proposed by Heimisson et al. (2020) and in eq. (25) these changes in effective normal stress are not honoured. Indeed, the effective normal stress σ 0 (used in the calculation of t a and the Coulomb stress) is assumed constant and equal to the initial effective normal stress at t = 0; and the friction coefficient (used in the calculation of the Coulomb stress) is also assumed constant and equal to the ratio between the initial shear and effective normal stress.

Potential missing modelling ingredients for the stress development
The lower Poisson's ratio of the C-area might be a direct consequence of its relatively higher reservoir porosity compared to the SW-area (Fig. 2). Still it remains to be understood why our inversion procedure favors such very low mean posterior Poisson's ratio (0.015) for the C-area. Poisson's ratio values inferred from smallscale laboratory experiments (e.g. Hol et al. 2018;see Buijze 2020 for an exhaustive summary) performed on Groningen core material are around 0.25, and range from 0.07 to 0.34. Thus, the very low Poisson's ratio inverted here seems to be very unlikely. This tends to indicate that our approach might be missing modelling components which will be discussed in this section. Clearly further efforts will have to be devoted to identify these missing modelling components. Nevertheless, we consider that our current approach captures the correct stress development even if an unreasonably low Poisson's ratio is needed, and advances the core goal of the present work to disentangle the relative contributions of (1) the induced stress development and (2) the fault frictional response. Interestingly, results of the inversion performed by Dempsey & Suckale (2017) with a different modelling approach for induced stresses and seismicity at Groningen also favored a very low Poisson's ratio.  It is also important to mention that the uncertainties in the initial stress state were not fully honoured in our current modelling approach. Indeed, the horizontal to vertical effective stress ratio, the initial minimum to maximum horizontal stress ratio, and the direction of the maximum horizontal stress were all considered as constants (Table 1). In a future implementation, these uncertainties could be honoured. However, even if the initial stress state were spatially variable at the Groningen scale, this variability would not influence the Coulomb stress rate and thus would not explain the very low Poisson's ratio inferred for the C-area.

Stress interactions between events
For a lower Poisson's ratio, the seismicity data requires a higher stressing rate. One missing ingredient that could lead to higher stressing rates without the need of a very low Poisson's ratio is the elastic stress transfers between induced events. Instead of tackling the difficult challenge of explicitly modelling the additional stress changes induced by each earthquake, we here focus on the statistical removal from the observed catalogue of the events triggered by these stress transfers. indicates that the C-area is more prone to earthquake triggering than the SW-area: the aftershock proportion is 15 per cent for the C-area and only 4 per cent for the SW-area. The relatively higher percentage of triggered aftershocks might explain why only the Carea is characterized by a very low mean posterior Poisson's ratio. This hypothesis can be further tested by comparing the modelled event rate of the C-area using a Poisson's ratio of 0.2 against the declustered catalogue (Fig. 19). As expected, using a Poisson's ratio of 0.2, the modelled cumulative number of events is significantly lower relatively to the one derived using the mean optimum inverted value of 0.015. More importantly, using a Poisson's ratio of 0.2, the modelled cumulative number of events is lower than the cumulative number of main shocks of the declustered catalogue. The difference in optimum Poisson's ratio between the C-area and SW-area should thus persist even after accounting for the elastic stress transfers.

Plastic compaction
An additional missing modelling ingredient is the plastic compaction of the reservoir. Indeed, recent laboratory experiments on core samples of the Groningen reservoir rock show that a large part Downloaded from https://academic.oup.com/gji/article/229/2/1282/6459727 by guest on 06 April 2022 Figure 16. Comparison of the spatial pattern of the models and data event density (events per km 2 ) for the C-area. For each year, the left column corresponds to the event densities and residues displayed in Fig. 12 for the inversion of the full Groningen with a constant Poisson's ratio of 0.2. For each year, top row: model, middle row: data, bottom row: residue. The residual map is the difference between data and model. Downloaded from https://academic.oup.com/gji/article/229/2/1282/6459727 by guest on 06 April 2022 Figure 17. Comparison of the spatial pattern of the models and data event density (events per km 2 ) for the SW-area. For each year, the left column corresponds to the event densities and residues displayed in Fig. 12 for the inversion of the full Groningen with a constant Poisson's ratio of 0.2. For each year, top row: model, middle row: data, bottom row: residue. The residual map is the difference between data and model.  (up to 50 per cent) of the deformation occurring during laboratorysimulated depletion is inelastic (Hol et al. 2015(Hol et al. , 2018Pijnenburg et al. 2018Pijnenburg et al. , 2019. Can the relatively lower optimum Poisson's ratio of the C-area compared to the one of the SW-area be explained by plastic compaction? In other words, can the higher propensity to plastic compaction of the C-area explain its apparently lower Poisson's ratio? Multiple factors can potentially control the relative proportion between elastic and plastic deformation, for example rock porosity, clay content, initial mean effective stress directly linked to the reservoir depth (Buijze 2020). These factors are contrasted between the C-area and SW-area; for example rock porosity (Fig. 2) and reservoir depth are both higher for the C-area. Interestingly Buijze (2020) indicates that the primary response of the Groningen stress path can be approximated with a linear elastic behavior. In other words, using an apparent Poisson's ratio, the laboratory-observed elasto-plastic stress path of the Groningen reservoir rocks can be mimicked. One can thus postulate that the difference in optimum Poisson's ratio between both areas of interest could be an apparent result of a difference in the relative proportion of elastic and plastic deformation.

Vertical elastic contrast between reservoir and burden
Another elastic parameter which could have been considered in addition to the Poisson's ratio is the Young's modulus. However, as mentioned previously, the MACRIS calculation (as well as any poroelastic solutions honouring the fault offsets and depletion heterogeneities) would lead to an identical stress path and thus Coulomb stress rate history when varying the Young's modulus (Fjaer et al. 2008;van Wees et al. 2018). Indeed, when elastic properties are considered identical between the reservoir and its surrounding, varying the Young's modulus has no effect on the stress path. However, refining the current MACRIS calculation in order to handle the potential contrast in Young's modulus between the reservoir and its surrounding, would lead to different stress paths. Using a finite-element approach van Wees et al. (2018) already showed that for a synthetic gas field mimicking the Groningen field, a reservoir with a relatively lower Young's modulus compared to its surrounding leads to a steeper stress path than in an elastically homogeneous reservoir-burden scenario. Van Eijs et al. (2006), screening multiple hydrocarbon fields in the Netherlands, also showed that the stiffness contrast between seal and reservoir rock is highly correlated with the occurrence of induced earthquakes.
As introduced previously for the elastically homogeneous reservoir-burden scenario, considering a laterally extensive depleting reservoir undergoing uniaxial compaction, the Coulomb stress change scales linearly with the stress path and thus the Poisson's ratio (see eq. 13). The lowering of the Poisson's ratio from 0.2 to 0.015 corresponds to an increase of the stress path from 0.75 to ∼1. This last 25 per cent increase of the stress path should thus lead to an equal increase of 25 per cent of the Coulomb stress rate. Interestingly, the Coulomb stress rates obtained with MACRIS for all the 'fault patches' of the C-area using a Poisson's ratio of 0.015 are also ∼25 per cent higher than those obtained with a Poisson's ratio of 0.2 (Fig. 20).
The objective is now to identify if a realistic contrast of Young's modulus between reservoir and its surrounding could lead to such a 25 per cent increase of Coulomb stress rate relatively to the elastically homogeneous reservoir-burden scenario. While keeping constant the Poisson's ratio (ν = 0.2), the relatively higher porosity of the reservoir layer at the C-area (Fig. 2) would lead to a relatively lower reservoir Young's modulus and thus a relatively higher contrast with the stiff surrounding. Considering a synthetic single-fault reservoir model (similar to van Wees et al. (2018) and mimicking the Groningen field), finite-element results screening a range of realistic contrast of Young's modulus between reservoir and burden demonstrate that a 25 per cent increase of Coulomb stress rate relatively to the elastically homogeneous reservoir-burden scenario can be reached (see Appendix A5). Indeed, keeping an identical Poisson's ratio of 0.2 for both reservoir and burden but using a Young's modulus between 10 and 20 GPa for the reservoir and 50 GPa for the burden, can lead to 25 per cent increase of the Coulomb stress change relatively to the elastically homogeneous reservoir-burden scenario.
Two independent modelling strategies, one honouring the effect of the contrast in Young's modulus between the reservoir and its surrounding, and our approach with MACRIS honouring the Poisson's ratio effect, thus point toward the same: a relatively higher Coulomb stress rate for the C-area. The reality is probably in between, and the relatively higher Coulomb stress rate for the C-area might be a consequence of both: (1) Young's modulus contrast and (2) low reservoir Poisson's ratio; both effects being associated to the relatively higher porosity of the C-area.

Concluding remarks
One key advantage of our forward modelling strategy (coupling MACRIS with Dieterich's seismicity rate theory) is that the effect on the seismicity predictions of each modelling component (stress development and frictional response) can be individualized and thus tracked down. Indeed, our prime achievement is the disentanglement of the relative contributions of each physical process controlling the occurrence of induced seismicity at the Groningen gas field. Without deploying an efficient Ensemble-Smoother formulation, which only requires a small prior ensemble to constrain the posterior distribution of each model parameter, this disentangling exercise would not be feasible.
Our integrated workflow (coupling a physics-based forward modelling and an efficient data assimilation procedure) demonstrates that even after honouring the potential spatial heterogeneity in stress development across the Groningen gas field, the spatial variability of the observed induced seismicity at Groningen requires spatial heterogeneity in fault frictional responses. A corollary is that any modelling strategy seeking to predict seismicity in the Groningen field must explicitly honour spatial heterogeneities in both (1) stress development and (2) fault frictional response.

A C K N O W L E D G E M E N T S
The first author dedicates this publication to the memory of Bogdan Orlic. Kay Koster is generously thanked for his help to make the visual on the geographical location of the Groningen gas field. We also gratefully acknowledge the anonymous reviewers whose suggestions improved our initial contribution.

DATA AVA I L A B I L I T Y
The earthquake catalogue used for this study can be found at the KNMI website at www.knmi.nl. The 4-D yearly induced stress changes and pore pressures at faults generated during this study are available from the corresponding author on reasonable request.

A P P E N D I X A1 Comparison MACRIS versus finite-element solution
The MACRIS stress solution has been already extensively compared against finite element code in van Wees et al. (2018van Wees et al. ( , 2019. The objective here is to follow the same benchmarking exercise for the specific case of Groningen as followed by van Wees et al. (2018). We use the commercial finite-element (FE) package DIANA FEA [http://dianaf ea.com, for more details about the FE modelling procedure the reader is referred to Orlic & Wassing (2013)]. The selected, and built for this exercise, plane-strain elastic FE model of a depleting gas reservoir includes two reservoir compartments, separated by a fault with a throw of half the reservoir thickness ( Fig. A1 and Table A1). This specific fault offset is one of the key ingredients of the induced stress developed along the Groningen faults (Van den Bogert 2015; Buijze 2020). The reservoir depth and thickness of the numerical model, as well as all the underlying parametrization (including the elastic properties) are representative of the inputs used for MACRIS when running the real Groningen case. The final spatially uniform depletion throughout the reservoir is 25 MPa, again mimicking the present depletion at Groningen and modelled via the reservoir flow simulation used as input for MACRIS. For MACRIS, we use a 3-D representation of the single-fault reservoir with a square geometry in map view, and with a lateral x, y dimension of 60 km (Fig. A2). This very large lateral dimension, mimicking an infinite lateral extension of the reservoir, is required to match the plane strain boundary condition of the 2-D FE model. The 3-D MACRIS model used same parametrization for initial stress conditions and elastic properties as the 2-D FE model (see Table A1), and the same spatially uniform depletion throughout the reservoir of 25 MPa.
For the comparison of both solutions, the MACRIS-based and the FE-based, the stresses along the centre-pillar of the 3-D MACRIS model (Fig. A2) are compared against the ones of the 2-D FE model. Fig. A3 demonstrates that the effective normal, shear and Coulomb stress calculated with MACRIS are well superposed on the ones of the FE model. The stress peaks at the reservoir edges are identically captured by both approaches.

A2 From Dieterich (1994) to Segall & Lu (2015)
This section presents the detailed derivation from the original equation of Dieterich (1994) [eq. (8) in the current manuscript] to the formulation proposed by Segall & Lu (2015) [eq. (10) in the current manuscript]. Dieterich (1994) initially proposed: where θ is a state parameter, A is a frictional parameter, σ n is the effective normal stress andṠ is the Coulomb stressing rate. The current seismicity rate is defined as: where r 0 is the background seismicity rate andṠ 0 is the background stressing rate. Segall & Lu (2015) derived an alternative formulation of the Dieterich's model to eliminate the state parameter θ. They defined a normalized seismicity rate as: Using (A3) and (A2), leads to where t a = Aσ n /Ṡ 0 . It follows: That is the eq. (10) of this current manuscript as proposed by Segall & Lu (2015).

A3 Resolution of the Poisson's ratio
The objective of this section is to demonstrate that a seismicity rate history obtained with one stressing history (that is generated by one Poisson's ratio) cannot be reproduced by a different Poisson's ratio adjusting the Dieterich's parameters. In other words, the objective is to demonstrate that our modelling approach can effectively constrain the Poisson's ratio.
Re-writing eqs (9) and (10) in dimensional form it reads: Eq. (A10) tells us that the same seismicity history can be fitted by another Coulomb stressing rate history of different amplitude but same time-dependence (i.e. x times the originalṠ) by adjusting the Dieterich's model parameters Aσ n andṠ 0 r 0 (multiplying both by the constant x). In other words, if varying the Poisson's ratio would  solely modify the amplitude of the Coulomb stressing rate, the observed seismicity rate history could be equally fitted/explained by multiple Poisson's ratio after adjusting the Dieterich's model parameters Aσ n andṠ 0 r 0 . Therefore, our inversion procedure would not help to constrain the Poisson's ratio. However, this conclusion only holds if changing the Poisson's ratio only modifies the amplitude of the Coulomb stressing rate but not its shape (i.e.Ṡ 0 normalized by its max value), which is only fulfilled when the friction μ is kept constant in the calculation of the Coulomb stressing rate (eq. 11).
In our case, as implemented by Rubin & Ampuero (2007), the changes in the friction (defined as μ = τ/σ n ) during the induced stressing history (that is the depletion history in our case) are honoured. Therefore, when the Poisson's ratio varies, both the amplitude and the time-dependent shape of the Coulomb stressing rate are modified. This can be demonstrated by considering the analytical solution for a laterally extensive depleting reservoir undergoing uniaxial compaction (e.g. Fjaer et al. 2008) where the changes in effective normal stress and shear stress induced by a depletion history P(t) are, respectively: and where ϕ = 90 − fault dip. Assuming initial conditions as τ init = 5 MPa and σ init = 13 MPa, and a depletion history representative of the Groningen gas field, Fig. A4 presents the Coulomb stressing rate histories and effective normal stress histories considering two different Poisson's ratios. Aσ n in eq. (A10) cannot be interpreted as a constant, since in our approach the temporal changes in effective normal stress are honoured when solving the ordinary differential eq. (10). Modifying the Poisson's ratio affects thus substantially the shape of both the Coulomb stressing rate history and the effective normal stress history (see Fig. A4). Our inversion procedure should thus effectively constrain the Poisson's ratio since an observed seismicity rate history can only be reproduced by one unique modelled seismicity rate history generated by one unique Poisson's ratio. In other words, a seismicity rate history obtained with one stressing history (that is generated by one Poisson's ratio) cannot be reproduced by a different Poisson's ratio adjusting the Dieterich's parameters.
To discuss the resolution of the Poisson's ratio, one can rearrange eq. (A10) as: a X (t, ν) + bY (t, ν) = 1,  ratio ν, the optimal values of a and b can be readily obtained by standard least-squares regression of eq. (A13). The objective function is defined as: where m = [a, b] , G is a matrix [X, Y ] and d is a vector of ones [1]. The least-squares solution is: To account for the covariances of X and Y , an orthogonal regression solver can be used. X and Y are combinations of model predictions (σ n andṠ) and synthetic data (Ṙ D and R D will be generated by the model with true/known model parameters). In order to identify the resolution of the Poisson's ratio, one can consider, as previously for Fig. A4, the analytical solution for the development of stresses (see eqs A11 and A12) induced along a fault embedded in a laterally extensive depleting reservoir undergoing uniaxial compaction. In our analytical example, the synthetic data (Ṙ D and R D ) are generated with a known set of parameters: ν = 0.05, a = 0.22 and b = 0.31. Note that the Dieterich's model parameters are similar to the inverted mean optimums for the C-area, and the Poisson's ratio of 0.05 is considered as representative of the low inverted mean optimum value of 0.015 for the C-area (see Table 2). For each member of the prior range of Poisson's ratio, the least-squares regression problem is solved and both the optimal values of a and b, and the two metrics for the quality-of-fit (eqs 22 and 23) are obtained. Fig. A5 demonstrates that, according to both metrics, the best model is effectively the one computed with a Poisson's ratio of 0.05 as the one of the data. In addition, the minimum values of both metrics (that is the best model) is rather sharp and well defined, meaning that small deviations from the optimum Poisson's ratio of 0.05 already result in suboptimal match between model and data. This last observation demonstrates that the Poisson's ratio can be well resolved following our inversion procedure. In other words, the seismicity history of the data can only be explained by one unique Poisson's ratio. Fig. A6 indicates that suboptimal models generated with Poisson's ratios higher than the one used for the data (that is, 0.05), are characterized by lower Dieterich model parameters a and b relatively to the ones used for the data (that is respectively 0.22 and Figure A4. Effective normal stress and Coulomb stress rate histories for two different Poisson's ratio (blue curves: ν = 0.2; cyan curves: ν = 0.05) computed with the analytical solution of a fault (with zero offset and a dip of 70 • ) embedded in a laterally extensive reservoir undergoing uniaxial compaction (Fjaer et al. 2008). Top left-hand panel: depletion history representative for Groningen. Top right: effective normal stress. Bottom left: Coulomb stress rate. Bottom right-hand panel: Normalized Coulomb stress rate, that is the Coulomb stress divided by its max value. 0.31). This same trend derived for the synthetic analytical case is also observed using MACRIS and the real seismicity catalog. The optimum Dieterich's parameters a and b derived in this study with the optimum Poisson's ratio of 0.015 are higher relatively to the ones obtained by Candela et al. (2019) where a fixed Poisson's ratio of 0.2 was used (see Table 2).
Our analytical exercise can be considered as representative of the behavior of each fault patch constituting the C-area. Since the modelled seismicity rate for the C-area results from the integration of the local seismicity rate for each fault patch, one can also expect the Poisson's ratio to be well resolved for the real case C-area as indicated by the posterior distribution of the Ensemble-Smoother procedure (Fig. 14). Both the performance metrics (R M S E and χ 2 N d ) of the real case and our analytical exercise, reveal that the Poisson's ratio (controlling the stress development during the Groningen depletion history) and the Dieterich's model parameters (controlling Downloaded from https://academic.oup.com/gji/article/229/2/1282/6459727 by guest on 06 April 2022 Figure A7. Coulomb stresses for a range of Young's modulus contrast between reservoir-burden and obtained with the FE single-fault reservoir model presented in Fig. A1 The right-figure is a magnified section of the left-figure. For the right-figure, the depth extent of both reservoir compartments is indicated on both sides of the fault in grey. The Poisson's ratio of both reservoir and burden are kept constant equal to 0.2. The Young's modulus of the burden remains constant equal to 50 GPa while the Young's modulus of the reservoir Er varies from 50 GPa (homogeneous reservoir-burden scenario) to 10 GPa (strong Young's modulus contrast). The dashed black curves correspond to the homogeneous reservoir-burden scenario (with a Young's modulus of 50 GPa) amplified of 25 per cent.