## 2 Forward Modeling Method
### Simulation of the Metal Well Casing  
The current source is simulated by combining individual electric dipoles. Each electric dipole source is orthogonally decomposed into \(\mathbf{X}\), \(\mathbf{Y}\), and \(\mathbf{Z}\) components, and each component is interpolated onto the nearest eight grid edges using cubic splines in the same direction. According to Equation (12), the total current density integrated around a specific edge in the cell is the sum of Equation (13) and all decomposed electric dipole  
Figure 3: Configuration of electromagnetic fields on grid cells: (**a**) the arrangement of the electromagnetic field on the grid cells; (**b**) the distribution of electric field components and magnetic field components.  
sources, where the current density from the sources appears as non-zero elements in \(\mathbf{J_{\mathrm{s}}}\). Equation (17) can be simplified into a large-scale linear system of equations:  
\[\mathbf{A}\mathbf{E}=\mathbf{S}, \tag{18}\]  
In our modeling, the bi-conjugate gradient stabilized method (BICGSTAB) is used to solve Equation (18). The coefficient matrix \(\mathbf{A}\) is a sparse, positive-definite, and symmetric complex matrix. Since most elements in the coefficient matrix are zero, for large sparse matrices, the memory usage can be reduced by storing only the non-zero elements of the matrix. The entire coefficient matrix is stored using Compressed Sparse Row (CSR) storage technology to achieve the assembly of the coefficient matrix.  
In this study, [PERSON]'s method is employed to divide the solution area into multiple layers of grids [62]. The fine grids are used to calculate local details, while the coarse grids are used to calculate global characteristics. Each layer employs the Gauss-Seidel iteration to smooth out high-frequency errors. The multi-grid solver is used as a preprocessor, combined with the BICGSTAB to achieve accelerated convergence. Geological modeling is constructed for the entire well section of the modeling area in forward modeling. To accurately simulate the fracturing anomaly body, grid densification is conducted for the target layer of the perforation intervals. Monitoring points are arranged near the ground projection point of the perforation intervals, enabling precise simulation of the surface response. The calculation efficiency of the model meets the actual interpretation requirements.  
The calculation formula for the observation data \(\mathbf{d}_{obs}\) at any position can be written as  
\[\mathbf{d}_{obs}=\mathbf{P}\mathbf{E}, \tag{19}\]  
where \(\mathbf{P}\) is a sparse matrix used for cubic spline interpolation operations from grid cells to observation positions.  
The cased well model is further simplified by treating it as a line excitation source composed of multiple electric dipoles. Current sources in the controlled-source frequency-domain electromagnetic method include vertical electric dipole sources, horizontal electric dipole sources, and long wire sources. The line source under a long metal casing can be simulated using the line source method. Since the radius of the casing is much smaller than its length, the long metal casing is simplified as a line excitation source, assuming the radius of the line source is infinitesimal, and the current density along the line is constant. A line current source consists of a series of electric dipoles, and the superposition of multiple electric dipole fields can effectively approximate the line source (Figure 4b). In the initial simulations in this study, the metal casing of a horizontal well is simplified into two parts (Figure 4a), a vertical line source and a horizontal line source, to investigate the influence of the simplified horizontal well metal casing on the formation.  
Figure 4: Schematic diagram of well trajectory model decomposition: (**a**) simplification of well trajectory model decomposition; (**b**) complex well trajectory model decomposition.