3-D DC Resistivity Forward Modeling Using the Multi-resolution Grid

We implemented a novel multi-resolution grid approach to direct current resistivity (DCR) modeling in 3-D. The multi-resolution grid was initially developed to solve the electromagnetic forward problem and helped to improve the modeling efficiency. In the DCR forward problem, the distribution of the electric potentials in the subsurface is estimated. We consider finite-difference staggered grid discretization, which requires fine grid resolution to accurately model electric potentials around the current electrodes and complex model geometries near the surface. Since the potential variations attenuate with depth, the grid resolution can be decreased correspondingly. The conventional staggered grid fixes the horizontal grid resolution that extends to all layers. This leads to over-discretization and therefore unnecessary high computational costs (time and memory). The non-conformal multi-resolution grid allows the refinement or roughening for the grid’s horizontal resolution with depth, resulting in a substantial reduction of the degrees of freedom, and subsequently, computational requirements. In our implementation, the coefficient matrix maintains its symmetry, which is beneficial for using the iterative solvers and solving the adjoint problem in inversion. Through comparison with the staggered grid, we have found that the multi-resolution grid can significantly improve the modeling efficiency without compromising the accuracy. Therefore, the multi-resolution grid allows modeling with finer horizontal resolutions at lower computational costs, which is essential for accurate representation of the complex structures. Consequently, the inversion based on our modeling approach will be more efficient and accurate.


Introduction
The direct current resistivity (DCR) method is one of the classical geophysical techniques, which is nowadays widely used in the mineral exploration (Oldenburg et al. 1997;Schoor 2005;Zhang et al. 2015), groundwater (Andrade 2011;Thompson et al. 2012), engineering (Chambers et al. 2014;Lysdahl et al. 2017), and environmental problems (LaBrecque et al. 1996;Rosales et al. 2012). The method is often used both for ground measurements and in the boreholes (Loke et al. 2013). Subsequently, the measured electric potential data are inverted to obtain a resistivity image of the subsurface to a depth depending on the separation distance between electrodes.
The forward modeling is an essential step of any inversion algorithm. Over the past few decades, the 1-D (O'Neil and Merrick 1984;Das and Verma 1980) and 2.5-D (Mundry 1984;Pidlisecky and Knight 2008) modeling of the DCR has been developed and used routinely. The 3-D DCR surveying became feasible with the development of the multi-electrode and multi-channel systems (Dahlin 2016;Loke et al. 2013). Therefore, the data in 3-D case can only be fully exploited using 3-D modeling and inversion algorithms. There is a number of algorithms developed so far for 3-D DCR modeling based on the integral equation method (Schulz 1985;Mendez-Delgado et al. 1999), finite-element method (Li and Spitzer 2002a;Rucker et al. 2006;Ren et al. 2018) and finite-difference method (Dey and Morrison 1979;Scriba 1981;Spitzer 1995;Zhang et al. 1995;Loke and Barker 1996;Zhao 1996;Wang et al. 2000;Wu et al. 2003a;Śebastien Penz et al. 2013). The integral equation approach is efficient to handle the models with a few anomaly bodies. On the contrary, the commonly used finite-element method (FEM) and finite-difference method (FDM) are more suitable to cope with models having an arbitrary number of structures (Li and Spitzer 2002b;Wu et al. 2003b). Despite the efforts to improve the efficiency, available 3-D modeling algorithms are still heavy in terms of computational requirements. Since the inversion requires numerous forward calculations, optimized modeling algorithm naturally leads to an efficient inversion.
The electrical field potentials induced by the current electrodes create a singularity effect around the source positions. In order to cope with it, the secondary potentials' formulation is often used (Lowry 1989). However, for both of the total and secondary field approaches, the simulated electric field potentials generally attenuate with depth, meaning that no significant variations can be expected at larger depth. Therefore, responses can be modeled on a grid using coarser resolution (discretization) in the deeper regions, without loss of accuracy.
The conventional staggered (SG) grid employs rectangular cells to discretize the model, which simplifies the gridding process. However, the conformal SG grid fixes the horizontal resolution, which extends to all depths. This may cause over-discretization of the deeper regions, hence leading to redundant computational requirements. The unstructured grid, generally adopted by FEM, commonly employs the tetrahedral cells to discretize the model, which allows local refinement and roughening to avoid the overdiscretization condition.
We present a new 3-D DCR forward modeling based on a finite-difference multi-resolution (MR) grid approach. The MR grid was initially developed for the electromagnetic forward problem and proved to be an efficient alternative to the conventional structured finite-difference approach (Cherevatova et al. 2018). The MR grid is a simplified implementation of the non-conformal grid, it resembles the approach suggested in the octree scheme of Haber and Heldmann (2007), but limited to only the horizontal resolution's refinement. The MR grid can be derived from a fundamental SG grid by horizontally combining adjacent cells. Thus, MR grid represents a vertical stack of several SG grids (sub-grids) with different horizontal resolutions. In this way, the forward modeling operators developed for SG grid are readily applied for each sub-grid.
The main difficulty in MR modeling is the definition of the forward operators at common interfaces between the adjacent sub-grids. Several approaches to defining differential operators at the interfaces were considered and tested in Cherevatova et al. (2018). Generally, all grid nodes and edges are separated into two groups: 'active' and 'inactive'. Inactive elements are not evaluated in the solution, but represented by their neighboring active elements through interpolation. This allows defining differential operators in the physically correct way at the common interfaces. The only difference between the approaches to handle operators on the interfaces lies in the definition of active and inactive grid elements and the interpolation scheme. As a result, different accuracy might be achieved, depending on the selected approach. Following the conclusion of Cherevatova et al. (2018), we selected the Coarse Active (CA) approach, which is shown to be the most efficient and accurate. Moreover, within the CA approach the discrete divergence operator is the adjoint of the discrete gradient operator, which leads to the symmetry of the coefficients matrix. In this case the linear equations set can be solved efficiently using the Preconditioned Conjugate Gradient method. As a result, the time required for solving the system of equations is nearly linear with respect to the degrees of freedom (DoF). The MR grid allows us to substantially reduce the amount of DoF leading to a significant speed-up of the forward modeling.
The paper is organized as follows. In the methods part, we briefly describe the standard staggered grid formulation of the DCR modeling problem, and its modification to define the problem on the multi-resolution grid. In order to verify our newly developed algorithm, in the following sections we present several synthetic examples to examine the speed and accuracy of the MR solution.

Methods
The DCR forward modeling problem can be viewed in terms of solving the Poisson-type equations for the electrical potentials u in a modeling domain X. The corresponding partial differential equations (PDEs) is: r Á ½rðx; y; zÞruðx; y; zÞ ¼ ÀIdðx À x 0 Þdðy À y 0 Þdðz À z 0 Þ; ð1Þ where rðx; y; zÞ is the electrical conductivity arbitrarily distributed in X; rÁ and r are the divergence and gradient operators, respectively. The right-hand side in Eq. (1) defines the source term, where I is the current intensity, ðx 0 ; y 0 ; z 0 Þ denotes the current electrode's coordinate, and dðÁÞ is the delta function.
In order to solve the PDEs, boundary conditions are required on the boundaries confining X into a finite space. In our approach, we implemented Neumann boundary condition on the top boundary C 0 to prevent the electrical current from flowing out. For the distant side boundaries C 1 , we simply employed the Dirichlet boundary condition (Mufti 1976;Tripp et al. 1984) in our tests.
Since the multi-resolution (MR) grid is based on the conventional staggered (SG) grid, first we briefly describe the DCR modeling on SG grid. Then we explain the implementation of MR grid approach and emphasize its differences.

Staggered Grid Approach
The modeling domain X on the SG is discretized into N x , N y and N z rectangular cells in the x, y and z directions, respectively. The entire modeling domain is divided into N c cells (¼ N x Â N y Â N z ) in total, and each cell holds a constant conductivity value r.
Nodes and edges are construction elements of the primary-grid. N n and N e denote the number of nodes and edges, respectively. Edge-lengths L E and cellvolumes V C are metric elements of the primary-grid. The dual-edges are defined as the lines between the adjacent cell-center points and constitute the dualgrid. Metric elements of the dual-grid include dualcell-volumesṼ C and dual-face areasÃ F , as shown in Fig. 1.
Using conventional SG grid finite-difference formulation, Eq. (1) can be represented in a discrete form as: where J is the source elements vector (N n Â 1) defined on the nodes. The injected current is assumed to be distributed overṼ C of the dual-cell around the 'lucky' node where the current electrode is placed (Scriba 1981;Pidlisecky et al. 2007). The elements of J are: On the left-hand side of Eq. (2), u denotes a vector (N n Â 1) of the unknown potential values defined on the grid nodes. G is the discrete gradient operator matrix (N e Â N n ), which maps from nodes to edges and is derived as: where G is a sparse topology matrix (N e Â N n ) of the first derivatives with AE1 non-zero entries mapping from nodes to edges; L À1 E is a vector (N e Â 1) with elements of reciprocal of the edge-lengths. The discrete divergence operator matrix G y (N n Â N e ) is an adjoint of G, which maps back from edges to nodes. G y can be derived using a topology matrix (N n Â N e ) combined with the metric weighting as well. It is readily verified that the transpose of ÀG works as a reasonable approximation to the topology matrix of G y (Cherevatova et al. 2018), see the spatial illustrations and matrix forms of G and ÀG T in Appendix. Thus, G y is defined as: Figure 1 Staggered grid discretization. Black points and solid lines are nodes and edges of the primary-grid, respectively. Dashed lines are dualedges of the dual-grid. The abbreviations denote: L E , the edgelength; V C , the cell volume (e.g. the blue region);Ã F , the dual-face area (e.g. the yellow face);Ṽ C , the dual-cell-volume (e.g. the red region) Vol. 177, (2020) 3-D DC Resistivity Forward Modeling Using the Multi-resolution Grid 2805 whereṼ C andÃ F are the dual-cell-volumes vector (N n Â 1) and dual-face-areas vector (N e Â 1), respectively. The conductivity parameters, initially defined at cell centers, need to be mapped to the cell edges to evaluate current densities from the Ohm's law. Hence, we need to estimate r which represents the averaged conductivity vector (N e Â 1) on the edges. This can be done using the volume-weighted averaging approach which maps from cells' conductivity to r as: whereŴ is a sparse averaging topology matrix (N e Â N c ) mapping cells to edges, see Appendix for further details. V C and r are the cell volumes and cell conductivity vectors (N c Â 1), respectively. Thus, the derived W is regarded as a volume-weighted averaging operator.
In short, Eq.
(2) can be expressed as a system of linear equations: where b (N n Â 1) represents the source term with boundary conditions; A, i.e. the product of G y diagðrÞG, denotes the derived non-symmetric coefficients matrix (N n Â N n ). However, by multiplying both sides of Eq. (7) with ÀdiagðṼ C Þ, i.e. the corresponding dual-cell-volumes, we obtain a new coefficient matrix e A, which is a symmetric and positive definite matrix, the source term is also changed to e b ¼ ÀdiagðṼ C Þb. Therefore the system of equations can be rewritten in the following form: A natural choice for solving a large symmetric system of linear equations would be one of the Krylov subspace iterative methods. We choose to solve Eq. (8) using the Preconditioned Conjugate Gradient (PCG) method, which takes advantage of the symmetry of e A and is superior to the generalized solvers (Spitzer 1995). As a preconditioner, the Symmetric Successive Over-Relaxation (SSOR) method is employed. The combination (SSOR-PCG) shows a fast and stable convergence in our tests. Later we refer to the left and right-hand sides of Eq. (8) as Au and b correspondingly omitting the tilde.

Secondary Field Formulation
The most significant potential variations appear around the current electrodes, causing the singular solutions at those points. Local refinement of the grid around the source region is one of the options to alleviate the numerical errors, but it requires very fine discretization around the current electrodes positions, and therefore, inevitably increases the computational costs. Alternatively, secondary field approach is commonly applied (Lowry 1989) to tackle this problem.
The total potential u can be decomposed into primary potential u p and secondary potential u s . Primary potential u p is the response of the background model, including the source. We use the homogeneous half-space as the background model, therefore u p can be calculated using the analytical solution (Lowry 1989) as: where r 0 and r Ã 0 are the positions of the actual and mirrored (against the surface) current electrodes, respectively; r b represents the conductivity of the background model.
Since u p is generated from the same source as u, meaning that the equations for the total potential and primary potential have the same right-hand side, we have: where A p is the coefficient matrix of u p that is derived in a similar manner as A, but with r b . By reorganizing Eq. (10), we obtain an expression for the secondary potential: where b s ¼ A p À A À Á u p is the source term of the secondary potentials. Since, u p is exactly evaluated 2806 J. Gao et al. Pure Appl. Geophys. using the analytical solution of Eq. (9), the new righthand side could provide a correction effect to the u s solution.

Multi-resolution Grid Approach
In the following, we demonstrate the implementation of the MR grid for the DCR modeling. MR grid can be represented as a vertical stack of several subgrids. Each sub-grid can be regarded as a standard SG grid. An example of MR grid is shown in Fig. 2a. We implemented the DCR modeling on both SG and MR grids within the ModEMM framework (Cherevatova et al. 2018), based on an object-oriented scheme (in Matlab). Thus, it allows the operations of the SG grid to be extended to the sub-grids in the MR grid.
The horizontal resolution of each sub-grid y ) is determined by the following rule: where i denotes the index of the sub-grid, N f x Â N f y is the finest horizontal resolution used in the whole MR grid, Cs ðiÞ is the Coarseness parameter (a non-negative integer) of the ith sub-grid. N ðiÞ x and N ðiÞ y between adjacent sub-grids can only change by a factor of two. It is noteworthy that Cs ðiÞ is not only limited to increase with depth but also allowed to decrease. Thus, N f x Â N f y is not fixed in the top sub-grid (although, we commonly implement it in this way), which could also be assigned to the lower sub-grids to describe complex geometries.
Based on the above rule, the MR grid can be derived from a fundamental SG grid with a discretization: is the vertical discretization of each sub-grid. Second, for a sub-grid with Cs ðiÞ [ 0, horizontally combine each set of 2 Cs ðiÞ Â 2 Cs ðiÞ cells to form a coarser cell. Following the above rule, the discretization of the MR grid is only controlled by the fundamental SG grid and the parameters: Cs ðiÞ ; N ðiÞ z ; Â Ã i¼1;N sg , which leads to a very simple gridding process. For instance, the MR grid in Fig. 2a with parameters: Cs ðiÞ ; N ðiÞ z ; Â Ã i¼1;2 ¼ ð0; 2; 1; 2Þ is derived from a 4 Â 4 Â 4 SG grid. Alternatively, the discretization can start from the coarsest SG grid by refining the sub-grids, however, we always describe sub-grids relative to the finest sub-grid.

Multi-resolution Grid Elements Classification
Since each sub-grid could be regarded as a standard SG grid, at the interior of each sub-grid, the spatial relations between the grid elements are the same as the SG grid, therefore the forward operators (G y , G and W) have the same structure as Eqs. (4), (5) and (6). However, the main challenge in the MR grid implementation is the definition of the forward operators on the overlapped interface between the adjacent sub-grids, e.g. the red interface in Fig. 2a.
On a common interface, part of the grid elements (nodes and edges) from the finer sub-grid overlap with the grid elements from the coarser sub-grid, which will lead to redundant elements. In addition, due to the changing of horizontal resolution, the remained grid elements from the finer sub-grid without the corresponding partner from the coarser sub-grid will be 'hanging' on the interface. See an illustration of the interface in Fig. 2b.
The nodes and edges in the MR grid are classified into two groups: active and inactive elements. The active elements are involved in the calculation explicitly, whilst the inactive elements are eliminated from the solution. All the nodes and edges at the interior of the sub-grids are classified as active, whilst the choice at the common interface will affect the overall solution accuracy. Three approaches were considered and tested in Cherevatova et al. (2018).
(1) Face Active (FA) case classifies the nodes and edges from the finer sub-grid as active and the elements from the coarser sub-grid as inactive similar to Haber and Heldmann (2007) results and provides operators symmetry; (2) Ghost Faces (GF) case is similar to FA, but it extends one 'ghost' interpolated finer layer into the coarser sub-grid, resulting in a second-order accuracy but lacks symmetry (Horesh and Haber 2011); (3) Coarse Active case (CA) is a reverse to FA that takes the nodes and edges from the coarser sub-grid as active. It has accuracy similar to GF but also maintains symmetry. Therefore, CA case is a preferable approach, which we used to implement the MR DCR modeling. Figure 2b shows the classification of CA case. Later we use N active n and N active e to denote the number of the active nodes and edges on MR grid, respectively.

Multi-resolution Grid Forward Modeling Operators
In order to construct the coefficient matrix A, the discrete gradient operator G on MR grid should be defined. In a similar manner to SG grid, G is represented by its topology matrix G mapping from active nodes to active edges and corresponding metric elements. As it was described before, inside the sub-grids, G is the same as the gradient operator on SG grid. However, we need to deal with the inactive elements on the common interfaces as well, this can be achieved by decomposing the topology operator G into three matrices as: Sparse matrix A coarse (N n Â N active n ) interpolates a set of active nodes to a full set of nodes (both the active and inactive nodes). At the common interfaces, A coarse has the different interpolation entries depending on the following conditions: 1. For an active node, A coarse acts as a self-mapping operator, therefore the corresponding interpolation coefficient in A coarse is 1. 2. For an inactive node from the finer sub-grid overlapping with an active node from the coarser sub-grid, e.g. the inactive and active nodes with label i1 and a1 respectively in Fig. 2b, the corresponding interpolation coefficient from a1 to i1 is 1 as well. 3. For an inactive node located on an active edge, the entries are defined by the linear interpolation. In Fig. 2b, the mapped inactive node i2 is interpolated from the adjacent active nodes a2 and a3. 4. For a 'hanging' inactive node located on a face of a coarser cell, it is linearly interpolated from four surrounding active nodes. For example, the inactive node with label i3 is interpolated by the active nodes with label a4 to a7 in Fig. 2b.
In general horizontally non-uniform case, the interpolation coefficients are derived by the corresponding metric elements.
Full gradient topology matrix G full is represented as a block diagonal matrix (N e Â N n ) with gradient topology operators G Thus, G full is the gradient topology matrix of the full MR grid that maps from the full node set to the full edge set, and the sparse selection matrix R coarse (N active e Â N e ) picks out the rows corresponding to the active edges only.
It is also readily verified that the transpose of ÀG is a good approximation to the topology matrix of the discrete divergence operator G y on MR grid, even around the common interface (Cherevatova et al. 2018), see the spatial illustrations and matrix forms of G and ÀG T in Appendix. Thus, G and G y on MR grid can be derived in a similar manner as Eqs. (4) and (5) for SG grid. However,Ṽ C (N active n Â 1) includes only dual-cell-volumes around the active nodes, L E andÃ F vectors (N active e Â 1) include only edge-lengths and dual-face-areas elements of the active edges, respectively. r (N active e Â 1) is only evaluated on the active edges and can be derived from Eq. (6) as well. The topology matrixŴ maps from the cells set to the active edges set. However, the rows ofŴ, which correspond to the active edges at the common interfaces, include more nonzero coefficients related to the extra cells from the adjacent finer sub-grid, see Appendix for further details.
We take the advantage of CA symmetry, i.e. G y is well approximated by the adjoint of G similarly to SG case, therefore, the derived coefficients matrix A is also symmetric. In the MR case, the system of equations can also be solved efficiently using SSOR-PCG iterative solver, but the degree of freedom (DoF) of the system can be significantly decreased compared to the original SG grid.

Results
In this section, several examples were presented to verify the accuracy and efficiency of our algorithm. All the modeling of the SG and MR grid approaches were computed on a PC with 2.60 GHz CPU and 16GB memory.

Example 1
The first example, we considered a 1-D model with three layers. The layers' resistivity and thickness were ðq; Dh; Þ ¼ ð100 Xm; 30 m; 300 Xm; 30 m; 10 Xm; $ Þ from top to bottom, as shown in Fig. 3. We calculated vertical electrical sounding data for Wenner-Schlumberger configuration. The electrodes spacing adopted a ¼ 20 m and n ¼ 1 to 8 levels (separation between the current and potential electrodes is n Â a).
The model was discretized by an SG grid with resolution 120 Â 120 Â 40 at first, the central region was divided by 5 m uniform-cells, see Fig. 3. An MR grid was derived from the SG grid and decomposed it into three sub-grids with parameters: Cs ðiÞ ; N ðiÞ z ; Â Ã i¼1;3 ¼ ð0; 9; 1; 21; 2; 10Þ, i.e. the horizontal resolution was gradually roughened downwards as 120 Â 120, 60 Â 60 and 30 Â 30. The common interfaces between adjacent sub-grids were located inside layers, as shown in Fig. 3. The apparent resistivity responses (q a ) of SG and MR grid solutions as well as the linear filter method solution (Ingeman-Nielsen and Baumgartner 2006) are shown in Fig. 4a, they all represent the characters of the layered model. Since the differences between them are nearly indistinguishable, we calculated the relative differences of the SG and MR grid solutions against the linear filter method solution, as shown in Fig. 4b. The maximum relative differences of the SG and MR solutions are 0:165% and 0:176%, which are nearly the same, meaning that the solutions are both accurate. The relative difference of the MR grid solution is slightly higher than the SG grid solution, which could due to the coarser discretization of the deeper regions, but it's negligible.

Example 2
Next, we considered a 3-D case to demonstrate the accuracy and speed of the MR grid approach against the SG grid approach. The model with 100 X m background included two anomalous blocks with the same resistivity of 1Xm and size of (a) (b) 60 m Â 60 m Â 40 m, where they were placed at different depths of 20 m and 80 m, as shown in Fig. 5a. The data were calculated for one 600 m poledipole profile with 30 m electrode distance (a) and n ¼ 1 to 8 levels, which includes 21 electrodes in total.
A larger SG grid 176 Â 176 Â 50 was created to fit the longer electrode offsets, and its central part was uniformly discretized with 5 m cells. Since the electric field implies the intensity of the potential variation (E ¼ Gu), the secondary electric fields (E s ) on SG grid generated by the current electrode at ðx; y; zÞ ¼ 0 m was calculated. The distribution of E s is shown in Fig. 5b to demonstrate the behavior of u s . Based on the intensity of E s , the most significant u s variations are observed around the shallow block, which is closer to the current electrode. On the contrast, the perturbations around the deeper block (with the same resistivity and shape) are much weaker. Therefore, it demonstrates a general tendency that u s is gradually losing the sensitivity with increasing the distance from the source. Since the potential variations attenuate with depth, the grid resolution can be decreased correspondingly.
The deeper block was placed in the middle sub-grid, which is one step coarser, see Fig. 5a. The q a pseudo-section of the MR grid solution is shown in Fig. 6a. Due to the closer distance to the source, the shallow block induced much stronger anomaly than the deeper block. The MR and SG grid responses were compared, and the relative difference between them is shown in Fig. 6b. Larger differences are mainly focused at the anomaly areas caused by the coarser discretization of the deeper block (-100 to -50 m of the profile coordinate). However, the maximum relative difference is merely 0:35%.
The u s on MR grid generated by the same current electrode at ðx; y; zÞ ¼ 0 m was calculated, and its distribution is shown in Fig. 7a. It is noteworthy that the contour map of each sub-grid was plotted separately without any interpolated connection between them. The contour lines of u s demonstrate that the sub-grids are well-connected and there are no artifices or interruptions due to the horizontal resolution changing between sub-grids. It means that we have sufficient accuracy at the common interfaces.
The relative difference of u between the MR and SG solutions is shown in Fig. 7b. The higher relative differences were induced by the deeper block with coarser discretization, but does not exceed 1:52%. Moreover, there are no obvious differences at the surface where data were measured. We compared the run-times and memory usage of the forward modeling in the SG and MR approaches. In order to avoid uncertainties, the modeling was tested for 5 times, and the averaged values are presented in Table 1.
The MR grid approach requires less time for solving the system of equations (43.33% of the SG grid), which is nearly linear with respect to the reduction in the amount of DoF. As solving the equations set is the most time-consuming part of the forward modeling, the MR approach allows us to greatly improve the efficiency of the forward modeling. The memory usage is also reduced accordingly.

Example 3
The topography has a profound effect in the DCR modeling. Hence, we present a 3-D model that incorporates a 40 m deep valley with 11 surface slope. The background of the model was 100 X m, and the air cells were assigned with 10 10 X m. Block 1 and 2 with the same size (40 m Â 40 m Â 30 m) and resistivity (500 X m) were placed at depths of 55 m and 65 m, respectively. In addition, block 3 with a bigger size (60 m Â 60 m Â 50 m) and higher resistivity (1000 X m) located at a deeper depth of 85 m. The data of two 400 m dipole-dipole profiles along x ¼ 0 m and y ¼ 0 m directions were calculated. Each profile adopted the geometry factors of a ¼ 40m, n ¼ 1 to 8 levels and included 11 electrodes. The electrode arrays followed the topography and located above the block anomalies. See the model and measurement configurations in Fig. 8.
Fine grid discretization is required near the surface to describe the topography. For the SG grid, the 2:5 m Â 2:5 m Â 2 m cells were employed in the central region of the topography layers. In the deeper layers (46 m z 225 m), the cells' thickness was gradually extended to 5 m. It makes the total discretization of the SG equal to 248 Â 248 Â 70. The MR grid divided the SG grid into four sub-grids with parameters: Cs ðiÞ ; N ðiÞ z ; Â Ã i¼1;4 ¼ ð0; 23; 1; 9; 2; 23; 3; 15Þ. The finest horizontal resolution (248 Â 248) was only used in the topmost sub-grid (0 m z 46 m) to discretize the topography portion, afterwards the horizontal resolution was gradually reduced in the lower sub-grids. Block 1 and 3 located 2D slices (central region, y ¼ 0 m) of the contour plots from example 2, the potential is generated by the positive current electrode at ðx; y; zÞ ¼ 0 m (red cone). Dashed lines show the outlines of the anomaly blocks. Thicker solid lines highlight the common interfaces between the adjacent sub-grids. a u s distribution, the contour map of each sub-grid was plotted by independent sections without interpolated connection between them. b The relative differences of u between the MR and SG solutions The abbreviations denote: DoF, the degrees of freedom; t Fwd , averaged time of forward modeling; t eq , averaged time for solving the system of equations (one source); RAM, maximum memory usage. The last row compares the ratio of the above parameters in percent 2812 J. Gao et al. Pure Appl. Geophys. in the second (horizontal resolution: 124 Â 124) and third sub-grids (62 Â 62) respectively, and the common interface (z ¼ 85 m) between them crossed block 2. The bottom padding layers were discretized with the coarsest horizontal resolution of 31 Â 31. See the MR grid discretization in Fig. 9.
In order to avoid the influences from using an improper background model in the secondary field approach (e.g. using the flat half-space model for this topography case), the total field approach was employed. The q a pseudo-sections of the MR grid response are shown in Fig. 10a, b. As it can be seen from the figures, the high q a values caused by the blocks can be observed, where the anomaly responses of two shallow blocks are more obvious. The topography caused perturbations in q a , which could heavily disturb the anomaly responses, especially for the response of block 3, which is hard to find out. Therefore, the topography effect should be taken into account when interpreting the data to avoid inadequate interpretations.
The MR and SG grid solutions were compared with the finite-element method (FEM) solution, which used unstructured grid (Ren et al. 2017). The relative differences between the SG and MR grid solutions against the FEM solution were calculated and presented in Fig. 10c-f. As one can see, the relative differences of the SG and MR grid solutions against the FEM solution are generally less than 1%. The deviations could mainly come from the discretization approaches with respect to the topography. The topography was discretized by SG and MR grids using a mass of fine rectangular cells to fit the surface in a staircase approximation manner. Contrarily, the tetrahedral cells used by FEM can preferably fit the topography surface, which is more suitable for coping with the topography case. Even so, the overall solutions are comparable. The maximum relative differences of SG and MR grids are 0:99% and 1:13%, respectively. The slightly larger differences occurred at n [ 5 part of MR grid solution, which can be explained by using the coarser discretization for the deeper regions of the MR grid comparing with SG grid. However, the relative differences of SG and MR grids generally resemble each other, thus, the MR grid achieved a proximate accuracy with the SG grid.
We present results of the speed-up test for example 3 in a similar manner as example 2, see Table 2. The fine horizontal resolution required to handle the topography leads to a higher number of DoF in the SG case and naturally caused higher computational demands ( Table 2). The MR grid roughen horizontal resolution downwards to avoid the over-discretization in the deep. Consequently, it reduced the DoF significantly to only 38:37% of the SG grid, and resulted in a smaller forward problem. Comparing run-times and memory usage, the MR grid reduced the computational costs to nearly one third of SG grid case, but without loss of the modeling accuracy. It means that the MR grid is more suitable than the SG grid for solving large modeling problem.

Conclusions
We implemented the multi-resolution (MR) grid approach to 3-D Direct Current Resistivity (DCR) modeling. The MR grid allows to adjust the horizontal resolution with depth, which is more flexible than the conventional staggered (SG) grid with respect to description of the shallow geometries, topography, and measurements configuration.
We applied a Coarse Active (CA) approach, which was carefully studied in Cherevatova et al. (2018), to handle the common interfaces between the sub-grids. Beneficial from the definition of the differential operators on the interface layers, the derived coefficient matrix remains in symmetric form as in the SG grid approach. Preconditioned conjugate gradient method was used as the solver, which provides fast and stable convergence due to the symmetry of the system of linear equations. We presented three synthetic studies to verify the accuracy and efficiency of the MR DCR modeling. The SG grid solution with fine discretization guarantees the accurate response, but leads to higher computational demands. The MR grid allows us to reduce the computational costs (CPU time and RAM usage) substantially without loss of accuracy. This is specifically important for inversion algorithms, which require numerous modeling calculations. Taking into account, that MR code can be parallelized to further accelerate calculations, the MR solver opens perspectives for large inversions, which are in high demand for geophysical problems.
b Figure 10 Pseudo-sections from example 3. a, c, e, b, d, f Represent the pseudo-sections of the profiles along y ¼ 0 m and x ¼ 0 m directions, respectively; a, b are the q a response of the MR grid solution; c, d are the relative differences of the SG grid solution against the FEM solution; e, f are the relative differences of the MR grid solution against the FEM solution Table 2 Comparison of the modeling efficiency between SG and MR grid approaches from example 3. The abbreviations denote: DoF, the degrees of freedom; t Fwd , averaged time of forward modeling; t eq , averaged time for solving the system of equations (one source); RAM, maximum memory usage. The last row compares the ratio of the above parameters in percent