Crosswell EM uses the principles of EM induction to interrogate the electrical resistivity between wells. As sketched in Figure 1B, a vertical magnetic dipole (VMD) induction coil placed in the transmitter well broadcasts sinusoidal EM signals into the formation. At the receiver well, the signals are detected, in our case via a double-coil receiver with two connected sensors separated by 5 m. Source and receiver usually traverse the wells with recordings triggered at regularly spaced intervals leading to a tomographic illumination of the interwell space.

Tomographic baseline crosswell EM data were collected over the depth section 200–320 m (Figure 1C) with a transmitter spacing of 2 m in observation borehole 2 (named OB2), thus leading to 59 source depth levels. The two-receiver string traversed 10 levels at intervals of 10 m in observation borehole 1 (OB1). Unfortunately, cross checks between signal amplitudes and receiver depth recordings revealed inconsistencies which rendered only 8 receiver sets per transmitter activation usable. While the cause remains unclear, we speculate that the bottom receiver of the double-sensor array might have buoyantly floated upwards, thus not keeping the constant reference separation (5 m) between the two sensors. Moreover, while the transmitter system operated at the two frequencies of 200 and 450 Hz, data noise turned out excessive at 450 Hz. Both issues reduced the usable data to a total of 454 complex magnetic-field data points which are the input to our field-data inversions presented below.

In addition to the data sparsity, another complication is caused by horizontal well deviations of the instrumented sections which max out at 7 m (OB2) and 0.5 m (OB1). Ensuing instrument rotations out of the horizontal plane are often neglected in academic crosswell EM studies. Given the dominating deviation of OB2, we incorporate source coil rotations in form of tilt and azimuth angles in the field-data inversions.

Croswell EM can be stymied by the signal-shielding nature of metal borehole casings, because large magnetic fields result from the transmitter coupling with hysteretic magnetic permeabilities of steel. If open holes (no casing present) cannot be realized, fiberglass casing avoids this coupling. At CaMI.FRS, the transmitter well (OB2) is steel-cased from the surface to 60 m and fiberglass-cased from 60 to 350 m, while OB1 is all steel-cased down to 350 m. The effects at frequencies below a few hundred Hz remain moderate for the received signals. Nevertheless, in the presence of weak geological effects on the data one wants to avoid any potentially overpowering casing effect. The casing effect at the receivers can be corrected for either in the course of data inversion, or it can be removed prior to the inversion process. Here we choose the latter in order to avoid the computationally demanding forward-modeling of casing effects during the inversion. The correction procedure essentially utilizes the fact that in the frequency domain one can approximate the casing effect by means of a source-signature correction with a complex factor (Wu and Habashy, 1994).

We outline the steps that correct for EM induction effects contained in receiver coils that are surrounded by steel casing. The main component is an appropriate reference model from which one can calculate non-casing reference data. The two resistivity logs in the two observation wells (Figure 2A) were employed to construct a 1D reference model. It involved fewer parameters than given by our 1D inversion presented below while being equivalent in terms of data fit. The predominantly horizontally layered geology at CaMI.FRS helps limit the error due to the approximate nature of the reference model.

For each common-receiver gather of the tomographic crosswell EM data, we extract the data pertaining to all transmitter depths that are in vertical proximity, here within 5 m of the depth of the receiver gather. These data are complex magnetic-field values containing the casing effect. Replicating their survey geometry with respect to the receiver depth, we calculate a corresponding non-casing data set through forward-modeling of the 1D reference model.

For each transmitter depth of a given receiver gather, we divide the non-casing reference data by its measured counterpart. This yields a set of estimated casing-correction factors which are also complex numbers. All correction factors of one receiver gather are averaged. The outcome is a best estimate for the magnetic-field attenuation and phase shift caused by the steel casing. The averaging is carried out for all individual receiver gathers, producing a set of location-dependent best-estimate correction factors. Making the correction depth-dependent intends to limit the impact on non-casing signal footprints one wants to retain as they contain for example effects due to local 3D structures.

The final step is the actual data correction. Each data point of a given receiver gather is multiplied by the gather’s correction factor, producing a common-receiver gather of corrected data. Applying the correction procedure to each receiver gather yields a tomographic data set of casing-corrected data that are also normalized for unit source moment. This set is the data input for the field-data inversions presented here. Image distortions due to this correction procedure are limited to the vicinity around the casing (Gao et al., 2008).

Compressional-wave velocity ( v p ) logs were acquired in OB1 and OB2 (Figure 2A). Cross-plotting the well measurements of velocity against resistivity (Figure 2B) highlights a degree of correlation that we deem as sufficient for using v p as reference for our EM inversions. In addition to the well logs, a high-resolution crosswell seismic survey with piezoelectric seismic sources was part of the baseline data collection at CaMI.FRS (Um et al., 2020). To achieve tomographic imaging, a total of 257 total source positions traversed the well from 216 to 344 m at depth intervals of 0.5 m. Figure 3A sketches the set of ray paths resulting from one source activation. Seismic signals were recorded in well OB2 using a string of broadband hydrophone receivers. The receiver interval of 5 m resulted in a vertical coverage of 95 m. The whole survey resulted in 51,400 ray paths. Corresponding travel times were picked for tomographic travel-time inversion (Vasco and Nihei, 2019). The travel time inversion was initialized by a 1D velocity model obtained through linear interpolation of the v p -logs from OB1 and OB2.

Both the v p -logs as well as the final 2D v p -image (Figure 3B) are the prior information that will contribute to the design of the parameter bounds for the field-data inversions presented below. First, we present a synthetic-data proof-of-concept in order to highlight the potential benefits of the method.

Our underlying premise is to employ bound constraints as a way of focusing the EM imaging by exploiting prior information. Our synthetic example simulates a case where prior information exists in form of a delineation of a focus zone, which could be the product of concurrent reservoir modeling. The bound constraints thus serve as a way of augmenting the geophysical model resolution in a zone that is of hydrogeological interest.

The zone of interest in our electrical-conductivity model will be a plume of ellipsoidal shape and with anomalous electrical conductivity relative to the background. We consider two cases, a resistive and a conductive plume. The ellipsoid’s semiaxes lengths (in meters) are (a, b, c) = (18.5, 10.5, 12.5) along the axes x, y, z, and it is centered at (x, y, z) = (22, 0, 257). White contour lines in Figure 4A outline the true target geometry. The colored ellipsoidal regions that occur incongruent to this actual geometry define the volume where lower and upper inversion parameter bounds are modified with respect to constant default bounds.

A prerequisite for designing modified bound constraints would be a reservoir flow model that approximates fluid saturations and thus plume evolution. While we do not carry out such flow-modeling, it is safe to assume that, at a minimum, a basic reservoir model will exist for the majority of GCS sites. However, depending on the degree of model refinement, there will always be some uncertainty associated with operational reservoir flow predictions. The incongruity between true and assumed plume boundaries is meant to simulate this uncertainty. To quantify the deviations, we assume that predictive flow modeling would deliver a plume estimate with the erroneous ellipsoid dimensions of (a, b, c) = (20.5, 12.5, 10.5). In addition, the center point is shifted to (x, y, z) = (18, 0, 260), which shall account for actual local preferential flow paths that would exceed the complexity of the available flow model.

Lower (Figure 4A) and upper (Figure 4B) σ -parameter bounds are designed according to Eqs 5–8. As illustrated in Figure 4C, preset default bounds [ a d e f , b d e f ] = [ 10 − 3 , 1.5 ] expand within the plume region, up to the maximal interval [ a m i n , b m a x ] = [ 10 − 6 , 100 ] at the center. Figure 4D exemplifies the parameter transformation (Eqs 3, 4) that would result for this widest interval. Generally, the asymptotic curve behavior indicates the logarithmic transformation’s effect of avoiding σ -parameter over- and undershoots during the inversion process. Widening the bound interval corresponds to increasing the parameter range in which the relationship between original and transformed parameter is close to linear.

First, we assume the simplistic case where a free CO₂ phase would act as an electrical insulator, thus causing a resistive plume. Figure 5A compares inversions using spatially constant bounds versus variable bounds. All inversions start with a homogeneous model of σ = 0.07  S/m ( ≈ 14.3   Ω m ). The first plot row shows the actual resistive anomaly to be imaged. It would be caused by an injection source at (x, y, z) = (22, 0, 257), creating the aforementioned ellipsoidal plume geometry (Figure 4A, white contour lines). Mimicking highest CO₂-saturation levels near the injector, σ gradually increases from a minimum of σ = 2.34 × 10 − 3  S/m ( ≈ 427   Ω m ) at the center towards the borders where it blends with the background. The background has a generally 1D structure with minor 3D anomalies and was obtained through transforming an image of seismic velocity v p (Figure 3B) into ρ through the relation lo g 10 ( ρ ) = m ⋅ lo g 10 ( v p ) + c , with the constants m = 3.88 and c = − 11 . This empirical relation is chosen here arbitrarily as it was found to explain ρ - v p correlations for near-surface materials (Meju et al., 2003).

The 3D model outcome in Figure 5A (second plot row) results from constant parameter bounds [ a d e f , b d e f ] = [ 10 − 3 , 1.5 ]  S/m, to be compared against the counterpart with variable bounds (third plot row). Most noticeable is the weak target anomaly reproduction in both images, which is in alignment with data sensitivity limitations. It is a general issue that VMD crosswell systems have limited capacity in resolving resistive targets that are relatively thin with respect to instrumented well sections (e.g., Grayver et al., 2014). At the same time, the inversion reproduces some of the conductive background structures, which is mostly visible in the x-z-sections of both images.

At low porewater salinities and injection levels that are shallow enough so that CO₂ remains in the gaseous state, dissociation of carbonic acid causes higher ion contents. The resulting electrically conductive plume is simulated here at the injection depth of 257 m. The corresponding model in Figure 5B (true model in first plot row) uses the same background conductivity and plume geometry as the resistive plume. Now, the conductivity decreases from a maximum of 1.7 S/m at the ellipsoid center towards the background values at the borders.

Compared to a resistive anomaly, the sensitivity of VMD systems increases for a conductive anomaly. Similar to its resistive counterpart, the inversion with constant bounds produces a plume image that remains faint (Figure 5B, second row). However, the usage of the variable bounds achieves a clearer plume delineation (third row).

We emphasize that the variable-bounds concept has the mere goal of improving model resolution in terms of boosting the NLCG search process. It cannot alleviate a deficient data sensitivity. Insufficient sensitivity would cause solution instability. Preceding sensitivity studies therefore ensure that model parameters of target zones are well-resolved, regardless of any type of inversion constraints. A basic sensitivity check compares H _z -field responses calculated for the 59 receivers due to a selected source position located at z = 259 m (Figure 6). Amplitude (upper plot) and phase (lower plot) data are calculated for the background, where no ellipsoidal anomaly is present, and plotted together with the data containing the signal due to the resistive (blue lines) and conductive (red lines) plume. Percentage differences (dashed lines) between background and anomalous responses confirm that the sensitivity is higher for the conductive target (red dashed lines). While the degree of sensitivity is generally low, it is well above the assumed noise level (here 0.5%), thus validating that the sensitivity is sufficient for this kind of anomaly. Moreover, percentage differences (dashed lines) that are above noise level appear over a slightly larger depth section for phase than for amplitude. This sensitivity behavior made us set the data weights pertaining to the imaginary field component twice as high as for the real component. Since this is an arbitrary choice, it is likely that the plume imaging might be further optimized by means of additional trial inversions with proper differential data weighting.

Figure 7A summarizes all four inversion instances (resistor and conductor, with constant and variable bounds) by showing the total RMS data fit against the model-updating iterations, the latter set to a maximum of 200. All inversions reach an RMS error below one. When facing generally low data sensitivities and such a similar general error convergence, the quality differences between different solutions are usually hardly discernible from individual data fits. Figure 7B exemplifies both initial and final fits produced by the variable-bound inversion for the conductive target. Here, maximal differences between the data fits at iterations 0 and 200 amount to 27% for amplitude and 31% for phase. Data fits for the corresponding constant-bound result (not shown here) are visually identical. These observations give rise to a verification of the solution consistency. Specifically, we want to ensure through two tests that 1) the anomalies are not an imaging artifact caused by the bound constraints, and 2) the solution is sufficiently independent of the starting model.

The first test for solution consistency is in form of an inversion for the background model. For this check, the inverted data is calculated from the model without any anomalous ellipsoid present, that is, the data contains only the footprint of the background. At the same time, the inversion uses the variable bounds. Figure 5C shows the result of this test, where the x-y section exhibits a minor artifact reflecting the zone of the augmented bounds. However, generally, the final image indicates neither a resistive nor a conductive anomaly.

The second test verifies solution consistency over a range of starting model values. Figure 8A summarizes a sequence of 24 variable-bound inversion realizations for the conductive plume target. The starting σ spans the range from (logarithmically) 0.013 S/m to 0.32 S/m. A simple RMS measure tracks the difference between each final model and the true model (true model shown by Figure 5B, upper row). A first observation is that the course of the RMS error does not appear erratic with changing starting σ , which would be the case for a more ill-posed inverse problem. Instead, the RMS becomes minimal within the interval [ σ 1 , σ 2 ] . Figure 8B shows the images for the two starting values σ 1 and σ 2 . Both images contain features of the conductive target, indicating independence from the starting σ over a range that exceeds an order of magnitude.

We outline a four-step procedure for the design of target-adapted lower and upper parameter bounds for a constrained 3D inversion. The main principle is the widening of a given default bound interval [ a d e f , b d e f ] → [ a d e f − Δ a , b d e f + Δ b ] based on the abovementioned inversion objective. The degree of widening, quantified by Δ a and Δ b , is non-zero in zones that are deemed anomalously conductive with respect to a given reference and increases with the degree of discrepancy. We reiterate that reliable prior information is important for this procedure. In our case, this information is provided by the v p and ρ well logs (Figure 2A) and a 2D image of seismic velocity (Figure 3B). Figure 10 illustrates the complete procedure of bounds construction with its steps described in the following.

Two comparative field data inversions employ the constant bounds [ a m i n , b m a x ] = [ 10 − 7 , 15 ] and the variable bounds created through the mapping scheme described above (Figure 10). The initial reasoning for choosing a very wide constant-bound interval was to let the inversion practically be unbounded in order to avoid any bias due to potentially too tight bounds. However, other constant-bound inversions not reported here experimented with narrower intervals which yielded no significant difference in the final image.

The initial model guess for both inversions is given by the final mean 1D model of Figure 9 (referred to as 1D model in the following). One could assume that the v p - σ lookup table that comes out of the GPR modeling provides a better starting σ -model (Figure 10, Step 2) as it combines the layered geology with local inhomogeneities. However, we employ this model merely for the identification of conductive anomalies and the corresponding bound modifications. Its values span a range that exceeds the σ -range of the resistivity logs and 1D model by an order of magnitude, which would result in a much greater initial RMS misfit.

The low final RMS misfit error of 1.78 achieved by the 1D model equals the starting RMS for the 3D inversions, confirming the dominating layered structure of the geology suggested by Um et al. (2020). Figure 11A summarizes the total error convergence. The variable-bound 3D inversion converges to a slightly lower final RMS and reaches the preset maximum iteration number (250), while the constant-bound inversion terminates after 173 iterations. However, final data fits are visually practically identical for both inversions, as evident in fits exemplified for one selected source depth (Figure 11B).

It is now of interest to what degree variable bounds contribute to a better target inversion. While synthetic-data trials make an assessment straightforward, the absence of a ground truth in field-data inversions makes us revert to examining model differences. Figure 12 shows relative σ -differences with respect to a reference. First, we want to visualize the difference Δ σ v a r = M v a r 3 D − M 1 D between the final variable-bound inversion result ( M v a r 3 D ) and the starting 1D model ( M 1 D ). Deviations from a horizontally layered structure that are within the center plane of the interwell region are most pronounced within a zone between z = 240 m and z = 260 m (Figure 12A). Moreover, the image Δ σ v a r strongly resembles Δ σ c o n = M c o n 3 D − M 1 D , where M c o n 3 D is the final constant-bound inversion result. Instead of Δ σ c o n , Figure 12B shows the difference Δ σ = M v a r 3 D − M c o n 3 D . Absolute differences are relatively small, indicating no significant bias due to the variable-bound constraints. Further, positive differences dominate within the well plane, which translates to σ ( M v a r 3 D ) > σ ( M c o n 3 D ) . These observations point to some degree of significance for anomalies characterized by Δ σ > 0 that deviate from a horizontally layered geology within the zone z = 240–260 m.