Abstract: Time-lapse magnetotelluric (MT) monitoring of reservoirs requires repeated three-dimensional (3D) forward modeling over the same region at different development stages. Due to the large model size and the large number of 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 minimum residual method (FGMRES) as the outer iterative solver, combined with an auxiliary-space Maxwell solver (AMS) 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 and phase show good agreement with the 1-D analytical solution. For a 3-D model containing a time-lapse resistivity anomaly representing a reservoir, convergence is achieved within 20 iterations at all frequencies, and the iterative results exhibit consistent apparent-resistivity and phase responses with those of the PARDISO direct solver, while peak memory usage is reduced by approximately 60%. These results demonstrate that the proposed iterative solution scheme significantly reduces memory consumption without compromising accuracy, and is capable of supporting the computational demands of multi-vintage repeated forward modeling, thus providing an efficient and viable approach for large-scale 3D time-lapse MT forward modeling of reservoirs.