Even though negative effects on the use of crude oil have surfaced over the past years, our energy matrix still largely relies on this energy source. The production of oil, therefore, plays an important role in our society. Unfortunately, the process of oil production is highly uncertain. There are uncertainties associated on the production strategy, e.g. where and how many wells should be drilled and how these wells should operate because of the uncertainties associated with the limited knowledge about the subsurface. In this thesis we are dealing with the uncertainty of the rock permeability distribution. Typically, rock permeabilities in the rock vary, but from the outside this can’t be perceived. If these rock permeabilities are estimated inaccurately, they will result in inaccurate pressure solutions. Then, this can lead to faulty decisions regarding the oil exploitation. To resolve this issue, a data assimilation technique may be applied to correct these model parameters based on mismatch of simulated data and observations. For this optimization technique, often gradient information is required. Since
in reservoir simulation the number of parameters generally is extremely high, computation of this information is computationally expensive. Therefore, a multiscale framework is employed to improve the
computational efficiency of the forward simulation. Multiscale methods are able to solve the model equations at a computationally efficient coarse scale and can easily interpolate this solution to the fine
scale resolution. Next, we use a Lagrangian set-up together with a multiscale framework to re-derive an efficient formulation for the derivative computation. However, as the multiscale method is prone to errors, this derivative computation formulation is recast in an iterative fashion, using a residual based
iterative multiscale method to provide control of these errors. In this thesis we show that this method generates accurate gradients. In contract to the high accuracy of the method, this method comprises a computationally heavy smoothing step. This issue can be resolved by making smart use of the Lagrange
multipliers, to re-derive an efficient iterative multiscale solution strategy. The multipliers are used to identify important domains of the region for which smoothing is required and for which regions we may neglect the smoothing. We show that the newly proposed iterative multiscale goal oriented method is computationally more efficient and we show that method is promising for efficient derivative computation, but that more work is required to fully demonstrate the benefit of this method.

Salt caverns formed by solution mining may cause soil subsidence, because the surrounding adapts to fill the void (cavity) created in the strata. The cavern volume is hence not only a function of salt extraction (estimated from mass balance) but also a function of salt creep. The cavern size can be measured by sonar, but a less costly way would be to correlate cavern volume to cavern compressibility. The size of the cavern can be obtained by a compressibility test. The compressibility test consists of injecting brine or fresh water into the cavern, causing a pressure build-up. The pressure increase depends on the cavern volume and compressibility. Bérest et al. (2006) [1] derived an equation that expresses the pressure response to injection of a volume V of fresh water. The equation includes creep that can be derived from a separate geometric creep model that relates the strain to the stress history using the Norton-Hoff law for describing an Arrhenius type of response model, and thermal expansion effects influenced by the geothermal temperature. It also considers Darcy flow and leakage through the permeable salt cavity. This thesis investigates the sensitivity related to the various mechanisms in the Bérest and Van Sambeek (2005) model. It also extends the model by introducing chemical dissolution of salt effects and taking into account the impact of the sump (insoluble precipitates at the bottom of the cavern). We state the model equations and solve them for an example of interest. We can also compare the calculated pressure response to the results of a pressure test. In this pressure test a given volume of fresh water is injected above the sump and produced at the top of the cavity. The produced salt concentration is also one of the parameters measured in the test. There are at least two widely used methods to obtain or predict cavern size of active salt solution mining caverns, namely sonar surveys and leaching simulators. However, these data can be flawed due to uncertainties of the methods and the accuracy of measuring devices like flow meters. The compressibility test can then be used as an additional survey method. The test execution is simple and inexpensive, making it possible to be performed more frequently than a sonar survey. The cavern volume is obtained from a compressibility test data by using a solution that includes not only the injection volume, but also volumes introduced by other phenomena. These other phenomena are caused by thermal, hydraulic, mechanical, and chemical influences. Previous work has attempted to introduced solutions incorporating these influences for idle caverns. The proposed solution incorporates these effects in an active leaching salt cavern. It also introduces the impact of sump (insoluble on cavern bottom) to be included in the model. The proposed solution is then tested against field data. The work here differs from previous work that it presents a comprehensive description of the various methods used in the literature and practice. Where current models predict cavity volumes with a mean absolute percentage error of 33% (Thiel’s method for BAS3O) to 127% (Bérest et al. for BAS4) for the studied caverns (BAS3O and BAS4), the model here, which includes the presence of insoluble particles, and a well-established creep model is able to predict cavity volumes to an accuracy of 33% (proposed model for BAS3O) to 12% (proposed model for BAS4).

The subsurface of Groningen gas field composed of several faults. The continuous production has been resulting in several micro-seismicity activities, particularly for the past decades. One of the reasons for the production-induced-seismicity is fault reactivation at depth. The fault that initially in the non-active state becomes active and starts slipping due to pore fluid extraction. The modelling of fault reactivation induced by production can be simulated by coupling of flow and geomechanics.

The project studies the coupling of flow and geomechanics for an idealistic subsurface model with fault. The flow occurs due to pore pressure depletion and is predicted by solving mass conservation. The corresponding deformation is estimated by momentum balance equation.

The observation will be limited to rock deformation that is quantified by effective stress and displacement on the fault due to a pressure change in the subsurface. The simulation is performed using Stanford's Automatic Differentiation General Purpose Reservoir Simulator (ADGPRS).

In the proposed model, the flow equation is discretised using finite volume method, and poromechanics is discretised by finite-element Galerkin's approach. Both problems share the same grid model, and there is no error associated with the information exchange. The coupled problem is being solved fully implicitly where flow and mechanics formulation are solved simultaneously.

The case study for the simulation is a simplified Groningen subsurface model with a fault zone. The impact of variability in production dynamic on fault reactivation is studied. The variations in production rate and production strategy are applied in different formation configurations. It is found that the correlation between the production rates and the stress state in the subsurface depends on the formation offset. However, the stress field does not depend on continuity of the production scheme. Several aspects including the presence of Gas-Water contact in the model was studied in this project as well.

Model based optimization of reservoir water flooding is an ill-posed problem where significantly different control strategies deliver near identical Net Present Values (NPVs). Discovering and exploiting the existence of "redundant" control strategies - particularly those in close proximity to an optimal strategy - is valuable since this offers operational flexibility in reservoir management. To identify such 'flexible strategies' this thesis proposes a workflow to characterize the space of feasible solutions. The feasible region or feasible solution space consists of all control strategies that deliver an NPV within some threshold from an optimal value. Ensemble-based optimization is performed with strong Wolfe line search to identify an optimal control strategy. The BFGS scheme is used to iteratively approximate the Hessian matrix. One dimensional exploration is performed along the singular vector directions that characterize the null space of the Hessian. A thorough exploration results in an accurate characterization of the feasible region. Such an approach however is computationally intractable in case of realistic reservoirs with multiple hundred controls. To address this, a high-dimensional polytope is first defined using the end points from the exploration step. Subsequently, an innovative cross-section constrained maximum volume ellipsoid is inscribed within this polytope to generate an ellipsoidal approximation of the feasible region. Validation results are then presented which show that even one ellipsoid centred at the optimum control vector coordinate provides a conservative description of the feasible solution space.

The oil industry is a high risk high reward venture. The capital and operating expenses runs into tens of millions of dollars with oil production being the primary source of revenue. During the initial development phase of the field, few wells are completed. Data collection during this period verifies assumptions made during the modelling phase and forms the basis for future development. History matching can play a significant role in these plans since it can determine uncertainties in future production.

History matching algorithms must be able to make an accurate estimates of uncertainty in future production while being computationally light. Towards this end, a number of history matching methods have been developed with emphasis being on the ensemble Kalman filter (EnKF) in recent times. While the variants of the EnKF seek to improve different aspects of the method, few have been successful in addressing all of these concerns.

The Distributed Gauss Newton (DGN) was developed with the same goal- accurate uncertainty prediction at low computational cost. It is not a variant of the EnKF but uses a sensitivity matrix determined through linear regression which decreases the computational load compared to existing gradient based techniques. In their tests, the authors report superior performance of the DGN compared to a Gauss-Newton scheme. This thesis aims to provide a detailed understanding of the method and its dependencies. This is followed up with a comparison of the DGN with an EnKF variant known as the ES-MDA.