<p>Galvanic distortion caused by lateral small-scale inhomogeneities near the Earth’s surface distorts the impedance tensor of the magnetotelluric (MT) method. The distorted impedance tensor can lead to a false image of the subsurface electrical resistivity structure in MT inversion. Although galvanic distortion can be estimated in addition to subsurface resistivity in the inversion, it requires the determination of the tradeoff parameters for the regularization terms for the distortion tensor as well as the resistivity structure. One of the most effective approaches for robustly estimating the resistivity structure with correction of galvanic distortion is using the response function of the Network-MT method, which uses metallic telephone lines to measure the electric field. The Network-MT response function is negligibly affected by galvanic distortion owing to the use of dipoles longer than several kilometers. However, the Network-MT response function has a disadvantage in its low resolution to shallow structures compared with the conventional MT response function because of the lack of high-frequency data. The combined inversion of the Network-MT and MT response functions can overcome this limitation. Therefore, we developed a three-dimensional inversion method that simultaneously inverts both of the response functions. The developed method is based on the finite element method and can reproduce the voltage difference of a long dipole in the Network-MT method by summing the voltage differences of its small segments. We implemented a combined inversion method for the tetrahedral and hexahedral meshes. By comparing the calculated and analytical response functions, we confirmed that the developed method computes the Network-MT response function with sufficient accuracy. By applying the developed method to synthetic and real-world datasets, we confirmed that combined inversion correctly estimates subsurface resistivity structures.</p> Graphical Abstract <p></p>