The vertical magnetic dipole source under uniform full space conditions is expanded in the right angle coordinate system by Equation 1 (Yu et al., 2017): H X x , y , z = M z x 4 π r 5 3 ϕ u − 2 / π u 3 + u 2 e − u 2 / 2 H y x , y , z = M z y 4 π r 5 3 ϕ u − 2 / π u 3 + u 2 e − u 2 / 2 H z x , y , z = M 4 π r 5 2 z 2 − x 2 + y 2 ϕ u − 2 z 2 − x 2 + y 2 1 + u 2 2 / π e − u 2 / 2 E x x , y , z = 2 / π M ρ y u 5 4 π r 5 e − u 2 / 2 E y x , y , z = 2 / π M ρ x u 5 4 π r 5 e − u 2 / 2 E z x , y , z = 0 (1) where H denotes the magnetic field strength component, in A/m; E is the electric field strength in m; (x, y, z) represents the spatial position and r is the in and out distance, r = x² + y² + z², in m; ϕ u = 2 / π ∫ 0 u e − u 2 / 2 dt is the probability integral; u = 2 π r / 2 π ρ t × 1 0 7 , and ρ is the uniform full space dielectric resistivity, in Ω-m. M denotes the emitted magnetic moment, M = IS, where I is the current emitted, in A; and S refers to the emitted area, in m².

The finite difference method is used to approximate the continuous variables by a system of different discrete equations containing a finite number of uncertainties. Therefore, it is necessary to establish a different format that preserves the essence of the main problem while remaining accurate. The basic for establishing the different equation is to discretize the variables and then approximate the differential quotient in the different equation and quotient. Thus, the different format is obtained for the different components of the magnetic and electric fields, as shown in Equations 2–7: e x n + 1 i + 1 / 2 , j , k = 2 γ − ∆ t n σ i + 1 / 2 , j , k 2 γ + ∆ t n σ i + 1 / 2 , j , k e x n i + 1 / 2 , j , k + 2 ∆ t n 2 γ + ∆ t n σ i + 1 / 2 , j , k × [ h z n + 1 / 2 i + 1 / 2 , j + 1 / 2 , k − h z n + 1 / 2 i + 1 / 2 , j − 1 / 2 , k ∆ y j ¯ − h y n + 1 / 2 i + 1 / 2 , j , k + 1 / 2 − h y n + 1 / 2 i + 1 / 2 , j , k − 1 / 2 ∆ z k ¯ ] (2) e y n + 1 i , j + 1 / 2 , k = 2 γ − ∆ t n σ i , j + 1 / 2 , k 2 γ + ∆ t n σ i , j + 1 / 2 , k e y n i , j + 1 / 2 , k + 2 ∆ t n 2 γ + ∆ t n σ i , j + 1 / 2 , k × [ h x n + 1 / 2 i , j + 1 / 2 , k + 1 / 2 − h x n + 1 / 2 i , j + 1 / 2 , k − 1 / 2 ∆ z k ¯ − h z n + 1 / 2 i + 1 / 2 , j + 1 / 2 , k − h z n + 1 / 2 i − 1 / 2 , j + 1 / 2 , k ∆ x i ¯ ] (3) e z n + 1 i , j , k + 1 / 2 = 2 γ − ∆ t n σ i , j , k + 1 / 2 2 γ + ∆ t n σ i , j , k + 1 / 2 e z n i , j , k + 1 / 2 + 2 ∆ t n 2 γ + ∆ t n σ i , j , k + 1 / 2 × [ h y n + 1 / 2 i + 1 / 2 , j , k + 1 / 2 − h y n + 1 / 2 i − 1 / 2 , j , k + 1 / 2 ∆ x i ¯ − h x n + 1 / 2 i , j + 1 / 2 , k + 1 / 2 − h z n + 1 / 2 i , j − 1 / 2 , k + 1 / 2 ∆ y j ¯ ] (4) b x n + 1 / 2 i , j + 1 / 2 , k + 1 / 2 = b x n − 1 / 2 i , j + 1 / 2 , k + 1 / 2 − ∆ t n − 1 + ∆ t n + 1 2 × [ e z n i , j + 1 , k + 1 / 2 − e z n i , j , k + 1 / 2 ∆ y j − e y n i , j + 1 / 2 , k + 1 − e y n i , j + 1 / 2 , k ∆ z k ] (5) b y n + 1 / 2 i + 1 / 2 , j , k + 1 / 2 = b x n − 1 / 2 i + 1 / 2 , j , k + 1 / 2 − ∆ t n − 1 + ∆ t n + 1 2 × [ e x n i + 1 / 2 , j , k + 1 − e x n i + 1 / 2 , j , k ∆ z k − e z n i + 1 , j , k + 1 / 2 − e z n i , j , k + 1 / 2 ∆ x i ] (6) b z n + 1 / 2 i + 1 / 2 , j + 1 / 2 , k + 1 = b z n + 1 / 2 i + 1 / 2 , j + 1 / 2 , k + ∆ z k × − b x n + 1 / 2 i + 1 , j + 1 / 2 , k + 1 / 2 − b x n + 1 / 2 i , j + 1 / 2 , k + 1 / 2 ∆ x i × − b y n + 1 / 2 i + 1 / 2 , j + 1 , k + 1 / 2 − b y n + 1 / 2 i + 1 / 2 , j , k + 1 / 2 ∆ y j (7)

The initial moments are selected with reference to the empirical equations proposed by Wang and Hohmann (1993), as shown in Equation 8: t 0 = 1.13 μ 0 σ 1 ∆ 1 2 (8) where μ ₀ is the vacuum permeability; σ ₁ denotes the air conductivity; and ∆1 is the grid step at the return line. With a non-uniform grid, a smaller spatial grid step is used near the return line source, anomalies, and measurement points to ensure computational accuracy and facilitate the extraction of observations. The grid step is gradually increased outward by a certain ratio. The air-ground boundary is treated by setting the air layer, while the Dirichlet boundary conditions are adopted as the other boundaries.

The grid subdivision of the model is shown in Figure 1. Generally, in TEM, the ratio of adjacent grid step sizes is required to be less than 1.2, and that between the largest and smallest grid sizes is limited to 20. However, the method proposed in this study can still achieve stable and reliable results when the ratio of adjacent grid step sizes exceeds 2, and no limitations have been observed on the ratio between the largest and smallest grid step sizes. Due to the allowance for larger ratios between adjacent grid step sizes and no restrictions on the maximum-to-minimum grid step ratio, a large model can be constructed with fewer grids. This makes the method highly suitable for applying Dirichlet boundary conditions (Li et al., 2022).

For the based large ground fixed source device form, the emission source is located at the surface, and an analytical solution as shown in Equations 9, 10 exists when the stratigraphic medium is a homogen geologic model for: ∂ B z ∂ T = I ρ a 3 3 ⁡ erf θ a − 2 π θ a 3 + 2 θ 2 a 2 e − θ a 2 (9) B z = I μ 2 a 3 π θ a e − θ a 2 + 1 − 3 2 θ 2 a 2 erf θ a (10) where I is the emission current, representing the apparent resistivity of the half space; a is the radius of the emission frame, denoting the permeability of the homogeneous half space; and T is the time, counting from the time when the current is switched off, where θ = μ 4 ρ T , erf θ a = 2 π ∫ 0 θ a e − t 2 d t are error functions.

In order to verify the accuracy of the algorithm, a simple model was built with the underground half space as a homogeneous medium with a resistivity of 100 Ω ^−m, the transmitting Transforming coil a 100 × 100 rectangular wireframes, and the supply current being 1A, as shown in Figure 2.

Numerical and analytical solutions were computed separately according to the model.

Figure 3 (left) displays a good agreement between the relationship curves of the numerical solution and the analytical solution, while Figure 3 (right) shows that the deviations between the computed results and the analytical solution are all within 0.01, and the error decreases as the time increases. The results indicate that the numerical solution is in excellent agreement with the analytical solution, demonstrating a high computational accuracy.

The geoelectric numerical simulation model was developed based on the aforementioned theory and the real geological conditions. The model includes a deep saline water layer, as depicted in Figure 4. The depth range of the deep saline water layer is between 950 m and 1150 m below the ground. As the depth increases, the resistivity of the ground also increases gradually. The cap layer above the deep saline water layer can be classified as a high resistance ground layer. Using this geoelectric model, a measuring line was established on the ground consisting of 11 measurement points, with each point spaced 50 m apart. A numerical simulation was conducted for this model, and the response characteristics were derived as depicted in Figure 4.

As can be seen in Figure 5 (left), that the multi-channel curve is horizontally dispersed without any irregularity undulation. The decay curve of inductive potential in Figure 5 (right) indicates that the inductive potentials of all the measurement points decrease progressively over time. At the later stage of decay, a small bulging section of the curve is clearly visible, indicated by the red box. Additionally, the decay rate of the inductive potential curve is diminished, which is a prominent feature of the low-resistive body’s response. The analysis of the apparent resistivity profile in Figure 6 reveals that the apparent resistivity is evenly distributed throughout the layer, with a progressive increase from the top to the bottom. According to the established model of the deep brackish water layer, the interface of the first layer of apparent resistivity is located at 350 m below the ground surface; that of the second layer is at 750 m below the ground surface; and the location of the low-resistive body is between 950 m and 1150 m below the ground surface. Due to the inherent low resistance of the deep brackish water layer, the apparent resistivity steadily declines from 750 m, resulting in a low resistance response in the overburden layer above the deep brackish water layer on the map.

Considering the complexity and diversity of the geological underground situation, where the strata is not possibly distributed completely horizontally, a geoelectric numerical simulation model was developed for the inclined deep saline layer. As shown in Figure 7, the model covers a depth range of 750 to,350 m. By using this geoelectric model, a reference line was established on the ground consisting of 11 points with a spacing of 50 m. Subsequently, numerical simulations were conducted by using this model. The response properties are depicted in Figure 8.

The multi-channel curves in Figure 8 (left) display a progressive shift from horizontal to inclined distribution as time progresses. This change is consistent with the inclined direction of the deep saline layer. On the other hand, the decay curves of the induced potentials in Figure 8 (right) demonstrate a gradual decay of the induced potentials at all measurement points. In the later stage of the decay period, a marginal bulge is observed in the decay curves, which differs from the horizontal model decay curve. Specifically, the bulge position of 11 measurement points is gradually backward, leading to a gradual increase in response time for the low resistance anomaly on the decay curve. This is illustrated in the red box. As opposed to the horizontal model decay curve, the 11 measurement points exhibit a progressive backward shift in their bulge positions. Consequently, the depth of the anomalies increases from left to right, leading to a gradual increase in the onset response time of the low-resistance anomalies on the decay curve, as depicted in the red box. The apparent resistivity profile in Figure 9, which aligns with the horizontal saline layer model, reveals an even distribution of apparent resistivity throughout the layer, rising gradually from the top to the bottom. Affected by the low resistance of the inclined saline layer, the interface of the first layer of apparent resistivity is located 350 m below the ground surface, while that of the second layer is situated 750 m below the ground surface. In contrast to the horizontal deep saline layer profile, the apparent resistivity contour does not move horizontally from a depth of 750 m onwards. Instead, it exhibits a downward tilt from left to right. Similarly, the apparent resistivity contour tilts downwards from left to right at depths of 1,100–1,400 m.

The greenhouse effect caused by the large quantity of CO₂ emission has attracted worldwide attention. Geological storage of CO₂ is an effective method to mitigate its impact. The process of CO₂ geological storage involves the compression of CO₂ into a supercritical fluid (35°C, 11 MPa) and its subsequent storage in media such as oil wells, coal beds, and deep saline aquifers at a depth of 800 m. In this study, the primary objective is to characterize the transient electromagnetic response of CO₂ injected into deep saline aquifers.

Physically, the resistivity of the deep saline layer is small. But as fluid CO₂ is injected into the deep saline layer, the water in the rock pores is replaced by fluid CO₂, which leads to a decrease in the electrical conductivity and an increase in the resistivity of the deep saline layer. There is a significant electrical difference with the deep saline layer where CO₂ has not been injected, and the sequestration process is changing with time. Therefore, it is possible to differentiate them by transient electromagnetic methods and to monitor the degree of CO₂ sequestration simultaneously.

In this study, CO₂ was disposed of in a liquid landfill by storing CO₂ in the form of a supercritical fluid under a low-permeability rock cap layer (e.g., mudstone, shale, etc.) and CO₂ in a “sealed canister”. It is required that the storage geological body, i.e., deep water-bearing saline aquifer, should be buried at a depth of more than 800 m (at which depth CO₂ is in the liquid state), the thickness of the underground saline aquifer should be more than 50 m, and that of the cap layer more than 100 m (Zhang et al., 2009). These key parameters guided the construction of the model, as shown in Table 1: 800–1,000 m below the ground is the saline aquifer, with a resistivity of 10 Ω ·m. Three overlaying layers cover the saline aquifer. The first layer, located at 0–100 m, has a resistivity of 200 m. The second layer, located at 100–500 m, has a resistivity of 150 m. The third layer, located at 500–800 m, has a resistivity of 100 m. The uniform half-space below the brackish water layer is bedrock, with a resistivity of 300 m.

Injecting CO₂ into the underground brackish water layer inevitably results in changes in resistivity due to the flow of supercritical CO₂. In order to simplify the analysis, the complex change of the resistivity of the brackish water aquifer is not considered for the time being. It is assumed that after the injection of CO₂, the resistivity changes to 80 Ω ·m due to the reduction of water saturation of the formation. Along with the injection of CO₂, the area with a resistivity of 80 Ω m gradually extends in all directions, as shown in Figure 10. In total, four CO₂ injection areas were set up to simulate the change of CO₂ volume with the rise of the injection time of CO₂ into the deep saline aquifer. The specific size of the four injection areas is as follows: 300 m × 300 m × 200m, 500 m × 500 m × 200m, 700 m × 700 m × 200m, and 900 m × 900 m × 200 m. Three measurement lines were set up, respectively, in the CO₂ injection body center directly above the center, 300 m to the right side of the center, and 500 m to the right side of the center. Each measurement line consisted of 11 measurement points with a point spacing of 200 m. Through numerical simulation, the transient electromagnetic response characteristics of the different volumes of the CO₂ injection body were explored, which in turn facilitated to monitor the geological storage of CO₂.

From Figures 11–15, it can be seen that the forward calculation of the initial model (i.e., the model without CO₂ injection) yields the apparent resistivity, which demonstrates good stratification and consistent data across the three measurement lines without any irregularity changes. As the quantity of CO₂ in the deep saline layer gradually rises, the contour lines of Measurement line 1 change. The contour lines between 800 and 1,000 m become concave, and the contour lines below 1000 m are affected by the high resistance, taking on an upward-convex shape with dense contour lines. The apparent resistivity in the middle position is higher than that at the two ends, which is a clear sign of a high resistance response. Measurement lines 2 and 3 are not situated on top of the anomaly or far away from it, resulting in a laminar response that is not immediately apparent.

When the volume reaches 500 m × 500 m, the contour of Measurement line 1 changes, displaying the same high-resistance response, but the up-convex and down-concave morphology is more powerful. The contour is dense, and the apparent resistivity in the middle position is higher than that at the two sides, which indicates an obvious high-resistance response characteristic. Due to the increased volume of the anomalies, the anomalous response begins to appear in Measurement lines 2 and 3. When the injection of CO₂ reaches 700 and 900 m in length and width, the anomalous response of Measurement line 1 intensifies, leading to further concavity in the apparent resistivity contours. Conversely, the high resistance more strongly affects the part of the apparent resistivity profile below 1,000 m. The value of the apparent resistivity gradually rises, which is shown by the profile’s color gradually becoming bluer and deeper. Simultaneously, the responses of Measurement lines 2 and 3 also gradually intensify, which can play a role in determining the boundary of the CO₂ injector.

As CO₂ is sequestered into the subsurface, the changes in time, strata movement, and the effects of anthropogenic activities lead to the emergence of fissures in the overlying layers of the deep brackish water layer. These fissures result in the escape of CO₂ outward, causing leakage. The escaped CO₂ may leak upward along the crack channels; and due to the pressure reduction in the overlying strata, its form would also transform from fluid to gas. At this point, CO₂ fills the pore water of the overlying formation, increases the resistivity of the rock layer, and creates a clear electrical difference with the original rock layer, thus allowing to detect CO₂ leakage with TEM.

Based on the information given above about CO₂ leakage, a numerical model was established. The saline aquifer that is 800–1,000 m below the ground has a resistivity of 80 Ω ·m. The overlay contains three layers of saline water: the surface layer, located between 0 and 100 m, has a resistivity of 200 Ω ·m; the second layer, located between 100 and 500 m, has a resistivity of 150 Ω ·m; and the third layer, located between 500 and 800 m, has a resistivity of 100 Ω ·m. This layer, along with the uniform half space below the bedrock, has a resistivity of 300 Ω ·m. The resistivity of the CO₂ leakage area is set to 1,000 Ω ·m. At a distance of 800m, the resistivity is set to 100 Ω ·m. The same layer of brackish water and the uniform half-space below, known as bedrock, have a resistivity of 300 Ω ·m.

Figure 16 displays the CO₂ leakage model which was established based on the relevant parameters mentioned above. Three measurement lines were set up at three different positions: directly above the center of the CO₂ injection body, 300 m to the right side of the center, and 500 m to the right side of the center. Each line consisted of 11 measurement points spaced 200 m apart. The transient electromagnetic response characteristic of CO₂ leakage were investigated through numerical simulation.

It can be seen from Figures 17, 18 that when there is no CO₂ leakage, the apparent resistivity profile displays a very good laminar distribution. As CO₂ leaks out, a clear high-resistance region appears, with the strongest high-resistance response right above the leakage location. The high-resistance irregularity response gradually decreases as it moves away from the leakage location. In this case, the location of the CO₂ leakage channel can be determined clearly.