Abstract: Time-lapse magnetotelluric (MT) monitoring of reservoirs requires repeated 3D forward modeling over the same region at different development stages. Due to the large model size and numerous frequencies involved, conventional vector finite element methods based on direct solvers such as LU decomposition suffer from memory consumption that grows rapidly with the number of degrees of freedom, making them impractical for large-scale time-lapse forward modeling. To address this issue, we employ the flexible generalized minimal residual (FGMRES) method as the outer iterative solver, combined with the auxiliary-space Maxwell solver (AMS) as a preconditioner to replace the direct solver. The complex-valued linear system is reformulated into an equivalent real-valued symmetric form, and a 2×2 block-diagonal preconditioner is constructed, in which AMS is applied to each real symmetric sub-block as the inner preconditioner, effectively accelerating the convergence of the Krylov subspace iteration. The algorithm is first validated against a layered model, where the apparent resistivity agrees well with the 1D analytical solution, with a maximum relative error of less than 2.5%. For a 3D model containing a reservoir-related time-lapse resistivity anomaly, convergence is generally achieved within 20 iterations over the tested frequency range, and the iterative results exhibit apparent-resistivity and phase responses consistent with those obtained by the Pardiso direct solver. For the model with 10989659 edges, compared with the Pardiso direct solver, the computational time is reduced from 47196.8 to 9780.2 s, and the peak memory usage is reduced from 320.3 to 122.2 GB, corresponding to reductions of approximately 79.3% and 61.9%, respectively. These results demonstrate that the proposed iterative solution scheme can support the computational demands of multi-vintage time-lapse forward modeling and provides an efficient and viable approach for large-scale 3D time-lapse MT forward modeling of reservoirs.