A framework for subsurface monitoring by integrating reservoir simulation with time-lapse seismic surveys

Reservoir simulations for subsurface processes play an important role in successful deployment of geoscience applications such as geothermal energy extraction and geo-storage of fluids. These simulations provide time-lapse dynamics of the coupled poromechanical processes within the reservoir and its over-, under-, and side-burden environments. For more reliable operations, it is crucial to connect these reservoir simulation results with the seismic surveys (i.e., observation data). However, despite being crucial, such integration is challenging due to the fact that the reservoir dynamics alters the seismic parameters. In this work, a coupled reservoir simulation and time-lapse seismic methodology is developed for multiphase flow operations in subsurface reservoirs. To this end, a poromechanical simulator is designed for multiphase flow and connected to a forward seismic modeller. This simulator is then used to assess a novel methodology of seismic monitoring by isolating the reservoir signal from the entire reflection response. This methodology is shown to be able to track the development of the fluid front over time, even in the presence of a highly reflective overburden with strong time-lapse variations. These results suggest that the proposed methodology can contribute to a better understanding of fluid flow in the subsurface. Ultimately, this will lead to improved monitoring of reservoirs for underground energy storage or production.


Introduction
Understanding fluid flow in subsurface reservoirs is crucial to predict underground processes related to the energy transition, such as geothermal energy [1], temporary storage of green gasses like hydrogen [2], and long-term storage of greenhouse gasses like CO 2 [3].Reservoir simulations allow us to accurately predict fluid flow inside a reservoir, based on a combination of geological, geophysical and borehole data [4,5].Geophysical methods, such as seismic monitoring, are able to observe time-lapse changes of dynamic properties like pressure and fluid saturation, everywhere in a three-dimensional subsurface.Seismic monitoring relies on the fact that changes in the reservoir will translate into changes in the seismic reflection response.The fluid flow inside the reservoir can then be imaged by comparing a baseline seismic survey with a monitor survey, recorded over the same location at a later point in time [6,7,8,9].Feasibility studies aim to assess the seismic detectability of fluid movement inside a hydrocarbon reservoir [10] or migration of injected CO 2 for CCS projects [11,12].These types of studies rely on reservoir simulations to predict the movement of the fluids in the reservoir.Although this methodology provides accurate estimates of the time-lapse changes inside the reservoir, it does not predict geomechanical changes in the overburden.However, these changes can have large effects on the repeatability of time-lapse experiments, as overburden time-shifts might be mistaken for changes inside the reservoir [13].Generally, an independent geomechanical model is used to compute the time-lapse changes in the layers above the reservoir [14,15].Recently, multiphase poromechanical models were introduced as an all-in-one solution to link fluid flow, transport and deformation in the subsurface [16].Traditionally, these models are used to predict induced seismicity due to fluid injection in the subsurface [17,18,19].Additionally, poromechanical simulations can, in theory, also be used to model both time-lapse changes inside the reservoir and overburden at once for seismic monitoring applications.In addition to time-lapse overburden effects, static overburden effects can also obstruct the reservoir signal in the baseline and monitor seismic surveys, due to the presence of highly reflective layers in the overburden.Both the static and dynamic overburden effects can be accounted for by isolating the reservoir response [20,21].This isolation is based on the 3D Marchenko equations that describe all orders of multiple scattering inside the medium [22,23,24].After this Marchenko-based isolation is applied to the seismic data, the reflections related to the reservoir are clearly visible in the seismic response.Next, the primary reflection from the top of the reservoir is used as a reference event that contains all the delays of the overburden.This reference event is then combined with events originating from the reservoir's base to retrieve time-lapse differences that are solely dependent on the changes inside the reservoir [21].In this work, a poromechanical simulator is proposed to model time-lapse changes in density and compressional velocity due to fluid injection in a subsurface reservoir.Since multiphase fluid flow as well as geomechanics are included in the formulation, the Cross-correlating these primaries cancels their common path, hence only the timelag between the two reflections remains (i.e. the traveltime through the reservoir).changes in the overburden and reservoir are modelled all at once.Next, the velocities and densities are computed at a number of time-steps during the simulation, which are used to model the seismic response for the seismic baseline and different monitor studies.Finally, time-lapse changes are retrieved.These changes are then independently assessed both before and after isolation of the reservoir response from the total seismic response (i.e. the response of the overburden, reservoir and underburden).In the next section, we will first discuss the governing equations behind the poromechanical model, the connection with seismic parameters and the retrieval of time-lapse traveltime differences.Secondly, the methodology is tested on a simple as well as a complex model.To conclude, we discuss the results and possible future improvements and extensions to the method.

Methods
This section discusses the background on how time-lapse changes can be extracted from a modelled reservoir.The constitutive equations related to poromechanics are first reviewed, then these equations are related to seismic properties, which can be used to model the seismic response at different times in the simulation.These responses are compared to one another to find seismic time-lapse traveltime differences between the different surveys.

Multiphase poromechanics
The geomechanical changes in an isotropic subsurface are represented by the conservation of momentum (Equation 1 [25]), and conservation of mass describes flow of immiscible fluids through a reservoir (Equation 2 and 3 [26]).This gives the following system of equations: Here, m, ρ, v and q, are the fluid mass per unit volume, density, velocity and sinks/sources, respectively.The subscripts α and β denote two different fluid phases.Furthermore, t denotes time, σ σ σ is the stress tensor, and f stands for the body forces acting on the system.Next, the stress tensor in Equation 1 is connected to the changes in fluid pressure and displacement according to Biot's theory of poroelasticity [27].Moreover, the mass per unit volume of each phase is related to it's saturation, density and the porosity (i.e.m = φ Sρ).Finally, Darcy's law is used to write the fluid velocity in terms of phase mobility and pore pressure.After applying all these conditions, the system of equations 1-3 now reads [16]: In Equation 4b denotes Biot's coefficient, ∇ s u = 0.5(∇u + ∇u T ) is the symmetric gradient operator operating on displacement u, I is a unit tensor and C dr the rank-4 drained elasticity tensor, which for isotropic linear elastic material reads: Here, δ i j is Kronecker's delta, λ Lamé's first parameter, and G Lamé's second parameter or the shear modulus.The subscript dr denotes that the elastic moduli are drained.p f in equations 4-6 represents the pore pressure, which in the absence of capillary forces is the same for each phase.In Equation 5 and 6 S indicates the saturation, for two-phase flow S β = 1 − S α .Moreover, λ λ λ depicts the phase mobility, which is equal to the rock permeability times the relative permeability over the viscosity (Kk r /µ).Lastly, φ is the porosity, which differss from a reference φ 0 due to the fluid pressure and volumetric strain ε v = trace(∇ s u) as [28]: with drained bulk modulus K dr = λ dr + (2/3)G dr .

Seismic parameters via fluid substitution
After the dynamic fluid and geomechanic quantities have been computed by the poromechanical simulation, they have to be converted into seismic parameters, namely density and compressional wave velocity.Note that, in this study, only compressional waves are considered for the forward seismic modelling, even though retrieving the shear wave velocity is trivial once all elastic parameters are calculated.This is due to current limitations of the Marchenko-based isolation of the reservoir response, as the Marchenko equations are not straightforwardly applied to elastic theory [29], but extensions are under investigation [30,31].
The saturated density can be calculated using the fluid saturation and density as well as the porosity and rock density: Next, the compressional wave velocity c p is computed using the elastic moduli K and G as well as density ρ: where the subscript sat denotes a saturated medium.Gassmann's equation describes how the saturated bulk and shear moduli can be found [32,33]: Again, K dr and G dr are the drained bulk and shear modulus, respectively [34].K 0 is the bulk modulus of the minerals of the rock that can be experimentally determined [e.g.35,36].K fl is effective bulk modulus of the fluid that can, for example, be calculated using the Reuss method for uniform saturation [37]: where κ is the compressibility of the fluid, which is equal to the inverse of the bulk modulus (κ = 1/K), and can be derived from the pressure, volume and temperature of the fluid [e.g.38].Note that Equation 11 and 12 are only valid at low frequencies (< 100 Hz), which makes them ideal for field-scale experiments such as in this study [34].

Extracting time-lapse traveltime differences
First, the subsurface is divided in three units, overburden a, target zone b and underburden c.The reservoir is located in the target zone.Next, the velocity and density explained in the previous section are used to model the seismic reflection response R abc (x R , x S ,t) of the entire subsurface, specified by the subscript abc.The source and receiver coordinates at the surface are indicated with x S and x R , respectively.Additionally, t denotes the seismic recording time, which is different from the flow simulation time in Equation 5 and 6.The seismic recording time is typically in the order of seconds, whereas the flow simulation time is in the order of hours to days.Reflections in the overburden can interfere with the signal from the reservoir, which can prevent the accurate retrieval of time-lapse differences.Therefore, the reservoir response has to be isolated from the full response.This isolated response R b (x R , x S ,t) is free from over-and underburden reflections, which allows accurate retrieval of the time-lapse differences inside the reservoir [21].Details of this isolation are discussed in the supplementary material accompanying this paper.
Correlations are a popular method to extract time-lapse traveltime differences from seismic data [39,40].The traveltime differences ∆t can be found by cross-correlating the same signal between a baseline and monitor survey and taking the argument of the maximum of the correlation: In Equation 14, x 0 represents the zero-offset coordinate, where x S = x R .Moreover, C denotes the event to be correlated of the baseline survey, and C the same event in the monitor survey.If this event is simply a primary originating from a reflector below the reservoir, all the time-delays of the overburden will be present in the calculated traveltime differences.Instead a control reflection from above the reservoir can be used to first compute the timelag inside the reservoir [21].In Figure 1 this idea is systematically depicted; primary 1 does not travel through the reservoir, while primary 2 does.By cross-correlating these two primaries the timelag in the reservoir is computed, while all time-delays from the overburden are cancelled out.This is akin to the idea of seismic interferometry [41].Equation 15 describes how correlation of primaries 1 (P1) and 2 (P2) returns the timelag (C P1P2 ) between these two events: If the timelag in Equation 15 is computed for both the baseline and monitor survey, the retrieved correlations can be inserted into Equation 14 in order to acquire time-lapse traveltime differences that only encompass the changes in the reservoir layer.
A similar procedure can be applied to multiples, selected from the isolated target response R b .Since multiples have traveled through the reservoir layer multiple times, they are more sensitive to time-lapse changes in the reservoir [21].In the following, only primaries are considered.

Numerical examples
A customised fully implicit multiphase poromechanics simulator was designed for this work.The poromechanics part of this simulator was benchmarked on the 1D Terghazi and 2D Mandel problems [42,43].Furthermore, the two-phase fluid flow simulation was validated using the DARSim Matlab simulator [44].The simulator also includes fluid substitution to find saturated elastic parameters and produce subsurface density and velocity models for seismic modelling.Forward seismic modelling was achieved with an existing finite-difference modeller [45], and Marchenko-based isolation of the target response was performed with publicly available algorithms [21].In this section two models are considered.The first example is a simple piston-like flow in a homogeneous reservoir with simple overburden.Second, a more heterogeneous model is considered with a highly reflective overburden.

Case 1: Simple model
The subsurface model for this numerical experiment is shown in Figure 2, and the values for the properties of each layer and fluid can be found in Table 1.To start the simulation fluid α is injected on the left side of the reservoir (initially filled with fluid     .In (a) the injection and production wells are marked with the green and red lines, respectively.In (b) the arrows mark primary 1 (red) and primary 2 (blue).These primaries delineate the reservoir and will be used to extract traveltime differences.
β ) with a constant pressure of 50 MPa, and the production wells, in line with the right side of the reservoir, have a pressure of 5 MPa, while the initial pressure in the reservoir is equal to 10 MPa.The fluid flow is constrained to the reservoir, and roller boundary conditions (i.e.zero normal displacement) are imposed on the four edges of the model.The total simulation time is 600 days, and a seismic survey is modelled at every 100th day.The forward seismic model utilizes a zero-phase wavelet with a flat spectrum between 5 and 70 Hz, a time sampling of 4 ms, and 401 co-located sources and receivers at a spacing of 10 m.Traveltime differences between the baseline and monitor study at different times of the reservoir simulation for the simple model.The dashed lines represent the actual differences, whereas the solid lines are the differences extracted from the modelled seismic response.To highlight the importance of removing overburden effects, the dotted line shows the traveltime differences after 200 days when primary 2 (below the reservoir) is immediately correlated in the baseline and monitor study.Opposed to this, for the solid lines the overburden effects are first removed using primary 1 (above the reservoir) as reference, and only thereafter the baseline and monitor study are correlated to get the time-lapse changes inside the reservoir.
Figure 3 shows the time-lapse change in P-wave velocity after 200 days compared to the baseline (a), as well as the zero-offset reflectivity modelled at this time step (b).In Figure 3a, a decrease in velocity is noted above the injector wells (caused by the time-lapse response of the soft layer on top of the reservoir), whereas the velocity increases above the production wells and inside the reservoir.Inside the reservoir an increase in velocity is noticed due to fluid β being replaced with fluid α.
The changes in P-wave velocity above the reservoir are caused by the pressure change, that is, the increase in pressure due to injection leads to a decrease in velocity; vice versa the pressure decrease above the production wells causes an increase in velocity.Furthermore, primary 1 and primary 2, at the top and base of the reservoir are clearly visible in the seismic reflectivity section, as indicated by the red and blue arrows in Figure 3b.Due to this clear visibility, no Marchenko-based isolation is necessary for the simple model.Next, reflectivity is modelled for every 100th day between 0 and 600 days.The initial reflectivity at day 0 is used as baseline study and the subsequent reflection responses are considered monitor studies.
After the simulation and forward seismic modelling is finished, the correlation of P1 and P2 (i.e. the reflection inside the reservoir) is first computed for each of seismic study using Equation 15.To further improve the resolution of the results, the correlations are interpolated to 0.5 ms, by padding the data with zeros in the frequency domain.These correlations are then used to find the traveltime differences between the baseline and monitor studies (Equation 14).The results of this numerical experiment are shown in Figure 4.The dashed lines in the figure show the reference result based on the time-lapse changes in velocity.The solid lines are time-lapse traveltime differences retrieved with the proposed method.These lines clearly mark the fluid front advancing from left to right in the reservoir.The lines do not perfectly coincide with the reference result due to limitations in the spatial and temporal resolution of the seismic responses, which were measured with a time sampling of 4 ms and a receiver spacing of 5 m.Finally, to illustrate the effect of the overburden, the dotted line shows the traveltime differences computed solely from reflection P2 of the baseline and monitor data after 200 days.This is in contrast to the solid lines where reflections P1 and P2 are first correlated, and this correlation is then used to compute the traveltime changes between the baseline and monitor study.This means that the dotted lines do not retrieve the correlation of Equation 15, but rather insert reflector P2 of the baseline and monitor study directly into Equation 14.All overburden changes are, therefore, included in the dotted line, hence the location of the fluid front can no longer be accurately observed (i.e. it deviates from the reference result).

Case 2: Complex model
The second numerical experiment examines a more complex model shown on the right in Figure 2, with all layer and fluid parameters listed in Table 2.This model contains a highly reflective overburden, designed to produce strong multiple reflections that interfere with the reservoir response [46].Additionally, this also means that the overburden yields a strong response due to geomechanical changes.Furthermore, the reservoir is no longer rectangular, instead it has a wave-like structure with variable permeabilities between 2.5e − 14 and 1e − 11 m 2 .This permeability is pseudo random, generated using Perlin noise that allows for a somewhat coherent distribution [47].As shown in Figure 2, a single injector well (50 MPa) is located at the top of the reservoir at 1500 m lateral distance, alongside with a production well at 3000 m with a pressure of 5 MPa.The boundary conditions are the same as in the simple model.Every 50th day of simulation the reflection response is modelled, the total simulation time is 300 days.The seismic modelling uses the same parameters as in the first numerical example, except for the spectrum of the source wavelet, which is flat between 5 and 50 Hz.The evolution of the saturation is shown in the left column of Figure 5.This figure also displays the changes in velocity for 6 monitor time steps, which shows that the time-lapse changes outside the reservoir overpower the changes inside the reservoir.Again, the velocity decreases and increases above the injection and productions wells, respectively.Figure 6 displays the zero-offset seismic sections of the initial baseline and final monitor studies at 0 and 300 days, respectively.Figure 6a and 6b show that the two primary reflections (marked in red and blue) are obscured by overburden events.As seen in Figure 6c the time-lapse differences from the reservoir are masked by overburden effects.Consequently, it is beneficial to isolate the reservoir response, contrary to the first example.The results of this isolation are shown in 6d, 6e and 6f.The desired events (i.e.P1 and P2) are now revealed in the reflection response.Similarly, the isolation is applied to the remaining monitor surveys as well as the baseline survey.Next, the correlations are interpolated from 4 to 1 ms via the frequency domain for improved resolution.Subsequently, the correlations between P1 and P2 are computed, which will be used to find the time differences of the reservoir layer.
Once again, the time-lapse traveltime differences inside the reservoir are computed according to Equation 14; the results of this procedure are plotted with solid lines in Figure 7.As before, the dashed lines depict the reference result calculated based on the velocity changes in Figure 5.However, their behaviour is a lot more complex than the previous results.The solid, coloured lines represent the traveltimes difference computed from the correlations of P1 and P2 of the baseline and monitor surveys.Lastly, the gray lines show the results of computing the time-differences from the full medium responses (i.e.without isolation, Figure 6a and 6b), these results clearly deviate from the reference solution due to overburden effects appearing in the selection window of the primaries.Even though the computed differences somewhat agree with the reference solution, the match is significantly poorer than for the simple model.Especially around 2200 m the results strongly differ from the reference solution, as all 4 solid lines underestimate the reference indicated with the dashed lines.This underestimation is also observed in the gray lines in the background, and could, therefore, indicate that some remainder of overburden effects is still present in the isolated response.Nevertheless, the results after applying the reservoir isolation (i.e. the solid coloured lines) show a clear improvement in recovering the trend of the fluid movement compared to the solid gray lines.Further improvements may be possible using a migration technique to collapse the diffractions to their true location, thus improving the spatial resolution of the data.

Discussion
In the previous section it was shown that poromechanical modelling can add valuable insights to seismic reservoir monitoring, specifically because overburden changes are predicted together with fluid flow.In this section possible improvements upon both the poromechanical and the seismic part of the method are discussed.Firstly, poromechanical simulations are computationally expensive, which sets practical limitations on their use.Currently, the simulation takes around 8 hours, when using 8 CPUs on the Delft High Performance Computer [48].However, this time dramatically increases on a grid with finer discretization.One solution to this problem is to apply a preconditioner to improve convergence of the equations, allowing to still achieve high resolution simulations [16].Alternatively, a multiscale approach could be used in order to limit the size of the problem, thus speeding up the simulations [49,50,51].Additionally, the simulator could be improved by extending the formulations to include fractures.This is especially relevant as fractures can unexpectedly block the fluid or bypass impermeable zones.Another feature that is currently missing from the simulator is the ability to model cyclic storage, which is required to accurately monitor temporary storage of hydrogen or other gasses in the subsurface [52].It would also be interesting to include a third fluid to the simulator, in order to account for

Conclusion
This work developed a multiphase poromechanical simulator that is tightly coupled with a seismic survey simulator, through the update of the seismic parameters.This integrated simulator allows to instantly resolve time-lapse changes not only inside the reservoir but also in its overburden.The simulator was used to test the feasibility of a novel methodology to extract time-lapse travel time changes after Marchenko-based isolation of the reservoir response.This methodology solves the repeatability issue of time-lapse seismic surveys by identifying and utilizing a reference reflection above the reservoir.Future developments should focus on inverting the seismic time-lapse changes back to the dynamic reservoir properties, in order to close the loop between reservoir simulations and seismic monitoring.These results are a significant step to achieve higher resolution monitoring of subsurface reservoirs.A better understanding of fluid flow in these reservoirs will improve predictions of underground processes related to geothermal energy, subsurface storage of gasses like hydrogen, and sequestration of CO 2 .

Figure 1 .
Figure1.Principle of seismic interferometry: a reservoir is situated between two reflectors that reflect primaries 1 and 2. Cross-correlating these primaries cancels their common path, hence only the timelag between the two reflections remains (i.e. the traveltime through the reservoir).

Figure 2 .
Figure2.Subsurface model for the first (left) and second (right) numerical example.The reservoir for the first numerical example (in red) has a constant permeability of 1e − 13 m 2 , and fluid is injected at the left border of the reservoir and produced at the right border.For the second example the reservoir (in colour) has a variable permeability between 2.5e − 14 and 1e − 11 m 2 .The green and red triangles indicate the location of the injection and production well, respectively (which are now a single point source or sink).Other relevant properties of the numbered layers can be found in Table1for the first model, and Table2for the second model.

Figure 3 .
Figure 3. Time-lapse difference in P-wave velocity between 0 and 200 days of simulation (a) and zero-offset reflection response after 200 days of simulating fluid injection (b).In (a) the injection and production wells are marked with the green and red lines, respectively.In (b) the arrows mark primary 1 (red) and primary 2 (blue).These primaries delineate the reservoir and will be used to extract traveltime differences.

Figure 4 .
Figure 4. Traveltime differences between the baseline and monitor study at different times of the reservoir simulation for the simple model.The dashed lines represent the actual differences, whereas the solid lines are the differences extracted from the modelled seismic response.To highlight the importance of removing overburden effects, the dotted line shows the traveltime differences after 200 days when primary 2 (below the reservoir) is immediately correlated in the baseline and monitor study.Opposed to this, for the solid lines the overburden effects are first removed using primary 1 (above the reservoir) as reference, and only thereafter the baseline and monitor study are correlated to get the time-lapse changes inside the reservoir.

Figure 5 .Figure 6 .Figure 7 .
Figure 5. Change in saturation (left column) and P-wave velocity (right column) over time.On the right side the blue and red colours indicate a decrease and increase in velocity, respectively.

Table 1 .
Properties for each layer and fluid for the first numerical example, layers are displayed on the left of Figure2.The asterisk indicates the reservoir layer.In this table the elastic parameters are represented by the Poisson ratio ν = λ /(2λ + 2G) and Young's modulus E = (G(3λ + 2G))/(λ + G).

Table 1
for the first model, and Table2for the second model.

Table 2 .
Properties for each layer and fluid for the second numerical example, layers are displayed on the right of Figure2.The asterisk indicates the reservoir layer.