The accuracy and efficiency of 3D forward modelling directly restrict 3D inversion. 3D forward modelling can be divided into direct and indirect methods. The indirect method is based on the time-frequency transformation technique (Newman et al., 1986), which transforms the frequency-domain response to the time domain. The integral equation (Cox et al., 2010), finite difference (Liu and Yin, 2014; Hoversten et al., 2015; Sun et al., 2021), finite-element (Li et al., 2016; Cai et al., 2017; Li et al., 2018; Qi et al., 2019), finite-volume (Jahandari and Farquharson, 2014; Ren et al., 2017; Liu et al., 2018; Zhou et al., 2018), and other numerical methods are applied to calculate the frequency-domain EM response. Widely used frequency-time transformation techniques include the sine (cosine) transform method, the GS inverse Laplace transform method, and the Hankel numerical filtering method (Knight and Raiche, 1982; Guptasarma, 1982; Chen Xiangbin et al., 2008; Sasaki and Cho, 2011; Yin et al., 2013; Li et al., 2016; Li et al., 2018). However, when using the frequency-time transformation method to calculate the TEM response, it is necessary to cover as wide a frequency range as possible. Furthermore, this leads to a large number of calculations and directly affects the calculation efficiency. In addition, precision differences in the calculations exist among these transformation methods, especially with complex geoelectric models, and the transformation methods will become unstable. A full-time TEM response with an arbitrary waveform has an inherent singularity at t = 0. To avoid this singularity, Yin et al. (2013) proposed a convolution operation using the step current and the derivative of the current. However, this requires a high sampling rate of time, further leading to a loss of calculation accuracy.

In terms of direct methods, some breakthroughs have been achieved with the development of the explicit time-domain finite-difference (FDTD) (Sun et al., 2013; Lu et al., 2021), implicit time-domain finite-element (FETD) (Yin et al., 2016; Li et al., 2020b), and other numerical methods. Wang and Hohmann (1993) used a combination of Yee’s staggered grid and a relatively stable and improved Du Fort-Frankel difference method to discretise the Maxwell equation in the time domain directly and alternately determined the electric and magnetic fields. Commer and Newman (2004) applied parallel computing to a TEM response solution based on Wang and Hohmann (1993). Owing to the strong versatility of the FDTD method, it does not need to solve large linear equations, and is easy to implement in parallel. Many scholars have applied this method to study 3D electromagnetic forward modelling (Sun et al., 2013; Commer et al., 2015; Lu et al., 2016). However, Wang and Hohmann (1993) and Commer and Newman (2004) demonstrated that the explicit method requires a strict selection of the time step to obtain a stable and accurate solution. Thus, the forward modelling efficiency is limited. In addition, because FDTD is based on a regular hexahedral mesh for discretisation, this type of mesh restricts the simulation of the TEM response with complex geoelectric models. Compared with the explicit method, the implicit method has a stable solution and is less limited by the time step, but it needs to solve a large linear equation system to obtain the TEM response of each time channel, so it requires high computer memory.

With improvements in computer performance and the development of solving techniques, the application of implicit methods has gradually increased (Movahhedi et al., 2007; Um et al., 2011;Lu et al., 2018; Zhou et al., 2018). Research on implicit methods is represented by the EM inversion team headed by Oldenburg and Haber. Oldenburg et al. (2013) combined a parallel direct solver with an implicit finite-difference space discretisation method to solve a multi-source TEM response. (Haber et al., 2007) developed codes for modelling the TEM response by combining the backward Euler scheme with the finite-volume method. Haber et al. (2007) applied the parallel method to improve the computational efficiency and realised 3D TEM forward modelling and inversion. Yin et al. (2016), Li et al. (2020) referred to Haber et al. (2007) developed the TEM forward modelling based on the unstructured finite-element method and backward Euler scheme, which is flexible to modelling complex earth model TEM response. Recently, another direct algorithm, the direct solution method based on the time-domain Krylov subspace ( Druskin et al., 1999 ), which is different from explicit and implicit methods, has been applied to 3D TEM forward modelling. This method starts with the governing equations in the time domain and independently solves the electromagnetic response at any time in the Krylov subspace. Liu et al. (2018) implemented the Krylov subspace finite-volume method to solve the semi-airborne TEM response and verified its effectiveness.

The integral equation and finite-difference methods are generally based on a regular hexahedral mesh, which is difficult to adapt to complex models. The finite-element method and finite-volume method can fit complex models through deformed hexahedral or tetrahedral meshes, but it is usually assumed that the EM field in each element changes linearly or constantly, so it is necessary to perform fine meshing to improve the accuracy. The detailed mesh subdivision directly affects the forward and inversion efficiencies.

To date, research on simulating dynamic changes in hydraulic fracturing monitoring based on conventional 3D TEM forward modelling is particularly limited. Generally, conventional modelling techniques (detailed mesh subdivision) and forward modelling methods (e.g., finite-element method, finite-volume method) have been used. Hoversten et al. (2015) implemented FDTD to conduct a 3D TEM forward modelling test for fracturing monitoring. For the isotropic thin layer model, a grid with a minimum size of 12.22 mm was used for discretisation in the fracture area, and the number of calculations was substantial. Curcio (2021) studied the TEM signals of a 2D anisotropic earth model based on the finite-element method applied to fracturing electromagnetic monitoring. Hoversten and Schwarzbach (2020) applied a 3D OcTree finite-volume code for modelling borehole-to-surface TEM responses with especially fine meshing around the wells and fractures. Liu et al. (2020) conducted 3D forward modelling for a controlled source electromagnetic method based on a mixed numerical method combining the integral equation method and finite-element method, assuming that the conductivity change caused by fracturing was uniform and did not match the actual situation.

The physical properties of the hydraulic fracturing region, including the fracturing media, considerable burial depth, small dimensions, and apparent anisotropy, cannot be ignored. The burial depth of the hydraulic fracturing target layer is generally greater than 2 km, and the width of each fracturing layer is several metres. The thickness of the SRV is within several centimetres. Additionally, there are high-conductivity casings with a diameter of only several centimetres in the fracturing formation. The multiscale and complex characteristics of the fractured formations cannot be ignored. Therefore, a detailed mesh subdivision containing a large depth scale and small fracture scale should be clearly described, or a 3D multi-scale TEM forward modelling applicable to fracturing monitoring should be adopted. Figure 1 show the mesh subdivision for the traditional earth modelling and multi-scale modelling, we can clearly see that the multi-scale mesh subdivision just uses small dimensional mesh for the small fracture scale, and applies large dimensional mesh for other modelling area.

For the calculation of multiscale physical problems, the solution of the coefficient matrix often consumes huge computing resources via a conventional algorithm (e.g., finite-difference, finite-element), and it is difficult to achieve optimal results, even with traditional iterative solutions. In recent years, the discontinuous Galerkin method has been studied (Angulo et al., 2015; Marras et al., 2016; Xu et al., 2018; Zhan et al., 2018) to solve the problem of multiscale physical model simulations, with increasing applications in fluid mechanics, elastic mechanics, computational electromagnetics, and other fields. In contrast to the traditional Galerkin method, this method adopts a similar approach to the finite-volume method to define continuous numerical fluxes of tangential field components at each interface and allows the solution space of discrete elements to be completely discontinuous at the boundary surfaces between adjacent elements. Furthermore, it can increase the flexibility of the conventional finite-element method. Therefore, the discontinuous Galerkin method effectively partitions and decomposes the overall computational domain and allows the solution of multiscale physical problems with large-scale degrees of freedom. Thus, the discontinuous Galerkin method provides a new technique for the multiscale TEM forward modelling problem for applications in fracturing monitoring.

It is generally assumed that the fracturing formation is isotropic in the analysis of TEM signals for fracturing monitoring. However, anisotropy is widespread (Liu et al., 2018). Different levels of stress on fractures in each direction cause differences in the extension of microscopic fractures. When the scale of the fracturing layer is sufficiently large, the cumulative effect of multiple fracturing causes prominent anisotropic characteristics in the macroscopic formation. This has been confirmed in microseismic monitoring and experimental studies of petrophysical electrical properties (Thiel and Stephan, 2017). It is vital to study the TEM anisotropic response characteristics and effectively identify these characteristics during the fracturing monitoring process. Therefore, electrical anisotropy is also a key factor in TEM forward modelling of hydraulic fracturing monitoring applications.

Additionally, the fracture system is an essential channel for oil and gas storage and transportation. Hydraulic fracturing opens and reforms oil and gas migration channels by constructing artificial fractures. It can then achieve efficient development of unconventional oil and gas. Li et al. (2020a) studied the electrical properties of fractured sandstone with a slightly increased equivalent dielectric conductivity model, indicating that small-scale fractures have a significant impact on the dielectric properties of sandstone. Shen et al. (2009) demonstrated that electrical anisotropy increases with increasing fracture density; in particular, when the formation pressure changes, the anisotropy characteristics of the fracturing formation change significantly. Yan et al. (2020) established an equivalent dielectric conductivity model for sloping anisotropic fracturing rocks, building a theoretical bridge between microfracturing anisotropic media and fracturing anisotropic formations. The scale of microscopic fractures in the fracturing layer is at the millimetre level, and the comprehensive action of different types of fluids and their concentration parameters in the fracturing rock determine the electrical characteristics of the hydraulic fracturing formation (Ahmadian et al., 2018). The influence of rock composition and microstructures on the electromagnetic response is typically ignored in conventional EM forward modelling (Hu et al., 2019). However, considering these physical properties is essential for evaluating the degree of fracture development and determining the fracturing fluid direction. Therefore, forward modelling using a macroscopic equivalent medium model that considers the physical properties of the microscopic fracturing medium is the premise and guarantee for obtaining a highly accurate TEM response of the fracturing zone.

Notably, the signals measured by the TEM-based fracturing monitoring contain feedback information of the entire subsurface medium (high-conductivity casing, non-fracturing formations, and fracturing formations), and the effect of the high-conductivity casing on the measurement signal cannot be ignored. In the process of multiple fracturing, the electrical characteristics of the fracturing stratum are often time shifted. In contrast, the physical properties of the casing and non-fracturing strata remained unchanged. Therefore, considering the non-fracturing formation and high-conductivity casing as the background conductivity and using the EM signal difference value between two adjacent multi-stage fracturing processes, the EM signal variations of the fracturing formation are extracted. Hu et al. (2022) proposed an efficient numerical approach to manage the discretion of well casings by assigning lumped conductive properties to the mesh edges.