In some geological formations, borehole resistivity measurements can be simulated using a sequence of 1D models. By considering a 1D layered media, we can reduce the dimensionality of the problem from 3D to 1.5D via a Hankel transform. The resulting formulation is often solved via a semi-analytic method, mainly due to its high performance. However, semi-analytic methods have important limitations such as, for example, their inability to model piecewise linear variations on the resistivity.

more »

... n, we develop a multi-scale finite element method (FEM) to solve the secondary field formulation. This numerical scheme overcomes the limitations of semi-analytic methods while still delivering high performance. We illustrate the performance of the method with numerical synthetic examples based on two symmetric logging-while-drilling (LWD) induction devices operating at 2 MHz and 500 KHz, respectively. The first commercial LWD device appeared in the 1970s. They were commercially used for formation evaluation, especially in high-angle wells. Nowadays, LWD devices are used both for reservoir characterization [6] [7] [8] and geosteering applications [9] [10] [11] . Modern borehole resistivity instruments can measure all nine couplings of the magnetic field, namely xx, xy, xz, yx, yy, yz, zx, zy and zz couplings (the first letter indicates the orientation of the transmitter and the second one indicates the receiver orientation) [13, 14] . Since the depth of investigation of LWD resistivity measurements is limited compared to the assumed thickness of the geological layers, it is common to approximate subsurface models in the proximity of the logging instrument with a sequence of 1D models [5] . In a 1D model, we reduce the dimension of the problem via a Hankel or a 2D Fourier transform along the directions over which we assume the material properties to be invariant [15] [16] [17] . The resulting ordinary differential equations (ODEs) can be solved (a) analytically, leading to a semi-analytic approach after the numerical inversion of a Hankel transform [15, [17] [18] [19] [20] , or (b) numerically, leading to a numerical approach [16] . Solving the ODEs analytically has some major limitations, for example, (a) the analytical solution can only account for piecewise constant material properties and other resistivity distributions cannot be solved analytically, which prevents to properly model, for example, an oil-to-water transitioning zone when fluids are considered to be immiscible; (b) a specific set of cumbersome formulas has to be derived for each physical process (e.g., electromagnetism, elasticity, etc.) anisotropy type, etc.; (c) analytical derivatives of certain models (e.g., cross-bedded formations, or derivatives with respect to the bed boundary positions) are often difficult to obtain and have not been published to the best of our knowledge [13] . Solving the resulting ODEs numerically is also possible. The resulting method also exhibits a linear cost with respect to the discretization size, since it consists of a sequence of independent 1D problems. An example of this approach can be found in [16] , where Davydycheva et al. use a 2D Fourier transform to reduce the dimension of the problem. Then, they employ a highly accurate 1D finite difference method (FDM) to solve the resulting ODEs. This method is relatively simple to code. However, this combined methodology requires a computational cost that is over 1000 times larger than that observed in semi-analytical methods. This occurs due to the elevated number of unknowns required to properly discretize the ODEs. If a common factorization of the system matrix based on a lower and an upper triangular matrix (the so-called LU factorization) would be precomputed for all source positions, the situation would only worsen: one would need a refined grid to properly model all sources, and since the cost of backward substitution is also proportional to the discretization size (as it occurs with the LU factorization), the large computational cost required to perform backward substitution would significantly increase the total cost of the solver. The use of other traditional techniques such as a finite element method (FEM) (see, e.g., [21] ) to solve the resulting ODEs would not alleviate those problems. Additionally, in [16] , the use of a 2D Fourier transform presents another burden on the performance of the solver, since the number of ODEs (and the total solver cost) increases quadratically with respect to the number of integration points of the 1D Inverse Fourier Transform (IFT), while in the case of a Hankel Transform, the cost only grows linearly with respect to the number of integration points on the 1D Inverse Hankel Transform (IHT). The main objective of this work is to overcome the above limitations of both semi-analytic and existing numerical methods by solving each 1D problem (associated with a Hankel mode) using an efficient multi-scale FEM. Additionally, we seek a method with the following features: (a) we can consider arbitrary resistivity distributions along the 1D direction, and (b) we can easily and rapidly construct derivatives with respect to the material properties and position of the bed boundaries by using an adjoint formulation, which allows us to compute numerically the derivatives forming the Jacobian matrix needed by the Gauss-Newton inversion method at (almost) no additional cost. Despite these advances, presently our proposed multi-scale method is slower than the semi-analytic one. However, it is approximately two orders of magnitude faster than a traditional 1D FEM or a 1D FDM, like the one presented in [16] . This speedup is essential for practical applications since one often needs