A stable inversion technique has been developed for threedimensional (3D) interpretation of magnetotelluric (MT) data. The inversion method was based on the Gauss-Newton (linearised least-squares) method with smoothness regularisation. Static shifts were also treated as unknown parameters in the inversion. The forward modelling was done by using the staggered-grid finite difference method. The Akaike Bayesian information criterion (ABIC) was applied to the search for the optimum trade-off parameters while the minimising the data misfit, model roughness, and static shifts. The inversion method was then applied to a large volume MT dataset obtained in the Ogiri geothermal area, southwestern Japan. The inversion was performed stably, resulting in a good fit between the observed and computed apparent resistivities and phases. The 3D resistivity structure recovered was generally similar to that obtained by a two-dimensional (2D) inversion in previous work, although the deeper portion of the 3D model seemed to be more realistic than that of the 2D one. The 3D model was also in a good agreement with the geological interpretation of the geothermal reservoirs. These results indicated the necessity for 3D interpretation in geothermal exploration and in other applications in complicated geological environments.