Numerical Simulation of Hot Dry Rock Fracture Monitoring by Time-Lapse Magnetotelluric Method

: The spatial distribution of rock fractures during water injection in hot dry rock (HDR) geothermal is essential information for evaluating geothermal well fracturing. The rock fractures produce anisotropic electrical characteristics, which may facilitate geophysical methods for monitoring fracture distribution. This paper studies the anisotropic magnetotelluric (MT) response in water injection fracturing based on the ﬁnite difference method (FEM). It uses the residual phase tensor (PT) method to calculate the characteristics of fracture strike and dip angle. The feasibility analysis of the technique under a different angle and resistivity conditions in geothermal reservoir investigation is discussed in detail. The results demonstrate that the MT can explain the characteristics of anisotropic resistivity structure changes from the water injection fracturing of HDR at the 5 km depth. The time-lapse MT monitoring method can provide a reliable scientiﬁc basis for fracture distribution direction and spreading characteristics.


Introduction
As a stable and sustainable renewable clean energy, hot dry rock (HDR) geothermal has attracted extensive attention worldwide because of its wide distribution and huge resource reserves [1,2]. The HDR is mainly composed of granite with low permeability and high temperature. Monitoring injected fluid during fracturing, and the characterization of the resulting permeability zone, is essential for developing HDR resources [3]. Time-lapse geophysical methods, including seismic, gravity, and electromagnetic methods, can be used to detect and evaluate the geothermal reservoir base on the differences in subsurface physical properties (velocity, density, and resistivity, et al.). Among them, microseismic is the primary method for the monitoring of hydraulic fracturing, which locates fracture sources via micro-earthquakes generated by fracturing to speculate the spatial distribution of fractures [4,5]. However, information on key parameters such as fluid transport and fracture connectivity cannot be obtained. Magnetotelluric (MT) is a passive-source electromagnetic method based on rock resistivity difference to achieve HDR structure and parameter inversion, which has an essential role in circling target areas and determining geothermal well location [6][7][8].
Time-lapse MT has proved to be an effective tool for monitoring the development of dynamic geothermal reservoirs. The apparent resistivity and phase difference are notable characteristics of many underground injected fluids [3,9,10]. Peacock et al. conducted time-lapse MT monitoring in the HDR area of the COSO geothermal reservoirs by calculating anisotropic azimuth and simulating phase changes due to fracture distribution [11]. Abdelfettah et al. used MT to monitor the Soultz geothermal in France for six months. The analysis concluded that the signal effectively monitored the changes during two water injections and the characteristics of fracture extension [12]. Rosas et al. performed a 3D probabilistic MT inversion of the Paralana geothermal to determine the fluid transport depth and flow direction by analyzing the time-lapse results [13]. Peacock et al. conducted MT monitoring before and after water injection at the Paralana geothermal and concluded that fractures have anisotropic characteristics and that the development angle migrates to the northeast [14].
The above studies proved the reliability of the time-lapse MT method in geothermal. To understand the signal response before and after fluid injection, we generated anisotropic resistivity models for different types of HDR reservoirs based on two-dimensional MT forward modelling. The residual phase tensor (PT) method is introduced, and its fracture sensitivity under different strikes and dips is discussed in detail. The results show that the time-lapse MT monitoring is sensitive to the changes in anisotropic conductivity structure and the differences between before and after fluid injection into the geothermal well. The approach could characterize the fracture distribution and development direction generated during the stimulation of hydraulic fractures in HDR and provide a basis for the temporal and spatial evolution of fluids.

Anisotropic Magnetotelluric Finite Difference Simulation
MT is an electromagnetic bathymetric method proposed by Tikhonov and Cagniard, which uses natural alternating electromagnetic fields as the source for studying electrical structure [15,16]. The general MT utilizes a 10 −4 Hz and 10 4 Hz frequency range. The electromagnetic field of the subsurface medium follows Maxwell's equations, and the harmonic field is obtained by Fourier transform: Conductivity tensor σ can be diagonalized, and its characteristics in three directions can be expressed by three central conductivity [10]. When the properties are expressed as isotropic, σ x = σ y = σ z . Upon adding the rotation matrix and rotating the basic Eulerian space by α S , α D , and α L angles [17], where R z and R x are the rotational quantities of the anisotropy, and α S , α D , and α L are the 'strike', 'dip', and 'slant' angles of the anisotropy (Figure 1). After calculation by a rotation matrix, σ is expressed in full tensor form, but the diagonal coefficients are not equal to the rotated diagonal coefficients, σ x = σ xx , σ y = σ yy , σ z = σ zz . The continuity equation for the current density still holds in the anisotropic, where the conductivity tensor's non-diagonal elements can change the current density's direc tion concerning the electric field [18,19]. Expanding the electromagnetic components in Equation (1), constraining / = 0 and replacing the right term with the full tenso yields the following: Equations (3) and (4) can be simplified into the second-order partial differential equa tions containing only and : (5) is essentially just the individual components of the second-order curl curl equation projected into 2D with the assumption / = 0, including the E-polariza tion mode of the electric field (TE mode) and the H-polarization mode of the magnetic This i discretized to form the coefficient matrix by staggering the order of the electrical indice on the grid and using the differential and appropriate averaging operators for approxi mation [19]. Boundary conditions are added to consider the ill-posed problem. The inter nal assumption of the model is the permeability under vacuum. The internal electric and magnetic components of the boundary are continuous, and the external Dirichlet bound ary conditions can be given. To ensure sufficient attenuation of the induced secondary field at the top boundary, the air layer with low conductivity is placed on the top of the model. The permutation equations contain coefficients and node information, and the un derground electric and magnetic components are obtained by solving the partial differen tial problem with the finite difference method.

Phase Tensor
In MT surveys, the inhomogeneity of the subsurface near conductivity distorts the The continuity equation for the current density still holds in the anisotropic, where the conductivity tensor's non-diagonal elements can change the current density's direction concerning the electric field [18,19]. Expanding the electromagnetic components in Equation (1), constraining ∂/∂x = 0 and replacing the right term with the full tensor yields the following: Equations (3) and (4) can be simplified into the second-order partial differential equations containing only E x and H x : Equation (5) is essentially just the individual components of the second-order curl-curl equation projected into 2D with the assumption ∂/∂x = 0, including the E-polarization mode of the electric field (TE mode) and the H-polarization mode of the magnetic field (TM mode), where A = (σxyσyz−σxzσyy) yz . This is discretized to form the coefficient matrix by staggering the order of the electrical indices on the grid and using the differential and appropriate averaging operators for approximation [19]. Boundary conditions are added to consider the ill-posed problem. The internal assumption of the model is the permeability under vacuum. The internal electric and magnetic components of the boundary are continuous, and the external Dirichlet boundary conditions can be given. To ensure sufficient attenuation of the induced secondary field at the top boundary, the air layer with low conductivity is placed on the top of the model. The permutation equations contain coefficients and node information, and the underground electric and magnetic components are obtained by solving the partial differential problem with the finite difference method.

Phase Tensor
In MT surveys, the inhomogeneity of the subsurface near conductivity distorts the electromagnetic signal of the underlying or deep conductivity structure. As the signal frequency decreases, the inductive effect by the near-surface structure gradually decreases and becomes negligible compared with the conductive response [20]. This phenomenon is known as electric field distortion. Eliminating the effect of electric field distortion can effectively recover the accurate electrical response of the subsurface structure. Caldwell (2004) proposed the phase tensor (PT) instead of the original swift method [21]. Taking the ratio between the real and imaginary part of the impedance tensor, such that where Φ is the phase full tensor, Φ can be expressed as a singular value decomposition of the rotation matrix multiplication in the form where R is the rotation matrix, Φ max and Φ min are the maximum and minimum values of the coordinate invariants, and the angles α and β have transformation processes concerning the full tensor and diagonal symmetry. Formally, where angle β and α are the rotation of the phase tensor asymmetry, and the three coordinate invariants define the phase tensor. Equation (9) has the same form as a square matrix's singular value decomposition (SVD) [22]. Thus, the Φ max and Φ min (maximum and minimum values) are the singular values of Φ. In general, the tensor invariants algebra is represented by trace, skew, and determinant: The invariant is expressed as a first-order function, where the maximum, minimum, and skew are Graphically, we may use an ellipse to represent the second-rank phase tensor ( Figure 2) [23], where (X 1 , X 2 ) is defined as the reference coordinate system of the tensor. The tensor principal axis is the major and minor axes of the ellipse, and the angle α − β is the major axis orientation, where the phase tensor is symmetric β = 0, and the major axis orientation is given by α. Heise et al. introduced the residual phase tensor (RPT) by combining the tensor information from two MT datasets to obtain the difference matrix [24].
where I is the identity matrix, Φ 1 is the pre-injection PT, and Φ 2 is the post-injection PT.
Calculating the three rotational invariants of ∆Φ max , ∆Φ min , and ∆β can indicate the change of the preferred transport direction of the current [11,25].
Energies 2022, 15, x FOR PEER REVIEW 5 of 17 orientation is given by . Heise et al. introduced the residual phase tensor (RPT) by combining the tensor information from two MT datasets to obtain the difference matrix [24].
where is the identity matrix, 1 is the pre-injection PT, and 2 is the post-injection PT. Calculating the three rotational invariants of Δ , Δ , and Δ can indicate the change of the preferred transport direction of the current [11,25].

Phase Tensor to HDR Reservoir Model Testing
This section uses typical HDR reservoir models to simulate the fracture isotropic and anisotropic characteristics by the MT resistivity and phase tensor response. The electrical structural state of the injection process is discussed by simulating fractures containing high conductivity layers to demonstrate that fracture distribution is closely related to electrical anisotropy.
The first HDR reservoir isotropic layer mode is shown in Figure 3, in which the model depth is 15 km, the background resistivity is 1000 Ω·m, and there is a high conductivity layer of 20 Ω·m. (x, y, z) in the MT measurement coordinate system. The impedance tensor is obtained by MT simulation and expressed as a second-order matrix with a non-diagonal equal. The results show that a unit circle represents PT. The primary and minor axis of the ellipse is elongated at 0.1 Hz due to the difference in electrical structure, and the primary axis value is positively correlated with the conductivity. The MT response propagated to the deeper part and the axes narrowed, indicating no structural changes.

Phase Tensor to HDR Reservoir Model Testing
This section uses typical HDR reservoir models to simulate the fracture isotropic and anisotropic characteristics by the MT resistivity and phase tensor response. The electrical structural state of the injection process is discussed by simulating fractures containing high conductivity layers to demonstrate that fracture distribution is closely related to electrical anisotropy.
The first HDR reservoir isotropic layer mode is shown in Figure 3, in which the model depth is 15 km, the background resistivity is 1000 Ω·m, and there is a high conductivity layer of 20 Ω·m. (x, y, z) in the MT measurement coordinate system. The impedance tensor is obtained by MT simulation and expressed as a second-order matrix with a non-diagonal equal. The results show that a unit circle represents PT. The primary and minor axis of the ellipse is elongated at 0.1 Hz due to the difference in electrical structure, and the primary axis value is positively correlated with the conductivity. The MT response propagated to the deeper part and the axes narrowed, indicating no structural changes. orientation is given by . Heise et al. introduced the residual phase tensor (RPT) by combining the tensor information from two MT datasets to obtain the difference matrix [24].
where is the identity matrix, 1 is the pre-injection PT, and 2 is the post-injection PT. Calculating the three rotational invariants of Δ , Δ , and Δ can indicate the change of the preferred transport direction of the current [11,25].

Phase Tensor to HDR Reservoir Model Testing
This section uses typical HDR reservoir models to simulate the fracture isotropic and anisotropic characteristics by the MT resistivity and phase tensor response. The electrical structural state of the injection process is discussed by simulating fractures containing high conductivity layers to demonstrate that fracture distribution is closely related to electrical anisotropy.
The first HDR reservoir isotropic layer mode is shown in Figure 3, in which the model depth is 15 km, the background resistivity is 1000 Ω·m, and there is a high conductivity layer of 20 Ω·m. (x, y, z) in the MT measurement coordinate system. The impedance tensor is obtained by MT simulation and expressed as a second-order matrix with a non-diagonal equal. The results show that a unit circle represents PT. The primary and minor axis of the ellipse is elongated at 0.1 Hz due to the difference in electrical structure, and the primary axis value is positively correlated with the conductivity. The MT response propagated to the deeper part and the axes narrowed, indicating no structural changes.   Figure 4 shows the anisotropic layered model with the inclusion of angular deviations. There are different resistivity and angle variations along the electrical axis, σ 1 = σ 3 = 20 Ω·m, σ 2 =1000 Ω·m. The isotropic layer produces a unit circle change with PT response. However, the elliptical shape shows an increase in the major axis and the angle change when passing through the mutation layer. The high conductivity difference in the x-axis is more sensitive than in the y-axis, so the TE mode has a higher apparent resistivity than the TM. Because of the angle change in y and z directions, electromagnetism is no longer parallel. The phase difference produces an anomaly greater than 45 • , making the PT deviate from the central axis.
Energies 2022, 15, x FOR PEER REVIEW 6 of 17 Figure 4 shows the anisotropic layered model with the inclusion of angular deviations. There are different resistivity and angle variations along the electrical axis, 1 = 3 = 20Ω·m, 2 =1000 Ω·m. The isotropic layer produces a unit circle change with PT response. However, the elliptical shape shows an increase in the major axis and the angle change when passing through the mutation layer. The high conductivity difference in the x-axis is more sensitive than in the y-axis, so the TE mode has a higher apparent resistivity than the TM. Because of the angle change in y and z directions, electromagnetism is no longer parallel. The phase difference produces an anomaly greater than 45°, making the PT deviate from the central axis. The results show that the TM and TE modes' phases and apparent resistivity curves overlap in isotropic layer media. The angular information about fracture distribution can not be observed. The deviation of the electrical axis in the anisotropic model causes the difference in the MT signals, and the TE and TM lag phases have been shown to contain much information related to the electrical angle of the formation. Therefore, using isotropic modelling will overestimate the measurement results. Perkor summarized that the anisotropy characteristics of MT data should be considered in most cases [11].
In Figure 5, an anisotropic layer is embedded into the isotropic space to simulate the fracture distribution. The response in the short period is mainly influenced by the isotropic layer, shown as a unit circle. As the period increases, the fracture produces structural differences, and the change in y-direction decreases the apparent resistivity, as indicated by the increased phase central axis. In addition, the TE is affected by the resistivity in the x-direction, and the phase remains constant. The central axis is larger than the minor axis and points in the y-direction. In low frequency, the interface anomaly no longer affects the surrounding electromagnetic, and the PT returns to the isotropic response. In the anisotropic model, the PT of TM and TE can effectively obtain critical information related to the angle of fracture distribution. The results show that the TM and TE modes' phases and apparent resistivity curves overlap in isotropic layer media. The angular information about fracture distribution can not be observed. The deviation of the electrical axis in the anisotropic model causes the difference in the MT signals, and the TE and TM lag phases have been shown to contain much information related to the electrical angle of the formation. Therefore, using isotropic modelling will overestimate the measurement results. Perkor summarized that the anisotropy characteristics of MT data should be considered in most cases [11].
In Figure 5, an anisotropic layer is embedded into the isotropic space to simulate the fracture distribution. The response in the short period is mainly influenced by the isotropic layer, shown as a unit circle. As the period increases, the fracture produces structural differences, and the change in y-direction decreases the apparent resistivity, as indicated by the increased phase central axis. In addition, the TE is affected by the resistivity in the x-direction, and the phase remains constant. The central axis is larger than the minor axis and points in the y-direction. In low frequency, the interface anomaly no longer affects the surrounding electromagnetic, and the PT returns to the isotropic response. In the anisotropic model, the PT of TM and TE can effectively obtain critical information related to the angle of fracture distribution.

Residual Phase Tensor HDR Fracture Monitoring
The PT method depends on electrical interface changes, and the complex structure caused by fractures misjudges the data. The residual PT monitors the moments before and after water injection and can capture a series of parameters such as the direction of fracture distribution and fluid-related electrical characteristics. Figure 6 shows a simplified model of the subsurface with a depth of 15 km, and the low resistivity presence at 3 km indicated the fracture. The background resistivity is 1000 Ω·m, and the fracture is changed by water injection, showing 200 Ω·m, 150 Ω·m, 100 Ω·m, and 50 Ω·m with time, which does not produce anisotropy. The results show that the residual PT is circular, and the color is expressed as the product of the major and minor axes. The decrease in resistivity indicates that fluid is being injected into the subsurface. Since there is no rotation and change of the electrical axis in any direction, the ellipse axes are equal and elongate with the increase of conductivity. The PT response corresponding to 3 km can be monitored from 0.1 Hz to 1 Hz, which indicates that the residual PT method is sensitive to the conductivity changes caused by water injection.

Residual Phase Tensor HDR Fracture Monitoring
The PT method depends on electrical interface changes, and the complex structure caused by fractures misjudges the data. The residual PT monitors the moments before and after water injection and can capture a series of parameters such as the direction of fracture distribution and fluid-related electrical characteristics. Figure 6 shows a simplified model of the subsurface with a depth of 15 km, and the low resistivity presence at 3 km indicated the fracture. The background resistivity is 1000 Ω·m, and the fracture is changed by water injection, showing 200 Ω·m, 150 Ω·m, 100 Ω·m, and 50 Ω·m with time, which does not produce anisotropy. The results show that the residual PT is circular, and the color is expressed as the product of the major and minor axes. The decrease in resistivity indicates that fluid is being injected into the subsurface. Since there is no rotation and change of the electrical axis in any direction, the ellipse axes are equal and elongate with the increase of conductivity. The PT response corresponding to 3 km can be monitored from 0.1 Hz to 1 Hz, which indicates that the residual PT method is sensitive to the conductivity changes caused by water injection. Figure 7 shows the PT response by varying the resistivity in the three electrical axes orientations. The fracture resistivity with 3 km depth is equal to 10 Ω·m in the x, y, and z orientations, respectively. The resistivity in the other orientation is 50 Ω·m.      Figure 7 shows the PT response by varying the resistivity in the three electrical axes orientations. The fracture resistivity with 3 km depth is equal to 10 Ω·m in the x, y, and z orientations, respectively. The resistivity in the other orientation is 50 Ω·m. The residual PT method allows the acquisition of the maximum electrical axis in the fracture (Figure 7). The fracture shows an anisotropic resistivity response in a certain axis direction during water injection. The PT of the x and y axes are elliptically distributed from 0.1 to 10 Hz, and the major axis points to 0° and 90° of the rotation coordinates. The  The residual PT method allows the acquisition of the maximum electrical axis in the fracture (Figure 7). The fracture shows an anisotropic resistivity response in a certain axis direction during water injection. The PT of the x and y axes are elliptically distributed from 0.1 to 10 Hz, and the major axis points to 0 • and 90 • of the rotation coordinates. The vertical axis anomaly cannot respond to the PT because the electromagnetic component is insensitive to the z-axis. Figure 8 shows the effect of depth on the residual PT. In the background resistivity of 1000 Ω·m, the fractures are buried at a depth of 2 km to 5 km, and the fractures develop in the X direction.

Effect of Anisotropic on the Residual Phase Tensor
Energies 2022, 15, x FOR PEER REVIEW 9 o vertical axis anomaly cannot respond to the PT because the electromagnetic componen insensitive to the z-axis. Figure 8 shows the effect of depth on the residual PT. In the background resistiv of 1000 Ω·m, the fractures are buried at a depth of 2 km to 5 km, and the fractures deve in the X direction. The fractures with different depths can still respond to the residual PT, and the m imum electrical axis in the X direction can be observed. Along with the increase of depth of the anomaly, the frequency of MT signal will decrease, and the resolution anomaly characterization will also decrease. This is reflected in the attribute color va of the observation ellipse.
Different strokes are also sensitive to the response of residual PT. Figure 9 shows anisotropic model, considering the change in the direction of the x and y axes. For modelled scenarios presented in Figures 10 and 11, an observable relationship exists tween the strike and the ellipse orientation, which causes the ellipse's long axis to ori tate itself parallel to the anisotropy strike. Note that with the increase of the strike, ellipse's major axis shortens, and the minor axis extends. When the strike angle is 45°, circle phase has the maximum attribute color, and the residual PT cannot distinguish ellipse deviation angle. There is the same phenomenon on the 45° of the y-axis. Thi because the rotation matrix assimilates the anisotropic difference in the horizontal dir tion, where the major and minor axes are equal at 45°. The y-axis variation has the sa characteristics as the x-axis, indicating that the residual PT ellipse possesses an observ response from 0 to 180°. The ellipse axes were normalized to visualize all observatio using the peacock method, highlighting the angular information [14]. The 45° strike an change is observed in the figure on the right, indicating that the original strike is includ in the ellipse calculation. The fractures with different depths can still respond to the residual PT, and the maximum electrical axis in the X direction can be observed. Along with the increase of the depth of the anomaly, the frequency of MT signal will decrease, and the resolution of anomaly characterization will also decrease. This is reflected in the attribute color value of the observation ellipse.
Different strokes are also sensitive to the response of residual PT. Figure 9 shows the anisotropic model, considering the change in the direction of the x and y axes. For all modelled scenarios presented in Figures 10 and 11, an observable relationship exists between the strike and the ellipse orientation, which causes the ellipse's long axis to orientate itself parallel to the anisotropy strike. Note that with the increase of the strike, the ellipse's major axis shortens, and the minor axis extends. When the strike angle is 45 • , the circle phase has the maximum attribute color, and the residual PT cannot distinguish the ellipse deviation angle. There is the same phenomenon on the 45 • of the y-axis. This is because the rotation matrix assimilates the anisotropic difference in the horizontal direction, where the major and minor axes are equal at 45 • . The y-axis variation has the same characteristics as the x-axis, indicating that the residual PT ellipse possesses an observed response from 0 to 180 • . The ellipse axes were normalized to visualize all observations using the peacock method, highlighting the angular information [14]. The 45 • strike angle change is observed in the figure on the right, indicating that the original strike is included in the ellipse calculation.     In further analysis of the sensitivity of the dip angle to the residual PT method, when rotating the x-axis, the y-axis has the maximum resistivity ( Figure 12). The phase of 0° is elliptical, the increase of dip angle gradually decreases the conductivity response of the yaxis, and the major axis is shortened and shows a circular. Phase attribute color decreases In further analysis of the sensitivity of the dip angle to the residual PT method, when rotating the x-axis, the y-axis has the maximum resistivity ( Figure 12). The phase of 0 • is elliptical, the increase of dip angle gradually decreases the conductivity response of the y-axis, and the major axis is shortened and shows a circular. Phase attribute color decreases with increasing dip angle. Note that it is not normalized, because the ellipse axes are independent of the strike angle variation but contain as much information as possible about the dip angle. Figure 13 shows the dip rotation at a strike of 30 • , which is consistent with Figure 12. The above result shows that the dip does not affect the change in the strike but changes the ellipse axis. Figure 11. The residual phase tensor for different strike angles of the conductivity y-axis.
In further analysis of the sensitivity of the dip angle to the residual PT method, when rotating the x-axis, the y-axis has the maximum resistivity ( Figure 12). The phase of 0° is elliptical, the increase of dip angle gradually decreases the conductivity response of the yaxis, and the major axis is shortened and shows a circular. Phase attribute color decreases with increasing dip angle. Note that it is not normalized, because the ellipse axes are independent of the strike angle variation but contain as much information as possible about the dip angle. Figure 13 shows the dip rotation at a strike of 30°, which is consistent with Figure 12. The above result shows that the dip does not affect the change in the strike but changes the ellipse axis.

Time-Lapse Residual PT Monitoring of Fractured HDR Reservoirs
In this section, other characteristics of the residual PT are analyzed jointly using the coupled field conversion to conductivity model with the addition of a reasonable anisotropic rotation. First, we simulated the coupled injection underground field, and obtained the temperature and pressure distribution for 60 months. The model size is 1 km × 2 km, with intersecting fracture network inside (5th year results in Figure 14). The water injection point is located at 1 km, and 20 °C water is continuously injected at an injection rate of 20 m 3 /h. As temperature and water saturation affect the conductivity, empirical formulas were used to obtain electrical parameters and calculate MT response (Figure 14). For

Time-Lapse Residual PT Monitoring of Fractured HDR Reservoirs
In this section, other characteristics of the residual PT are analyzed jointly using the coupled field conversion to conductivity model with the addition of a reasonable anisotropic rotation. First, we simulated the coupled injection underground field, and obtained the temperature and pressure distribution for 60 months. The model size is 1 km × 2 km, with intersecting fracture network inside (5th year results in Figure 14). The water injection point is located at 1 km, and 20 • C water is continuously injected at an injection rate of 20 m 3 /h. As temperature and water saturation affect the conductivity, empirical formulas were used to obtain electrical parameters and calculate MT response ( Figure 14). For the specific work details of this part, please refer to Bai et al. [26]. Figure 13. The fracture anisotropy model and the residual phase tensor for different dip angles in the 30° strike angles.

Time-Lapse Residual PT Monitoring of Fractured HDR Reservoirs
In this section, other characteristics of the residual PT are analyzed jointly using the coupled field conversion to conductivity model with the addition of a reasonable anisotropic rotation. First, we simulated the coupled injection underground field, and obtained the temperature and pressure distribution for 60 months. The model size is 1 km × 2 km, with intersecting fracture network inside (5th year results in Figure 14). The water injection point is located at 1 km, and 20 °C water is continuously injected at an injection rate of 20 m 3 /h. As temperature and water saturation affect the conductivity, empirical formulas were used to obtain electrical parameters and calculate MT response (Figure 14). For the specific work details of this part, please refer to Bai et al. [26]. Add the anisotropy change of 30° strike to the formation after water injection to obtain the residual PT ( Figure 15). Electrical structure changes and elliptical deviations of 23 to 32° were monitored at 100 Hz. This indicates the anisotropy development, which coincides with the fracture angle caused by injection. Add the anisotropy change of 30 • strike to the formation after water injection to obtain the residual PT ( Figure 15). Electrical structure changes and elliptical deviations of 23 to 32 • were monitored at 100 Hz. This indicates the anisotropy development, which coincides with the fracture angle caused by injection.
The phase ellipse has three invariants and adding attribute colors can display more feature information, including the major, minor axis, and ellipse angle. The α − β value is the rotational of the ellipse, which is directly related to the anisotropic axis or the electric field distortion axis. ∆Φ max and ∆Φ min is the scalar value of ellipse axis, which represents the dimensional change of electrical structure. The ellipse color away from the injection becomes lighter, indicating that the electrical difference decreases. Generally, the horizontal anisotropy is related to the ellipse axes. When the conductivity increases, the major and minor axes elongate.
∆Φ max × ∆Φ min and ∆Φ max /∆Φ min are combined invariants. Previous work defines the former as having the information of structure, and the multiplication of ellipse axes can represent the comprehensive characteristics of horizontal anisotropy. The latter is defined as the elliptic flat invariants or the variance in subsurface horizontal anisotropy. When ∆Φ max /∆Φ min is larger, it means that the difference in the conductivity axes is more prominent. In the isotropic, ∆Φ max ∆Φ min = 1. The phase ellipse has three invariants and adding attribute colors can display more feature information, including the major, minor axis, and ellipse angle. The − value is the rotational of the ellipse, which is directly related to the anisotropic axis or the electric field distortion axis. and is the scalar value of ellipse axis, which repre sents the dimensional change of electrical structure. The ellipse color away from the injec tion becomes lighter, indicating that the electrical difference decreases. Generally, the hor izontal anisotropy is related to the ellipse axes. When the conductivity increases, the major and minor axes elongate.
× and / are combined invariants. Previous work de fines the former as having the information of structure, and the multiplication of ellipse axes can represent the comprehensive characteristics of horizontal anisotropy. The latter is defined as the elliptic flat invariants or the variance in subsurface horizontal anisotropy When / is larger, it means that the difference in the conductivity axes is more prominent. In the isotropic, = 1.

Discussion
The calculation of the PT is essential for solving different anisotropic, correcting electric distortions, and monitoring the subsurface time-lapse distribution. The PT variation response and subsurface structure can provide evidence for subsequent HDR monitoring. The residual PT uses the MT at different moments to calculate the difference, effectively eliminating the natural frequency and regular noise during the survey and monitoring the anisotropic angle of fractures caused by injection. This method has a good application in the time-lapse monitoring of MT. This paper discusses the sensitivity of PT to anisotropy and its application in HDR fracture systems in detail through the model test.
In conclusion on the results, the electrical variation in the subsurface was generally reflected in the anisotropic or isotropic change. The difference lay in the inequality or equality of the three conductivity axes (x, y, z). The residual PT was circular for isotropic media, which was proportional to the ellipse axes. Anisotropic media caused elliptical changes, with the major axis pointing to the maximum conductivity axis. Note that residual PT was sensitive to the horizontal axis (x,y), while the vertical axis (z) had no resolution. The more significant the difference in horizontal anisotropy, the flatter the ellipse.
Generally, there was an angle difference between the subsurface anisotropy and the surface measurement system, which indicated that the MT survey lines were not aligned with the conductivity axis. A rotation matrix was applied to transform the conductivity tensor to obtain the dataset in the measurement system. Figure 16 shows the sensitivity analysis of the change in the strike to the residual PT, where the curves indicate the ratio of the major axis to the minor axis. Here, ∆Φ max /∆Φ min > 1 suggests that the phase is elliptical, and the strike can be observed. Tests found that the residual PT had angular sensitivity over a range of horizontal variations from 0 to 180 • . An angle greater than 180 • was considered to produce a rotation. The residual PT makes an isotropic response when the survey system differs from the conductivity axis by 45 • or 135 • , showing a circle with equal ellipse axes. This phenomenon can be solved by normalization.
the time-lapse monitoring of MT. This paper discusses the sensitivity of PT to and its application in HDR fracture systems in detail through the model test.
In conclusion on the results, the electrical variation in the subsurface wa reflected in the anisotropic or isotropic change. The difference lay in the in equality of the three conductivity axes (x, y, z). The residual PT was circular f media, which was proportional to the ellipse axes. Anisotropic media cause changes, with the major axis pointing to the maximum conductivity axis. Note ual PT was sensitive to the horizontal axis (x,y), while the vertical axis (z) had tion. The more significant the difference in horizontal anisotropy, the flatter th Generally, there was an angle difference between the subsurface anisotro surface measurement system, which indicated that the MT survey lines were with the conductivity axis. A rotation matrix was applied to transform the c tensor to obtain the dataset in the measurement system. Figure 16 shows the analysis of the change in the strike to the residual PT, where the curves indica of the major axis to the minor axis. Here, / > 1 suggests that t elliptical, and the strike can be observed. Tests found that the residual PT h sensitivity over a range of horizontal variations from 0 to 180°. An angle greate was considered to produce a rotation. The residual PT makes an isotropic resp the survey system differs from the conductivity axis by 45° or 135°, showing a equal ellipse axes. This phenomenon can be solved by normalization. There are strict assumptions on the response of the dip angle to the MT shows the sensitivity analysis of the change in the dip to the residual PT. When line is parallel to the axis of maximum anisotropy and is rotated along it, the m the PT does not produce a response. The major axis records the rotating axis with the maximum anisotropy axis. In contrast, the change of the minor axis d the joint action of another horizontal axis and vertical axis. Dip changes are al to residual PT. In Figure 17, when y > z , the major ellipse axis shortens as th angle increases. When z > y , the major axis elongates. Note that the dip cha produce a deviation from the elliptical orientation. There are strict assumptions on the response of the dip angle to the MT. Figure 17 shows the sensitivity analysis of the change in the dip to the residual PT. When the survey line is parallel to the axis of maximum anisotropy and is rotated along it, the major axis of the PT does not produce a response. The major axis records the rotating axis coinciding with the maximum anisotropy axis. In contrast, the change of the minor axis depends on the joint action of another horizontal axis and vertical axis. Dip changes are also sensitive to residual PT. In Figure 17, when σ y > σ z , the major ellipse axis shortens as the tendency angle increases. When σ z > σ y , the major axis elongates. Note that the dip changes do not produce a deviation from the elliptical orientation. The HDR monitoring needs to consider the MT data before injection as background for the residual PT to obtain critical information about fracture d and subsurface structural changes by observing MT after injection. The selecti ence data is crucial, which may be accompanied by an anisotropic electrical st related to fracture. Figure 18 shows the response of the residual PT when cons anisotropic in the subsurface before injection. Assuming that the background low conductivity before injection, the maximum electrical distribution exist horizontal axes. The model after the injection has the same parameter settings 9. The results show that anisotropic background affects the residual PT afte When the strike coincides with the maximum horizontal axis, the / smaller than the anisotropic ( Figure 16). This can be explained by the decrease i difference between the two moments. However, the fracture strike can still b in the results. The HDR monitoring needs to consider the MT data before injection as a reference background for the residual PT to obtain critical information about fracture distribution and subsurface structural changes by observing MT after injection. The selection of reference data is crucial, which may be accompanied by an anisotropic electrical structure unrelated to fracture. Figure 18 shows the response of the residual PT when considering the anisotropic in the subsurface before injection. Assuming that the background possesses a low conductivity before injection, the maximum electrical distribution exists along the horizontal axes. The model after the injection has the same parameter settings as in Figure 9. The results show that anisotropic background affects the residual PT after injection. When the strike coincides with the maximum horizontal axis, the ∆Φ max /∆Φ min will be smaller than the anisotropic ( Figure 16). This can be explained by the decrease in the phase difference between the two moments. However, the fracture strike can still be observed in the results. The HDR monitoring needs to consider the MT data before injection as a reference background for the residual PT to obtain critical information about fracture distribution and subsurface structural changes by observing MT after injection. The selection of reference data is crucial, which may be accompanied by an anisotropic electrical structure unrelated to fracture. Figure 18 shows the response of the residual PT when considering the anisotropic in the subsurface before injection. Assuming that the background possesses a low conductivity before injection, the maximum electrical distribution exists along the horizontal axes. The model after the injection has the same parameter settings as in Figure  9. The results show that anisotropic background affects the residual PT after injection. When the strike coincides with the maximum horizontal axis, the / will be smaller than the anisotropic ( Figure 16). This can be explained by the decrease in the phase difference between the two moments. However, the fracture strike can still be observed in the results.

Summary
This paper presents the numerical simulation study of HDR fracture time-lapse monitoring based on the anisotropic MT method. Typical model test results demonstrate that the residual PT can obtain reliable results in the injection process. The sensitivity and re- Figure 18. Sensitivity analysis of the anisotropy model before injection to the residual phase tensor.

Summary
This paper presents the numerical simulation study of HDR fracture time-lapse monitoring based on the anisotropic MT method. Typical model test results demonstrate that the residual PT can obtain reliable results in the injection process. The sensitivity and response characteristics of the residual PT are discussed in detail by testing the fracture model with different depths and parameters and calculating the time-lapse of the coupling field. The PT ellipse axes and angles are related to the subsurface conductivity and anisotropy angle, respectively. Compared with other geophysical methods, MT can indirectly interpret the distribution of fractures caused by the subsurface injection and critical characterization, such as development direction and dip angle, providing a technical means for the exploitation and subsequent evaluation of HDR geothermal. In the future, 3D-modelling and measured data will be studied to obtain actual site results, providing a basis for geothermal fracturing.