Open Access

ISSN: 2381-8719

+44 20 3868 9735

Research Article - (2014) Volume 3, Issue 5

*
*

Most renewable energy sources are intermittent and need buffer storage (e.g., compressed air energy storage, CAES) to bridge the time-gap between power supply and demand peaks. Replacing pore brine with CAES causes changes in electro-elastic properties and density, and justifies applications of multi-geophysical approach. In this numerical study we apply techniques of the elastic full waveform inversion (FWI), electric resistivity tomography (ERT), transient electromagnetic induction (TEM) and gravity to detect and monitor CAES in deep reservoirs and possible leakages in shallow groundwater aquifers of North Germany. For different subsurface model scenarios of CAES reservoirs and leakages, synthetic data sets are generated and inverted using constraints on the initial model. Results reveal principally the capability of our applied approach to resolve the CAES plume in deep saline reservoirs and shallow groundwater aquifers. The ERT resolution for leakages is highly enhanced for the combined surface-borehole survey compared to the individual surface and borehole surveys. The applied gravity technique is highly sensitive to the mass deficit caused by CAES plume. The detect ability limit of the technique is determined by the least CAES volume causing an anomaly with amplitude just above the accuracy range of modern micro-gravimeters. The FWI technique can map the shallow CAES leakage by anomalies in the reconstructed ΔVp, ΔVs and Δdb tomograms within the background aquifer. However, these tomograms contain inversion artifacts and smearing effects related mainly to the dominance of the Rayleigh wave in the data. Obviously, applied multi-techniques complement and confirm each other. CAES plumes cause strong mass deficits and moderate resistivity highs and thus are more sensitive for gravity and FWI methods. Applying constrained inversion minimizes interpretation ambiguities and helps recovering almost realistic electro-elastic parameters that can be applied in adequate petrophysical equations to quantify CAES saturations.

**Keywords:**
Renewable energy, Compressed air/gas energy storage (CAES), Elastic full waveform inversion, (FWI), Electric resistivity tomography (ERT), Gravity method, Transient electromagnetic induction (TEM), Petrophysical relations and parameters

1D, 2D, 3D: 1-, 2-, 3-dimension of spatial mapping with respect to the coordinates x, y, and z; a: Constant of Archie equation; α, β, γ: Tripotential electrode configurations; B: Magnetic field induction; C: Current electrode; CAES: Compressed air/gas energy storage; d: Bulk density; dbr: Brine density; dg: Gas density; Δgx, Δgy, Δgz Residual component of gravity field with respect to the background; EMI: Electromagnetic induction; ERT: Electric resistivity tomography; FWI: Full waveform inversion; GHG: Greenhouse gases; gx, gy, gz Gravity field components; gzx, gzy, gzz Gradients of gz; L1, L2 Minimization method of robust blocky and least squares normalization, respectively; L-BFGS Limited-memory Broyden-Fletcher-Goldfarb- Shanno of the quasi-Newton technique; m: Cementation factor of Archie equation; n: Saturation exponent of Archie equation; Ω: Ohm, electrical resistance unit; p: Pressure; P: Potential electrode; Φ: Porosity; R: Gas constant = 287.058 J/(kgK); RMS: Root mean square; ρ: Electric resistivity; ρa: Apparent electric resistivity; Sg: Gas saturation; t: time; T: Temperature; TDS: Total dissolved solids; TEM: Transient electromagnetic induction; Tx, Rx: Transmitter and receiver coil of EMI method; Vp Velocity of primary seismic waves; Vs Velocity of secondary seismic waves

One unprecedented challenge facing the human being is the energy resources, and its coupling with global climate change and warming from greenhouse gases, GHG. The anthropogenic GHG stems from burning fossil fuels which has a finiteness supply. Mitigating GHG demands developments of viable alternative of renewable energy resources including hydroelectric, biomass, solar, wind, marine (wave/ tides) and geothermal sources. Most of these sources are vulnerable to the climate vagaries, i.e., intermittent, and need buffer storage, such as a compressed air/gas energy storage (CAES). The daily electricity demand is constantly changing while the surplus renewable generation is put into storage. The basic CAES approach is as follows: When electricity generation is greater than demand, the energy surplus is used to compress air in the storage. When generation is less than demand, pressurized air is withdrawn

from storage and used to drive electricity turbines. Advantages of renewable energy storage are (1) balancing power demand and intermittent renewable energy production, (2) bridging temporal mismatch (disparity) between renewable energy production (off-peaks) and demand (peaks), i.e., storing off-peak energy supply to use it during peak demand periods, and (3) offering large buffer capacity to meet any disruptions in energy supply. The underground geology offers an adequate option for short- and long-term storage (hourly to seasonally) of huge energy quantities such as CAES, e.g., [1,2]. The North German Basin delivers favorable conditions (geological, geochemical) for underground space utilization including CAES in porous media (e.g., depleted gas fields or saline aquifers) and salt caverns (natural and artificial). However, this large scale use of this geologic energy storage induces secondary effects in the geologic subsurface, including pressure increase, migration of reservoir fluids, geochemical and biological changes, and geo-mechanical stresses as well as possible leakage of injected fluids or gases into drinking water aquifers [3]. In geologic energy storages, the gas replaces the pore brine causing strong changes in elastic moduli, density and electric resistivity. These physical contrasts justify the application of integrative geophysical techniques for monitoring this geo-storage. Since some years ago Germany practices a turnaround in the energy policy (German “Energiewende”) and is currently leading in the production of solar and wind energy [4]. Germany decided to phase out nuclear energy generation until 2022, which further accelerates the switch to renewable energy sources. Wind energy is produced mainly at the coastal areas (on-shore and off-shore) of North Germany which is characterized by high wind velocities (**Figure 1**). We started 2012 an interdisciplinary joint research project ANGUS+ dealing with impacts of using geologic subsurface as a thermal, electric or material storage in context with alternative energy resources [3]. Using numerical simulations we show here some results of the feasibility study by applying geophysical techniques of elastic full waveform inversion (FWI), electric resistivity tomography (ERT), transient electromagnetic induction (TEM) and gravity. The goal is to explore the capability of these multi-techniques to map and monitor deep CAES reservoirs and potential CAES leakages in the shallow groundwater aquifer. At first we present briefly basics of the applied geophysical techniques then show and discuss examples of mapping and monitoring results from deep CAES reservoirs and a CAES leakage plume in shallow groundwater, and finally give a summary and conclusion.

**Figure 1:** Average wind velocity in Germany and the study site [5].

Here we review briefly the geophysical methods of FWI, ERT, gravity and TEM applied in this numerical study.

**Seismic Full Waveform Inversion (FWI)**

Classical seismic travel time based tomographic approaches are limited in resolution to the first Fresnel zone. In the (spatial) FWI, the resolution can be significantly enhanced by incorporating phase and amplitude information of the whole recorded wave field (in addition to travel time information of specific waves) in the inversion process. The FWI aim is to minimize the data residuals δu = u^{mod} - u^{obs} between the modeled datau^{mod} and the field data u^{obs} to deduce high resolution models of the elastic material parameters in the underground. To solve this nonlinear optimization problem an appropriate objective function E has to be defined. Similar to Asnaashari et al. [5], the following objective function is used:

(1)

Where the term denotes the L2-norm of the data misfit and a weighted L_{2}-norm of the difference between the model parameters *m *and prior model information* mprior* used for model regularization. The parameter λ_{1} balances the contributions of the data misfit and the model regularization term, while the spatial variable weighting factor W_{m} defines the updated parts of the model during the inversion process. Like Asnaashari et al. [6], the magnitudes of the spatial weighting factors are based on elastic reverse-time migration results to spatially restrict model updates to the storage formation. The objective function Equation 1 can be minimized by iteratively updating the model parameters m_{n} (P-wave velocity V_{p}, S-wave velocity V_{s}, density d) at iteration step* n*, starting with an initial background model *m*_{0}, using the Newton method [7]:

(2)

with the gradient of the objective function Equation 1 with respect to the elastic model parameters

(3)

and H_{m} as the second derivative of the objective function (Hessian). As the explicit calculation of the Hessian in the time-domain requires high computation time, the quasi-Newton L-BFGS (Limited-memory Broyden-Fletcher-Goldfarb-Shanno) technique is used [7,8], where the product of the inverse Hessian with the gradient G_{m} is iteratively approximated by finite differences.

The effective calculation of the time-domain gradient directions with the adjoint method for different model parameterizations are described in, e.g., [9,10]. The step length τ_{n} is estimated by a linesearch satisfying the Wolfe-conditions [7] to assure a fast and accurate convergence of the L-BFGS algorithm. To assure model smoothness, a weak wave number domain filter is applied to the estimated search directions at every iteration step. In time-lapse FWI, the data residuals are modified according to

δu=({u^{mod}(t_{1}) − u^{mod}(t_{0})} − {u^{obs}(t_{1}) − u^{obs}(t_{0})})(4)

where δu denotes the difference between the modeled and the field data at time steps t_{0} (baseline model) and t_{1} [11]. This redefinition of the data residuals leads to a stronger focusing of the model updates at reservoir level.

**ERT technique**

At first we introduce briefly the approach for optimized electrode configurations applied here for surface, borehole and combined surface-borehole surveys. For these surveys ERT data acquisition can be conducted in the tripotential 4-pole configurations α (CPPC, C = current electrode, P = potential electrode), β (CCPP) and γ (CPCP). For an N collinear multi-electrode array, a whole comprehensive data set consists of [(N-1) (N-2) (N-3)/8] independent non-reciprocal quadrupole configurations [12]. The effective comprehensive data set results from excluding the redundant configurations of less stable inversions from the whole set, i.e., γ configurations and those of very large geometric factors [13]. The resulting comprehensive data set is still huge, e.g., a pair of 32 borehole electrodes yield >106 data points. It can map subsurface targets with the highest possible resolution but at very long acquisition times, i.e., poor temporal resolution and high costs. Therefore, we apply here an optimization approach based on the model resolution matrix, e.g., [14,15]. It searches for electrode configurations that maximize the resolution of survey results [16]. We generated optimized data sets of practical sizes of up to 2% only of the comprehensive data set (i.e., ≤ 20,000 data points) but with almost the same spatial resolution. Comparative studies of diverse configurations (standard and non-standard) show the superiority of the optimized configuration results, e.g. [17]. Pore storage formations consist of a highly resistive matrix (e.g., sandstone) and conductive pore brine saturant. The bulk resistivity (p) resulting from the air/gas displacing the brine is predicted by Archie’s law [18]:

Where pbr is the brine/water resistivity, ϕ the porosity, S_{g} the gas saturation and a, m and n are Archie constants. The separate phases (matrix, brine and gas) are assumed without any interaction. The electrical conductivity of the storage formation is caused mainly by the electrolytes of its pore brine. In the North German Basin, temperature (T), pressure (p) and particularly the salinity or total dissolved solids (TDS) increase with depth. Increasing TDS and T causes a dramatic decrease in the resistivity, e.g., [19]; Schlumberger log interpretation charts, 2009. The TDS rise increases the number of ions carrying electrical currents. The T rise increases the salt solubility and decreases the brine viscosity which in turn enhances the ion mobility. The p increase, on the other hand, causes a slight increase in the resistivity due to the closure of cracks. However, this effect decreases with increasing depth and is negligible at pressure >0.3 GPa, e.g., [20]. In formations of the North German Basin, the average vertical gradient (with depth) approaches 100 mg/L/m for TDS of brine, 0.03°C/m for T and 22.6 kPa/m for p, e.g., [21]. In our numerical simulations, ERT surveys are conducted by placing electrode arrays at the surface and in boreholes. The storage targets host these vertical borehole arrays at an aspect ratio mostly of 1, i.e., array lengths are equal to their offset, to achieve a reasonable spatial resolution of the survey results. Sometimes the aspect ratio is increased (up to 2) to increase lateral coverage on the expenses of the resolution and vice versa [17]. The electric bulk resistivity models are parameterized by applying almost realistic values for the input parameters of Archie Equation 5 that are prevailing in the North German Basin. A 2.5D forward and inverse ERT modeling is carried out using modern codes (RES2DMOD, RES3DMODx64 and RES2DINVx64) based on algorithms by, e.g. Loke et al. [22]. The forward modeling code is applied to generate synthetic data sets using optimized electrode configurations for the different surveys. These data sets are generated after CAES injections in the storage reservoir. The data quality (0.6% average simulation error) is confirmed by results of tests on a homogeneous model with a constant p value. In the ERT inversions, diverse setup constraints (mainly regularizations) are applied. These include the minimization methods of least squares (L_{2}) or robust blocky normalization (L_{1}), and initial models of a constant homogeneous resistivity or an approximate inverse model, e.g., [23]. To minimize the ambiguity problems of non-uniqueness in the solution, ERT data are inverted with incorporating mapping data of the subsurface. For geostorage monitoring purposes this information may be well known from prior (baseline, seismic, logs) surveys, e.g., [24].

**Gravity modeling method**

Gravity studies are conducted here by 3D density modeling to test the technique capability for detecting and monitoring CAES plumes and their temporal variations in the underground (reservoirs and leakages). Rock densities depend on the mineral composition, porosity and its content (fluid and gas), p, T deformations etc. The injected CAES in saline reservoirs displaces pore brine and causes a drop of the bulk density which in turn causes a decrease in the gravity components and gradients. The bulk reservoir density (d) of partial gas saturations is given by [25]:

d=(1 − ϕ) + [d (1 − S_{g}) + d_{g} S_{g}](6)

where d_{g}, d_{m} and d_{br} are the gas, rock matrix and brine/water density, respectively, and S_{g}the gas saturation. Dry air density is calculated as a function of depth (i.e., p and T) by:

d_{g}=p/(RT)(7)

Where R is the gas constant=287.058 J/(kgK) for dry air. Obviously dbr depends on the water density, TDS, p and T and is calculated as described in Batzle and Wang [26]. For calculating these densities, we consider almost realistic values for TDS, p and T as a function of depth within the North German Basin. The Stratigraphy and average bulk densities applied in the modeling study at CAES storage and leakage sites are based on borehole measurements and data base of Geotectonic Atlas of north-western Germany [27-29]. These almost realistic density values of rock formations show wide variations from 2050 kg/m^{3} for top Tertiary formations to 2710 kg/m^{3} for the bottom pre-Permian formations, e.g. [21]. We use here the software IGMAS+ (Interactive Gravity and Magnetic Application System) designed for 3D gravity, gravity gradient and magnetic modeling, e.g., [30,31]. The model is extrapolated outside the volume of interest in all directions (about twice the model length) to avoid any edge effects. We calculated the gravity field components (g_{x}, g_{y} and g_{z}) and gradients (g_{zx}, g_{zy} and g_{zz}) before and after CAES injection, respectively, as well as their difference (residual) anomalies (Δg_{x}, Δg_{y} and Δg_{z}). Here we show only the vertical component Δg_{z} maps calculated at the surface which reflect the strongest anomalies with respect to CAES reservoirs and leakages. We will discuss these gravity anomalies resulting from the CAES plumes and the least resolvable anomaly measurable by the modern microgravimeter with accuracy of 3-5 μGal.

**TEM method**

The transient electromagnetic induction (TEM) method can be applied on land, in the air and in the sea. Applying an electric current to a transmitter coil (Tx) and a rapid turn-off ramp (transient EM) induces a primary magnetic field with a strength depending mainly on the Tx-dipole moment and ground conductivity. This field generates Tx eddy currents within ground conductors which decays with time and in turn induces a secondary magnetic field (B). The temporal change in B (dB/dT) is measured as a decaying voltage by a receiver coil (Rx) at different time gates. This response curve can be transferred to the apparent resistivity (p_{a}) (dB/dT) directly proportional to used for a qualitative interpretation in a finite subsurface model. A conductive layer (e.g., saline aquifer) at shallow depths enhances response signals but decreases the penetration depth, and at greater depths enhances the penetration depth and TEM becomes insensitive to the underlying basement. Loops of high dipole moments and low frequencies are used to resolve the subsurface layers at great depths. Compared to surface ERT, TEM has advantages of surveying large areas without coupling problems to the ground, but limitations of the susceptibility to anthropogenic noises. Here we use the AarhusInv code for the forward generation of synthetic 1D TEM soundings resulting from a subsurface model of CAES leakage scenario and their inversions to recover this model [32]. Analogously to ERT, a hydro-geologic model resulting from multiphase flow modeling with the CAES plume are parameterized in a similar way (by using Archie Equation 5) to transfer them to resistivity models. This inversion relies on adaptive iterative fitting the calculated response of 1D models to the synthetic data. The code allows incorporations of constraints on the initial model (e.g., layer parameters) to minimize interpretation ambiguities.

**Monitoring pore CAES reservoir**

Here we present exemplary monitoring results of two growing CAES plumes injected in pores of thin brine reservoirs within the North German Basin. We introduce at first the synthetic site, and then the monitoring results of ERT and gravity techniques.

**Synthetic pore CAES reservoir site**

To approach more realistic modeling scenarios in this numerical study we selected the Wagrien site within the North German Basin to simulate pore CAES reservoirs (**Figure 1**). The site geology is reconstructed from data base of Geotectonic Atlas of Northwest Germany including seismic and borehole data [28,29]. This site has a model block of 29x28x5.5 km3. The subsurface consists of a thick succession of 14 sedimentary layers ranging in age from pre-Permian to Tertiary. Wagrien stratigraphy shows nearly horizontal layering with a gentle anticline fold (**Figure 2**). The succession includes two thin brine reservoirs of porous sandstone (5-30 m thickness), namely the shallow Rhaet formation (1-1.5 km depth) and the deep Quickborn formation (2-2.5 km depth). Both formations are potential pore reservoirs for CAES and are considered in this study. However, the stratigraphy shows an unconformity within the depth rang (0.7-1.0 km) of the shallow potential reservoir, where formations of the Lower Cretaceous, Lias and Rhaet disappear inside the anticline crest. The caprock for these storage formations consists mainly of siltstone, mudstone and salt layers which are nearly impervious representing seals for these potential CAES reservoirs. This CAES is injected in the top Rhaet/Keuper formation at three locations representing horizontal layering, anticline crest and anticline flank (**Figure 2b**, locations 1-3). In the deep Quickborn formation CAES is injected in the anticline crest (**Figure 2b**, location 4).

**Figure 2:** Study site of Wagrien with: (a) 3D block below the Quaternary- Cretaceous overburden and (b) 2D Stratigraphy section perpendicular to the main anticline structure showing locations of CAES reservoirs in the Rhaet/Keuper (1-3) and quickborn (4) formations, [28].

**ERT monitoring results of pore CAES reservoir**

Here we explore the capability of ERT technique to monitor the CAES plume in the Rhaet/Keuper layer below Wagrien site within the North German Basin (**Figure 2b**, locations 1-3). We consider here six different scenarios of varying geologic settings and injection times. The CAES is injected in a nearly horizontal layer, and at the crest and flank of the gentle anticline. These CAES plumes are observed at early, intermediate and late stages of injections. We may note that due to the disappearance of the Rhaet at the anticline crest (unconformity), the CAES plume is simulated in the underlying Keuper layer at this location. The geologic model is transferred into geoelectric model by applying the Archie Equation 5 with almost realistic parameterization, i.e. applying parameter values prevailing in this basin. In this equation we applied for the sandy reservoir 0.2 for ϕ, 0.08 Ωm (TDS about 100 g/l) for p_{br}, and 0.8 for S_{g}. The Archie constants a, m and n are substituted by values of 1, 2 and 2, respectively, which are commonly used in porous sandstone reservoirs [18]. Based on this petrophysical modeling, the CAES injection increases the bulk resistivity of the brine reservoir from 2 Ωm (without CAES) to 4-50 Ωm for saturations 30- 80% [17]. The ERT survey is conducted for this geoelectric model in a series of borehole electrode arrays within 600-1900 m depth range of the main storage targets (reservoir layer with overlying caprock and underlying aquitard). For these anomalous models caused by the injected CAES, synthetic data sets of apparent resistivity are generated using the code RES2DMOD and optimized electrode configurations for all applied CAES scenarios. These synthetic data sets are inverted using the code RES2DINVx64 with the a priori incorporation of the layer interfaces (without the unknown plume interfaces) to reproduce the initial subsurface model. For monitoring purposes, it is more likely to suppose that considerable subsurface information (including layer interfaces) is already available from previous (baseline) surveys such as seismic mapping and well logs. Our ERT data inversions reconstruct directly the true subsurface resistivity tomograms including the study CAES plume, i.e., no model differencing between tomograms with and without gas plume is necessary (unlike seismic). **Figure 3** shows only the best-fitting tomograms of all differently independent inversions using different regularization parameters [17]. Every best-fitting tomogram shows least root mean square (RMS)-errors of <0.5% and iteration number almost of 5, and is optimized with the blocky L_{1} norm for sharp interfaces. This L_{1} norm yields significantly more accurate results (less RMS-error) than L_{2}norm. During inversion iterations, the applied L_{1} norm allows the model resistivity to vary smoothly within each layer and abruptly across the layer boundaries. This is in accordance with our initial (input, true) subsurface resistivity model with abrupt changes at sharp target boundaries (cf., Loke et al. [33]). The low RMS-error value is explained by the good convergence of the synthetic data set toward the final solution. Obviously, our applied ERT technique reconstructs accurately the different single storage targets (the caprocks, aquitards and particularly CAES plumes) for all applied six scenarios. The technique is also capable to characterize these single scenarios from each other which is of an utmost importance for a successful monitoring strategy. The permanently installed (borehole) electrodes help to maximize the reliability of monitoring data. Modeled anomalies minimize the background effect and thus maximize the time-varying response caused here by CAES. A quantitative evaluation confirms the high resolution capability of applied ERT inversion codes as proven by its low values for the region of index and model residual relative to the input model [24]. Also adding random noise (1-5%) to the synthetic data increases the RMS-error values by a factor of 2-9 but slightly decreases the mapping capability of the techniques [16,34]. In conclusion our applied borehole ERT technique with optimized electrode configurations may be capable to monitor all these applied scenarios of thin CAES reservoir developments accurately if the inversion is geometrically constrained by a priori data on subsurface stratigraphy.

**Figure 3:** Geoelectric monitoring of the CAES injected into the Rhaet/ Keuper storages (Figure 2b, locations 1-3) with structures of horizontal layering, anticline crest and flank (columns), and injections stages of early, intermediate and late times (rows). Continuous lines mark interfaces, dots borehole electrodes and red color CAES plumes.

**Gravity monitoring results of pore CAES reservoirs**

Here we have simulated the growth of CAES plume with increasing injection time in the deep pore storage Quickborn under Wagrien site (**Figure 2**, location 4). The temporal development of the CAES plume is simulated by different successive time-lapses of increasing injection time, i.e., volume and depth to the plume bottom. **Figure 4** (top) and **Table 1** shows the depth to the base of the CAES plume and plume thickness at five different time-lapses. Usually the super light CAES propagates radially into the storage formation following the injection induced pressure gradient and accumulates mainly at the top of the anticline structure after the injection stop due to buoyancy effects [35, 36]. As mentioned before a realistic parameterization is considered for this location of the North German Basin. These subsurface scenario models are transferred to 3D density models by applying typical density values of the whole stratigaphy from literature [e.g., 25,28,29]. 3D density of CAES reservoir was calculated from Equation 6 by substituting ϕ by 0.15 and d_{m} by 2600 kg/m^{3}, and calculating d_{g} from Equation 7. By applying the 3D gravity modeling program of IGMAS+, we calculated for these different time-lapses Δg_{z} anomalies relative to the background without CAES. Resulting maximum Δg_{z} amplitudes are listed for all time-lapses (**Table 1**) and only the measurable Δg_{z} anomalies (≥ 5 μGal, accuracy of micro-gravimeters) are shown in **Figure 4**, bottom. Results display the capability of the gravity method to monitor the temporal growth of CAES plume as clear negative anomalies starting from the size of time-lapse C and later. Time lapse C defines the least detectable plume scenario with 15 m thickness at this anticline crest. Early timelapses (**Figure 4**,** A** & **Figure 4B**) have Δg_{z} amplitudes below 5 μ Gal and thus cannot be detected with this technique. Results reflect the continuous increase of the amplitude and size of the mapped plume with increasing the injection time. Distinguishing the successive time-lapses from each other is hardly between C and D but clearly between D and E due to their corresponding maximum Δg_{z} differences of 2.9 and 10.2 μ Gal, respectively (**Table 1**). This good resolution is owed mainly to the large mass deficit resulting from our assumption of the full replacement of the pore reservoir brine by the CAES of 100% saturation and the absence of any temporal effects on the gravity. Obviously, time-lapse measurements do not require most corrections (e.g., free-air, elevation, Bouguer). However, some near surface changes (e.g., fluctuations of the water table) highly affect the gravity readings. Our gravity modeling on a water table at 10 m depth below the study site yields a measurable 5 μ Gal anomaly (microgravimeter accuracy) already for 0.5 m fluctuations only. Therefore, such fluctuations should be observed in wells to remove their effects from gravity readings.

**CAES leakage in groundwater**

Here we apply the different geophysical techniques to study CAES leakage in the near surface groundwater aquifer based on hydrogeological multiphase flow simulation. In the following we present at first the hydro-geological model, then results of potential field techniques (ERT, TEM and gravity) and finally FWI results.

**Multiphase flow model**

We have modeled a typical leakage scenario for CAES plume migrating upwards from a deep pore storage formation into a shallow fresh water aquifer within the North German Basin. In this multiphase/ multi component flow simulation we used TOUGH2-MP (EOS3) and PetraSim as graphical user interface for preprocessing to import the geological structure and post processing to analyze results [37]. Here CAES seeps for 10 years with 1 kg/s leakage rate into 500 m thick Quaternary aquifer consisting mainly of sandy sediments (**Figure 5**). Leaked CAES searches for release pathways along preferentially permeable zones (e.g., faults and buried channels). The upward CAES migration is driven by the lateral pressure in the underlying reservoir, buoyancy and density difference. In the aquifer CAES plume spreads with groundwater flow and accumulates below a surface cover of low permeable aquitard layer. The hydro-geological model after 10 years of the CAES leakage is used here for studying the applicability of the different geophysical methods to detect and map this leakage.

**Mapping CAES leakage by potential field methods**

We present at first mapping results of ERT and TEM techniques, and finish with forward modeling results of gravity technique. The 2D CAES model simulated in the near surface of the North German Basin (**Figure 5b**) is transferred to resistivity model (**Figure 6a**) by applying the petrophysical Archie Equation 5. The same values of water and air saturations resulting from the hydro-geological modeling are considered here. We applied typical values for water mineralization or resistivity (average 25 Ωm) and ϕ-aquifer (0.20) prevailing at shallow groundwater aquifers of North Germany (e.g., [38]). ERT surveys are conducted using surface and borehole electrode arrays at 20 and 10 m intervals, respectively. Synthetic datasets using the optimized electrode configurations are generated for three survey designs conducted from the surface, in boreholes and from the combined surface-borehole. These data sets are inverted without any constraints and with incorporating information about the aquifer away from the unknown CAES leakage. **Figure 6b**-**Figure 6e** shows some examples of the resulting ERT inversion models. Reconstructed models show the capability of the ERT techniques (of different surveys and inversion constraints) to resolve the resistive anomalies of CAES leakage within the conductive aquifer. Obviously, the resolution of the inverted models is enhanced by: (1) applying the combined surface-borehole survey (compared to the single surface or borehole survey, cf., **Figure 6b**, **Figure 6d**, & **Figure 6e**), and (2) considering constrained inversions (compared to unconstrained inversions, cf., **Figures 6b** and **6c**). Hagrey and Petersen [39] found that the mapping resolution of the near surface zone is enhanced strongly by adding borehole electrodes to the surface electrode array. This resolution enhancement (particularly with depth) is even superior to that resulting from applying the optimized electrode configuration instead of standard ones. We can see that the resolution of the surface survey decreases with increasing depth where the deep CAES anomaly show a smeared oversize and reduced amplitude. This is in accordance with the transverse equivalence for this type of a resistive anomaly inside the conductive medium of groundwater aquifer, where the technique is unable to resolve correctly the anomaly resistivity from its thickness (see below). Regarding TEM, this technique is applied on the same leakage model derived from the hydrogeologic CAES model of **Figure 5b**. Using the ERT input model (**Figure 6a**) three different 1D models were digitized representing the separate shallow and deep plumes of CAES leakages and their combination. Here AarhusInv code is used for forward modeling and inversion [32]. For these different models, 1D datasets (vertical electrical soundings) of apparent resistivity (response) curve were generated for ground- and air-based TEM surveys. These datasets are then inverted by posing constraints on the initial model (layer parameters). Prior knowledge on the model layers with the distribution of resistivity highs and lows may be guest from the qualitative interpretation of apparent resistivity curve obtained from the forward modeling. **Figure 6f** shows that the inverted (output) models reproduce their corresponding initial (input) models with a satisfactory misfit. This reflects the capability of TEM techniques (from ground and aero-surveys) to detect such shallow and deep resistive anomalies and to resolve them from each others. As anticipated the resolution decreases with depth, where the deep resistive anomaly shows more smearing effects (larger thickness and lower amplitude) compared to the shallow anomaly. Obviously, both ERT and TEM techniques are able to detect this resistive CAES leakage and to resolve its shallow and deep anomalies within the conductive groundwater aquifer. However, the resolution of these resistive air plumes in this conductive medium is governed by the equivalence principle of the transverse resistance (p_{h}=constant, h=layer thickness), where smearing effects increase the thickness on expenses of the resistivity amplitude, i.e., the amplitude is here underestimated (**Figure 6**). As we have shown applying the constrained inversion minimizes strongly this smearing effect and may help in recovering the true resistivity which in turn may yield more realistic CAES saturations by applying the petrophysical Equation Finally, the 3D gravity technique of forward modeling using IGMAS+ code is applied on 3D leakage model of **Figure 6a**, [30,31]. This leakage model is transformed to 3D density model using Equation 6. For the sandy aquifer we substituted ϕ by 0.20, d_{m} by 2600 kg/m^{3}, d_{g} by 1.23 kg/m^{3}and applied the same values for water and CAES saturations resulting from the multiphase simulation. Resulting Δg_{z} anomaly relative to the background (**Figure 7**) shows a high sensitivity to this shallow leakage in the groundwater aquifer. The anomaly amplitude depends strongly on the air saturation and leakage size. For this leakage, a negative anomaly approaches amplitude of up to 200 μGal which is far higher than the least measurable value of 5 μGal for microgravity. As mentioned before, this good resolution is owed mainly to the large mass deficit resulting from the replacement of the pore groundwater by the super light air and the assumed absence of any temporal effects on the gravity. Time-lapse data for monitoring purposes do not require most corrections but temporal shallow fluctuations of the water table highly affect the gravity readings. In conclusion, the applied 3D gravity technique shows high sensitivity to the mass deficit resulting from the leakage of the gas phase and its saturation changes, but the data should be corrected for temporal fluctuations of the groundwater table.

**Figure 6:** Mapping CAES leakage in groundwater aquifer of North Germany. It shows 2D ERT start model derived from flow model of Figure 5b (a) and reconstructed resistivity models from surface (b, c) borehole (d) and surface borehole surveys (e), as well as 1D TEM sounding models (f) at locations 1-3 (a). The input leakage plume is recovered by unconstrained (b) and constrained inversions (c-e). Models of surface-borehole and constrained inversion show enhanced resolution compared to their other corresponding models.

**Mapping CAES leakage by WFI**

Based on the results of the multiphase flow simulation (**Figure 5b**) changes of the elastic material parameters due to variations of the gas saturation are calculated via an appropriate petrophysical relation. This is based on a Gassmann fluid substitution model [40] assuming a patchy gas-distribution [41], when the gas partially drains the fully water saturated aquifer. The resulting changes of the elastic material parameters are shown in **Figure 8a- Figure 8c**. A synthetic reflection seismic survey along the transect of 6 km length and 130 m depth is computed by solving the 2D isotropic elastic equations of motion with a time-domain finite-difference (FD) technique on a Cartesian grid [42,43]. The reflection seismic acquisition geometry consists of 500 vertical component geophones located at the surface. For the synthetic dataset 100 shots are recorded using a vertical impact source. The source signature is a 20 Hz Ricker wavelet. A free surface boundary condition is assumed on top of the model, while convolution perfectly matched layers [44] are used at all other boundaries. Due to the free surface boundary condition the data is dominated by the Rayleigh wave, which highly increases the nonlinearity of the inverse problem [45]. The synthetic seismic sections are the input data for the FWI. The initial model for the time-lapse waveform inversion at each time-step is the true elastic baseline model before the CAES leakage. No constraints to the time-lapse data, like sequential frequency/offset filtering or time windowing, are applied and all elastic model parameters are inverted simultaneously. The inversion results of the seismic data (**Figure 8d- Figure 8e**) are compared with the true changes of V_{p}, V_{s} and bulk density d_{b} (**Figure 8a- Figure 8c**). Due to the dominance of the Rayleigh wave in the time lapse data only V_{s} model variations could be recovered with some success, while the changes in V_{p} and d_{b}are underestimated. This is one of our preliminary results and further technique refinements are currently under investigation.

Renewable energy resources are intermittent and need buffer storage to bridge time-gap between production and demand peaks. The North German Basin has favorable conditions (geology and geochemistry) and a very large capacity for CAES in saline formations. Geostorage of gas in porous formations causes strong changes in electro-elastic properties and density, and justify applications of multi-geophysical approach. These CAES targets are characterized by increased resistivity and impedance contrasts, and mass deficit. Using numerical simulations we study here the feasibility of the geophysical techniques of ERT, TEM, gravity and FWI in monitoring deep pore CAES reservoirs, and detecting their leakages in shallow groundwater aquifers below synthetic sites of North Germany. We generated synthetic data sets for different subsurface model scenarios of CAES reservoirs and leakages. Datasets are inverted using different constraints on the initial model. Results of this numerical study reflect the capability of our geophysical techniques to detect and monitor CAES plumes and leakages in the groundwater aquifer below the North German Basin. The constrained ERT inversions of optimized data sets acquired from surface, borehole and combined surface-borehole surveys can resolve well the slightly to highly resistive CAES plumes in the deep pore brine reservoir (extremely conductive) and the shallow groundwater aquifer (moderately conductive). Also the capability of ERT technique in resolving shallow leakages is highly enhanced for the combined surface-borehole survey compared to the separate surface and borehole surveys. The permanently installed (borehole) electrodes help to maximize the reliability of monitoring data by minimizing the background effect, i.e., maximizing the timevarying response caused here by CAES. Resistive CAES anomalies in shallow aquifers can be detected even by the TEM technique which is usually more sensitive to conductive anomaly. The applied gravity technique can detect pore CAES reservoirs in the deep underground, and can determine its lower boundary value of detectability, i.e., the least CAES quantity causing anomaly amplitude in the range of the accuracy of modern micro-gravimeters of 3-5 μ Gal. Moreover, the gravity method is sensitive to CAES leakages in shallow groundwater and the anomaly amplitude depends on the leakage size and CAES saturation. Also the FWI technique can map the shallow CAES leakage in groundwater aquifer by anomalies in the reconstructed Δ V_{p}, Δ V_{s} and Δd_{b} tomograms. However, these tomograms contain inversion artifacts and smearing effects related mainly to the dominance of the Rayleigh wave in the data. We may conclude that applied multi-geophysical techniques complement and confirm each other. CAES plumes cause strong impedance and mass deficits, and moderate resistivity highs, i.e., they are more sensitive for gravity and FWI methods. Applying the constrained inversion minimizes strongly the smearing effects and artifacts of inversions, and helps recovering the true electrical and elastic parameters which in turn may yield more realistic CAES saturations by applying adequate petrophysical equations. All results show that the resulting model resolution is highly enhanced by applying our concept of: (1) the multi-geophysical approach, (2) the optimized approach for surveys and data acquisitions, and (3) the constrained inversion with the a priori use of all available data. Currently we focus on refining the tomography technique of elastic FWI and the constrained inversions of all applied techniques. Our goal is to develop an integrative multi-geophysical approach for detecting mapping, monitoring and quantifying CAES and leakages at high resolution. For this purpose we conduct detailed sensitivity and resolution analyses for all these single techniques with respect to: (1) CAES dimensions and saturations, (2) adequate petrophysical approach (parameters and laws), (3) geologic and tectonic settings of storage sites, (4) post injection fracturing, (5) field setups and data acquisition techniques, and (6) various constrained inversion parameters.

Special thanks for E. Auken and A. Christiansen for providing AarhusInv program, and H.-J. Götze and S. Schmidt for providing the software IGMAS+. We thank D. Wehner, and B. Weise for computer work. This study has been carried out within the framework of ANGUS+ research project funded by the German Federal Ministry of Education and Research (BMBF).We thank project partners for efficient collaborations.

- Crotogino F, Mohmeyer K-U, Scharf R (2001) Huntorf CAES: More than 20 years of successful operation. Solution Mining Research Institute Meeting, Orlando, Florida, USA.
- Succar S, Williams RH (2008) Compressed air energy storage: Theory, resources, and applications for wind power. Princeton Environmental Institute, Princeton University, Princeton 81.
- Bauer S, Beyer C, Dethlefsen F, Dietrich P, Duttmann R, et al. (2013) Impacts of the use of the geological subsurface for energy storage: an investigation concept. Environmental Earth Sciences 70: 3935-3943.
- IEA (International Energy Agency) (2013) World energy outlook special report 2013: redrawing the energy climate map. OECD/IEA, Paris.
- Asnaashari A, Brossier R, Garambois S, Audebert F, Thore P, et al. (2013a) Regularized seismic full waveform inversion with prior model information. Geophysics 78: R25-R36.
- Asnaashari A, Brossier R, Garambois S, Audebert F, Virieux J (2013b) Target-oriented time-lapse imaging using FWI with prior model information, 75th EAGE Conference & Exhibition, London, Expanded Abstracts, We P07 06.
- Nocedal J, Wright SJ (2006) Numerical Optimization, Springer, New York.
- Brossier R (2011) Two-dimensional frequency-domain visco-elastic full waveform inversion: parallel algorithms, optimization and performance. Computers & Geosciences 37: 444-455.
- Tarantola A (1986) A strategy for non-linear elastic inversion of seismic reflection data. Geophysics 51: 1893-1903.
- Kohn D, De Nil D, Kurzmann A, Przebindowska A, Bohlen T (2012) On the influence of model parametrization in elastic full waveform tomography. Geophysical Journal International 191: 325-345.
- Denli H, Huang L (2009) Double-Difference Elastic Waveform Tomography in the Time Domain, SEG Technical Program Expanded Abstracts 2302-2306.
- Noel M, and Xu B (1991) Archaeological investigation by electrical resistivity tomography: A preliminary study. Geophysical Journal International 107: 95-102.
- Loke MH, Wilkinson PB, Chambers JE (2010) Fast computation of optimized electrode arrays for 2D resistivity surveys. Computers & Geosciences 36: 1414-1426.
- Stummer P, Maurer H, Green AG (2004) Experimental design: Electrical resistivity data sets that provide optimum subsurface information. Geophysics 69: 120-139.
- Wilkinson PB, Meldrum PI, Chambers JE, Kuras O, Ogilvy RD (2006) Improved strategies for the automatic selection of optimised sets of electrical resistivity tomography measurement configurations. Geophysical Journal International 167: 1119-1126.
- Hagrey SA (2012a) 2D optimized electrode arrays for borehole resistivity tomography and CO2 sequestration modelling. Pure and Applied Geophysics 169: 1283-1292.
- Hagrey SA (2012b) 2D model study of CO2 Plumes in Saline Reservoirs by Borehole Resistivity Tomography. International Journal of Geophysics.
- Archie GE (1942) The electrical resistivity log as an aid in determining some reservoir characteristics. Transactions of the American Institute of Mining Engineers 146: 54-62.
- Arps JJ (1953) The effect of temperature on the density and electrical resistivity of sodium chloride solutions. Transactions of the American Institute of Mining Engineers, AIME 98: 327-330.
- Brace WF, Orange AS, Madden TR (1965) The Effect of Pressure on the Electrical Resistivity of Water-Saturated Crystalline Rocks. Journal of Geophysical research 70: 5669-5678.
- Magri F, Bayer U, Pekdeger A, Otto R, Thomsen C, Maiwald U (2009) Salty groundwater flow in the shallow and deep aquifer systems of the Schleswig-Holstein area (North German Basin). Tectonophysics 470: 183-194.
- Loke H, Acworth I, Dahlin T (2003) A comparison of smooth and blocky inversion methods in 2D electrical imaging surveys. Exploration Geophysics 34: 182-187.
- Claerbout JF, Muir F (1973) Robust modeling with erratic data. Geophysics 38: 826-844.
- Hagrey SA, Strahser M, Rabbel W (2013) Seismic and Geoelectric Modelling studies of parameters controlling CO2 geostorage in saline reservoirs. International Journal of Greenhouse Gas Control 19: 796-806.
- Gasperikova E, Hoversten GM (2008) Gravity monitoring of CO2 movement during sequestration: Model studies. Geophysics 73: WA105-WA112.
- Batzle M, Wang Z (1992) Seismic properties of pore fluids. Geophysics 57: 1396-1408.
- Inselmann R (1985) Dichte-Messungen in der Bundesrepublik Deutschland. Erdöl-Erdgas 101: 76-78.
- Baldschuhn R, Frisch U, Kockel F (2001) Geotektonischer Atlas von Nordwest- Deutschland und dem deutschen Nordsee-Sektor - Strukturen, Strukturentwicklung, Paläogeographie. Geol Jb A153, Hannover.
- Hese F (2012) 3D Modellierungen und Visualisierung von Untergrundstrukturen für die Nutzung des unterirdischen Raumes in Schleswig-Holstein. PhD thesis, University of Kiel, Kiel.
- Götze H-J, Lahmeyer B (1988) Application of three-dimensional interactive modeling in gravity and magnetics. Geophysics 53: 1096-1108.
- Schmidt S, Plonka C, Götze H-J, Lahmeyer B (2011) Hybrid modelling of gravity, gravity gradients and magnetic fields. Geophysical Prospecting 59: 1046-1051.
- Christiansen AV, Auken E (2009) Presenting a Free, Highly Flexible Inversion Code: Near Surface 2009, Dublin, Ireland, Expanded Abstract, B09.
- Loke MH, Dahlin T, Rucker DF (2014) Smoothness-constrained time-lapse inversion of data from 3D resistivity surveys. Near Surface Geophysics 12: 5-24.
- Ramirez AL, Nitao JJ, Hanley WG, Aines R, Glaser RE, et al. (2005) Stochastic inversion of electrical resistivity changes using a Markov Chain Monte Carlo approach. Journal of Geophysical Research 110: 1-8.
- Graupner BJ, Li D, Bauer S (2011) The coupled simulator ECLIPSE – OpenGeoSys for the simulation of CO2 storage in saline formations. Energy Procedia 4: 3794-3800
- Benisch K, Kohn D, Hagrey SA, Rabbel W (2013) Comparative study of CO2 mass quantification using process modeling and geophysical techniques. Expanded abstract, EAGE, 2nd Sustainable Earth Sciences Conference & Exhibition (SES 2013) Pau, France 5.
- PruessK (1991) TOUGH2: A general purpose numerical simulator for multiphase fluid flow. Report LBL-29400, Lawrence Berkeley Nat Lab.
- Kunkel R, Hannappel S, Voigt HJ, Wendland F (2002) Die natürliche Grundwasserbeschaffenheit ausgewählter hydrostratigrafischer Einheiten in Deutschland. Endbericht eines FuE-Vorhabens im Rahmen des LänderfinanzierungsprogrammsWasser und Boden” der Länderarbeitsgemeinschaft Wasser 97.
- Hagrey SA, Petersen T (2011) Numerical and experimental mapping of small root zones using optimized surface and borehole resistivity tomography. Geophysics 76: G25-G35.
- Gassmann F (1951) Uber die Elastizität poröser Medien. Vierteljahrsschrift der Naturforschenden Gesellschaft in Zürich 96: 1-23.
- Mavko G, Mukerji T (1998) Bounds on low-frequency seismic velocities in partially saturated rocks. Geophysics 63: 918-924.
- Virieux J (1986) P-SV wave propagation in heterogeneous media: velocity-stress finite difference method. Geophysics 51: 889-901.
- Levander A (1988) Fourth order finite-difference P-SV seismograms. Geophysics 53: 1425-1436.
- Komatitsch D, Martin R (2007) An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation. Geophysics 72: SM155-SM167.
- Brossier R, Operto S, Virieux J (2009) Seismic imaging of complex onshore structures by two-dimensional elastic frequency-domain full-waveform inversion. Geophysics 74: WCC63-WCC76.