Monitoring temperature changes in geothermal applications is crucial to ensure sustainable heat production and storage operations. This work focuses on a geothermal project situated on the campus of Delft University of Technology in the Netherlands. In connection to deep low-enthalpy geothermal reservoir exploration, an aquifer thermal energy storage system for the purpose of seasonal shallow heat storage is planned. To enable monitoring changes in the electrical resistivity distribution due to heat injection and extraction operations in the shallow subsurface, a new 480 m deep borehole will be equipped with an electrode setup. Measuring the vertical component of the electric field in the borehole using a frequency-domain surface-to-borehole controlled-source electromagnetic setup is particular effective for monitoring reservoir changes. Due to corrosion effects, conventional electrodes have rather limited lifespans, which may not be sufficient for the multi-decade operational plan for this geothermal application.
We use a new approach integrating capacitive electrodes in composite borehole casings. Tests in shallow boreholes have shown comparable results to standard electrodes. Integrating the capacitive sensors in the composite borehole casing is rather time- and cost-intense requiring pre-drilling installation and specially designed electronics. Therefore, we want to optimise the electrode placement along the borehole trajectory. We simulate vertical electric fields at closely spaced receivers along the borehole trajectory for different subsurface scenarios originating from resistivity changes introduced by injecting and extracting hot fluids. Applying the Ramer–Douglas–Peucker algorithm, we determine which electrode locations and combinations are optimal to record resulting variations in the vertical electric field component. The proposed methodology for optimal electrode placement is promising to improve monitoring efficiency in geothermal applications, ensuring sustainable and effective operation over extended periods.@en

Electromagnetic induction measurements from multi-coil configuration instruments are used to obtain information about the electrical conductivity distribution in the subsurface. The resulting inverse problem might not have a unique and stable solution. In that case, a local inversion method can be trapped in a local minimum and lead to an incorrect solution. In this study, we evaluate the well-posedness of the inverse problem for two and three-layered electrical conductivity models. We show that for a two-layered model, uniqueness is ensured only when both in-phase and quadrature data are available from the measurements. Results from a Gauss–Newton inversion and a lookup table demonstrate that the solution space is convex. Furthermore, we demonstrate that for even a simple three-layered model, the data contained in such measurements are insufficient to reach a correct or stable solution. For models with more than 2 layers, independent prior information is necessary to solve the inverse problem. The insights from the numerical examples are applied in a field case.

@en

The data assimilation process for geothermal reservoirs often relies on well data, which primarily offer insights into the immediate vicinity of the borehole. However, integrating geophysical methods can provide valuable information beyond well proximity, possibly enhancing reservoir predictions. Current methods of monitoring geothermal reservoirs struggle to maintain a good signal-to-noise ratio for deep reservoirs. Diffusive electromagnetic (EM) methods can be sensitive to the decreasing conductivity from heat extraction in geothermal reservoirs and offer promising additional value. To test their potential effectiveness, numerical examples are simulated. A scheme to incorporate diffusive EM observations into a data assimilation process for geothermal reservoirs is presented and implemented in this study. First, an ensemble of prior models representing the reservoir uncertainty is used to determine the moments of the resulting temperature field using a forward geothermal simulation. Subsequently, a conductivity model is calculated from the temperature field using an empirical relationship. The expected electric field response can then be simulated using an EM forward model. EM sources are placed on the surface around the expected cold plume location. The receiver is placed at reservoir depth. To assimilate the data, the ensemble smoother with multiple data assimilation method is used. The findings demonstrate that the incorporation of EM data provides valuable information regarding the temperature field. This improves the accuracy of the temperature forecast of the entire reservoir when combined with the localized data from the production well and, therefore, helps to resolve the complex migration of the cold front. These results highlight the monitoring potential of EM observations for geothermal reservoirs.@en

In agriculture, there is a demand for new methods to monitor the dynamics of fresh rainwater lenses overlaying on saline seeping groundwater. For this purpose, integrating different geoelectrical measurements is a non-invasive and low-cost approach to obtaining subsurface information. Geoelectric methods such as electromagnetic induction (EMI) and electrical resistivity tomography (ERT) have proven effective in characterizing subsoil electrical properties, which can be correlated to petrophysical properties such as fluid salinity. These methods have different sensitivities and can provide complementary information about the electrical conductivity and geometry of the subsurface. This study explores the effectiveness of a methodology that combines EMI measurements with laterally constrained inversion as prior information for ERT inversion. We investigate the usefulness of the method using synthetic data and data from a coastal Dutch polder system. The findings are promising, demonstrating improved delineation of changes in electrical conductivity, potentially linked with salinity fluctuations in the subsoil. This methodology proves effective in mapping in-depth variations in electrical conductivity. It could facilitate the impact assessment of level-controlled drainage systems on augmenting shallow rainwater lenses and mitigating salinization in Dutch polders.@en

Marchenko-type integrals typically relate so-called focusing functions and Green's functions via the reflection response measured on the open surface of a volume of interest. Originating from one dimensional inverse scattering theory, the extension to two and three dimensions set in motion various new developments regarding imaging in complex materials. This extension, however, is based on wavefield decomposition inside the volume and a truncated medium state, i.e. a version of the medium that is reflection-free underneath the focusing location, suggesting that evanescent, refracted and diving waves cannot be included in the representation. We elaborate on a new derivation for Marchenko-like integrals that (i) extends the concept of wavefield focusing by using a generalised homogeneous Green's function, (ii) is based on partial differential equations and thereby allows for additional insights and a new physical intuition for Marchenko equations, (iii) unifies wavefield focusing for open and closed boundary systems, (iv) does not require wavefield decomposition or a truncated medium state, thus including the full wavefield Green's function, (v) enables using forward modelling to obtain, e.g., Marchenko-type, time-compact focusing functions. We place a particular focus on the latter point, illustrating and investigating how to solve the underlying partial differential equations for various types of focusing functions. This paves the way for a deeper understanding of focusing functions as well as advanced full wavefield Marchenko schemes. While the derivations are generally presented for the 3D case, we show numerical examples in 1D.

@en

Green’s functions and propagator matrices are both solutions of the wave equation, but whereas Green’s functions obey a causality condition in time (G = 0 for t < 0), propagator matrices obey a boundary condition in space. Marchenko-type focusing functions focus a wave field in space at zero time. We discuss the mutual relations between Green’s functions, propagator matrices and focusing functions, avoiding up-down decomposition and accounting for propagating and evanescent waves. We conclude with discussing a Marchenko-type Green’s function representation, which forms a basis for extending the Marchenko method to improve the imaging of steeply dipping flanks and to account for refracted waves. @en

Rigid boom electromagnetic surveys that use coil-coil configurations are often used to obtain information about the subsurface conductivity. Semi-analytic solutions help to simulate electromagnetic induction measurements for a large number of horizontally layered models, which can then be stored and used as a lookup table. This procedure is performed once and then used to find the corresponding model that produces the best data fit, eliminating the need for running numerous simulations in every minimization step of an inversion scheme for large field datasets. We apply this methodology to a numerical example and field data acquired in The Netherlands. Our results from both cases using the global search demonstrate its ability to estimate electrical conductivity distributions in two-layered models in a fast and accurate manner. Furthermore, we apply the workflow using a lookup table based on low induction number approximation-derived measurements. The outcome of implementing this methodology using the low induction number lookup table shows poor accuracy in the electrical conductivity estimations for both the numerical example and the field data in comparison to the semi-analytical approach.@en

Modeling ground penetrating radar (GPR) reflection data on glaciers with methods that require the discretization of the full subsurface domain is extremely computationally expensive because of the combination of a large domain size and the comparatively short wavelength of the signal. To address this issue, we build on and extend a previously proposed approach based on the assumption of a homogeneous background medium (ice) in which various scattering objects (e.g., crevasses, channels, boulders) are embedded far from each other such that multiple scattering can be ignored. The glacier bed, below which no scatterers are assumed to exist, represents the lower limit of the modeling domain. With this method, the two-way propagation of the radar waves through the ice is simulated in a semi-analytical way, whereby scattering surfaces are represented with a set of planar elements of different electric and reflective properties, allowing a wide range of objects to be simulated. As we also take the antenna radiation pattern at the air-ice interface into account, this simple algorithm is able to produce realistic 3D GPR data in a fast and memory efficient way. In this study, we validate the presented algorithm with an analytical solution for a layered model, and we simulate radar data for a model of the Otemma glacier in Switzerland featuring a realistic topography of the glacier bed and a subglacial channel.

@en

Natural gas hydrates have been an unconventional source of energy since the beginning of this century. Gas-hydrate-filled reservoirs show higher resistivity values compared with water-filled sediments. Their presence can be detected using marine controlled-source electromagnetic methods. We classify acquisition configurations into stationary and moving receiver configurations, which are described in terms of the design group, the operational details, and where they have been used successfully in the field for natural gas hydrate exploration. All configurations showed good numerical results for the detection of a 700 m long gas hydrate reservoir buried 200 m below the seafloor, but only the stationary configurations provided data that can be used to estimate the horizontal boundaries of the resistive part of the reservoir when the burial depth is known from seismic data. We discuss the operational steps of the configurations and provide the steps on how to choose a suitable configuration. Different CSEM configurations were used together with seismic data to estimate the edge of the gas hydrate reservoir and the total volume of the gas hydrates, to optimize the drilling location, to increase production safety, and to improve geological interpretations. It seems that CSEM has become a reliable method to aid in the decision-making process for gas hydrate reservoir appraisal and development.@en

It is not practical to obtain a large number of labeled data to train a supervised learning network in tunnel lining nondestructive testing with ground-penetrating radar (GPR). To decrease the dependence of supervised learning on the number of labeled data, an improved self-supervised learning algorithm—self-attention dense contrastive learning (SA-DenseCL)—is proposed and incorporated with a mask region-convolution neural network (Mask R-CNN), which is trained by unlabeled and labeled GPR data. The proposed SA-DenseCL adds a self-attention-based relevant projection head to the DenseCL architecture of self-supervised learning, capturing the spatially continuing information between adjacent GPR traces. In the workflow, some unlabeled GPR images are used to pre-train the SA-DenseCL network for feature extraction and obtaining the backbone weights, which is superior to the conventional pre-training methods of supervised learning pre-trained by ImageNet images. The weights of the pre-trained backbone are then used to initialize the Mask R-CNN through transfer learning. Subsequently, a limited number of labeled GPR images are used to fine-tune the Mask R-CNN for automatically identifying the locations of the reinforcement bars and voids and estimating the secondary lining thickness. The experimental results show that the average precision reaches 96.70%, 81.04%, and 94.67% in identifying reinforcement bar locations, detecting void defects, and estimating secondary lining thickness, respectively, which outperform the conventional methods that use ImageNet-based supervised learning or GPR image-based DenseCL for initializing the Mask R-CNN backbone weights. It is observed that the improved self-supervised learning-based framework can improve the detection and estimation accuracy in GPR tunnel lining inspection.

@en

In hydrocarbon drilling, mud filtrate penetrates permeable formations and alters the pore fluid characteristics in the immediate vicinity of the borehole. Typically, the prevailing in situ pore fluids are displaced by the invading mud filtrate, which leads to gradually changing distributions of the fluid and electrical properties. Understanding this invasion process is crucial for the interpretation of logging data and associated reservoir evaluations. Conventional logging methods tend to be inadequate for this purpose as their resolution is too low. We find that invasion depth can be determined from borehole radar data using an optimized antenna configuration and time-lapse measurements. A series of parametric sensitivity analyses provide information about the effects of variations of the rock and fluid properties on the identification and extraction of borehole radar signals reflected from the invasion front. Our results suggest that by embedding the radar antennas in cavities filled with an absorbing dielectric material, it is possible to minimize the interference arising from the metal components of the logging tool. In the simulated reservoir scenario, a time-lapse measurement mode with a time interval of at least 6 h can reliably extract the radar signals reflected from the invasion front, and the proposed borehole radar has a lateral detection range from 0.15 to 1 m. A comprehensive range of parametric sensitivity analyses indicates that the signals reflected from the invasion front are principally influenced by oil viscosity, porosity, and mud and formation water salinity, as well as by molecular diffusion coefficient and cementation exponent. These properties and parameters should be carefully explored and assessed when applying borehole radar to evaluate mud invasion information in a reservoir environment.

@en

Tracking temperature changes by measuring the resulting resistivity changes inside low-enthalpy reservoirs is crucial to avoid early thermal breakthroughs and maintain sustainable energy production. The controlled-source electromagnetic method (CSEM) allows for the estimation of sub-surface resistivity. However, it has not yet been proven that the CSEM can monitor the subtle resistivity changes typical of low-enthalpy reservoirs. In this paper, we present a feasibility study considering the CSEM monitoring of 4–8 Ω·m resistivity changes in a deep low-enthalpy reservoir model, as part of the Delft University of Technology (TU Delft) campus geothermal project. We consider the use of a surface-to-borehole CSEM for the detection of resistivity changes in a simplified model of the TU Delft campus reservoir. We investigate the sensitivity of CSEM data to disk-shaped resistivity changes with a radius of 300, 600, 900, or 1200 m at return temperatures equal to 25, 30, …, 50 °C. We test the robustness of CSEM monitoring against various undesired effects, such as random noise, survey repeatability errors, and steel-cased wells. The modelled differences in the electric field suggest that they are sufficient for the successful CSEM detection of resistivity changes in the low-enthalpy reservoir. The difference in monitoring data increases when increasing the resistivity change radius from 300 to 1200 m or from 4 to 8 Ω·m. Furthermore, all considered changes lead to differences that would be detectable in CSEM data impacted by undesired effects. The obtained results indicate that the CSEM could be a promising geophysical tool for the monitoring of small resistivity changes in low-enthalpy reservoirs, which would be beneficial for geothermal energy production.@en

This study presents a novel methodology to transform 1D resistivity data into layered resistivity models without prior information by using the concept of cumulative reference models. The proposed methodology involves deriving an error function that transforms apparent resistance measurements into a cumulative resistance, which is then transformed into a layered resistivity model. We applied the methodology to simulated data from various 1D models with different physical parameters, and the results demonstrate that our method can be used to directly transform the data into a layered resistivity model without requiring prior information. This methodology provides a valuable alternative to inversion methods when one local model is available and multiple measurements are available over an area with similar physical parameters. Furthermore, the retrieved rescaled model can be used as a reference model for the inversion process, reducing computational and economic costs. This study highlights the potential of cumulative reference models for subsurface characterization, providing a new paradigm to study the subsurface with increased efficiency.@en

Inspired by recent developments in wave propagation and scattering experiments with parity-time (PT) symmetric materials, we discuss reciprocity and representation theorems for 3D inhomogeneous PT-symmetric materials and indicate some applications. We start with a unified matrix-vector wave equation which accounts for acoustic, quantum-mechanical, electromagnetic, elastodynamic, poroelastodynamic, piezoelectric and seismoelectric waves. Based on the symmetry properties of the operator matrix in this equation, we derive unified reciprocity theorems for wave fields in 3D arbitrary inhomogeneous media and 3D inhomogeneous media with PT-symmetry. These theorems form the basis for deriving unified wave field representations and relations between reflection and transmission responses in such media. Among the potential applications are interferometric Green’s matrix retrieval and Marchenko-type Green’s matrix retrieval in PT-symmetric materials. @en

By solving a Marchenko equation, Green’s functions at an arbitrary (inner) depth level inside an unknown elastic layered medium can be retrieved from single-sided reflection data, which are collected at the top of the medium. To date, it has only been possible to obtain an exact solution if the medium obeyed stringent monotonicity conditions and if all forward-scattered (non-converted and converted) transmissions between the acquisition level and the inner depth level were known a priori. We introduce an alternative Marchenko equation by revising the window operators that are applied in its derivation. We also introduce an auxiliary equation for transmission data, which are collected at the bottom of the medium, and a coupled equation, which is based on both reflection and transmission data. We show that the joint system of the Marchenko equation, the auxiliary equation and the coupled equation can be succesfully inverted when broadband reflection and transmission data are available. This results in a novel methodology for elastodynamic Green’s function retrieval from two-sided data. Apart from these data, our approach requires P- and S-wave transmission times between the inner depth level and the top of the medium, as well as two angle-dependent amplitude scaling factors, which can be estimated from the data by enforcing energy conservation. @en

Green’s functions in an unknown elastic layered medium can be retrieved from single-sided reflection data by solving a Marchenko equation. This methodology requires a priori knowledge of all forward-scattered (non-converted and converted) waveforms. Moreover, the medium should satisfy stringent monotonicity conditions, which are often not met in realistic scenarios. In this contribution, we show that the situation is significantly less cumbersome if two-sided reflection and transmission data are recorded (for instance in laboratory settings). A novel methodology is presented to retrieve elastodynamic Green’s functions from such data. Apart from the two-sided reflection and transmission responses, our methodology requires knowledge of the direct non-converted PP- and SS-transmissions (a priori knowledge of forward-scattered converted waveforms is not needed). We demonstrate the success of our methodology by conducting a numerical experiment in an elastic layered medium that violates the monotonicity conditions of the Marchenko equation for single-sided reflection data. The limitations of the methodology and the sensitivity to errors in our initial estimates require further investigation.@en

Surface waves are widely used to model shear-wave velocity of the subsurface. Surface wave tomography (SWT) has recently gained popularity for near-surface studies. Some researchers have used straight-ray SWT in which it is assumed that surface waves propagate along the straight line between receiver pairs. Alternatively, curved-ray SWT can be employed by computing the paths between the receiver pairs using a ray-tracing algorithm. The SWT is a well-established method in seismology and has been employed in numerous seismological studies. However, it is important to make a comparison between these two SWT approaches for near-surface applications since the amount of information and the level of complexity in near-surface applications are different from seismological studies. We apply straight-ray and curved-ray SWT to four near-surface examples and compare the results in terms of the quality of the final model and the computational cost. In three examples we optimise the shot positions to obtain an acquisition layout which can produce high coverage of dispersion curves. In the other example, the data have been acquired using a typical seismic exploration 3D acquisition scheme. We show that if the source positions are optimised, the straight-ray can produce S-wave velocity models similar to the curved-ray SWT but with lower computational cost than the curved-ray approach. Otherwise, the improvement of inversion results from curved-ray SWT can be significant.

@en

Delft geothermal project (DAPwell) is a planned geothermal well doublet, where relatively cold water is going to be injected through one well into a low enthalpy geothermal reservoir to produce hot water from the other well. The volume of the cold water around the injection well will increase over time and, in the end, result in a thermal breakthrough. Thus, it is essential to trace the time-lapse change in the volume of the cold water to monitor the DAPwell efficiently. The invaded reservoir volume by the cold water is associated with a decrease in the pore fluid temperature and salinity. This increases the electrical resistivity of the geothermal reservoir, where the cold front is located. Hence, estimating the time-lapse change in the electrical resistivity of the geothermal reservoir can be used to identify the distribution of the cold water. From a theoretical point of view, the controlled-source electromagnetic (CSEM) method can be used to get information about the change in the electrical resistivity within the geothermal reservoir. In this study, we investigate the feasibility of monitoring a geoelectric model of the DAPwell using land CSEM forward modelling. The optimal survey design is investigated as well as the influence of cold water volumetric changes on the time-lapse electric field response. The impact of measurements undesired effects on time-lapse CSEM response is analysed and then synthesized. A subsurface model of the DAPwell is illuminated by a horizontal electric dipole source, which emits a sinusoidal field with several frequencies. Based on the numerical experiments, surface measurements do not pick up sufficient time-lapse signal to use them for field applications. On the other hand, the difference in the z-component of the electric field, determined along a depth section, allows for a successful detection of the electrical resistivity changes within the geothermal reservoir. The correlation between the spatial distribution of the cold water and the difference in time-lapse electric field responses is clarified. Finally, it is noticed that the difference in time-lapse signal is measurable in the presence of the different sources of noise. @en

Joint inversion of different geophysical methods is a powerful tool to overcome the limitations of individual inversions. Body wave tomography is used to obtain P-wave velocity models by inversion of P-wave travel times. Surface wave tomography is used to obtain S-wave velocity models through inversion of the dispersion curves data. Both methods have inherent limitations. We focus on the joint body and surface waves tomography inversion to reduce the limitations of each individual inversion. In our joint inversion scheme, the Poisson ratio was used as the link between P-wave and S-wave velocities, and the same geometry was imposed on the final velocity models. The joint inversion algorithm was applied to a 2D synthetic dataset and then to two 2D field datasets. We compare the obtained velocity models from individual inversions and the joint inversion. We show that the proposed joint inversion method not only produces superior velocity models but also generates physically more meaningful and accurate Poisson ratio models.

@en

Natural gas hydrates production tests over the last two decades has sown that production is not without risks. Indirect effects in the sedimentary rocks of phase changes are changes in porosity, permeability, and saturation. From a field production test site, porosity changes in the range of 15% to 19% and saturation from 5% to 60% were reported. Monitoring is in principle possible using an electromagnetic survey with a downhole vertical electric source and a horizontal electric field receiver on the seafloor. Computed model responses over a wide frequency range and for many depth locations of an electric current source show that both changes can be detected. Best detectability occurs when the current source is below the reservoir layer in case of changes differences can be detected above, inside and below the reservoir layer at frequencyies below 10 Hz. At a source operating frequency of 0.1 Hz maximum response difference between the two values in saturation occur when the source is 20 m above the top of the reservoir layer unil 100 m below the bottom. Only below the top of the reservoir there is almost no difference in the electric field amplitude between the two saturation levels below 10 Hz.

@en