Flow in a multiaquifer porous system can be simulated by the so-called "quasi three-dimensional" models. When heterogeneous and/or aquitards with nonlinear hydrogeologic behavior are considered, a fully numerical approach is required for the model solution. If the finite element method is used to integrate the partial differential flow equations, the final solution of large systems is required. In the present article, an original iterative solution strategy is developed. The global system is decoupled into a number of smaller subsystems consistent with the geologic structure (aquitards and aquifers) of the multiaquifer system. The aquifer and the aquitard equations are solved separately with the modified conjugate gradient and the Thomas algorithms, respectively, while the final coupled solution is obtained with an iterative procedure equivalent to a Block Jacobi scheme. The procedure can be efficiently implemented on a parallel super-computer distributing the computational load, so that two successive blocks (related to an aquifer and the underlying aquitard) are solved on the same processor. The procedure is analyzed with linear porous media, where the convergence is theoretically ensured. The results obtained with a realistic linear multiaquifer system, employing a massively parallel computer like the CRAY T3D, have pointed out the high degree of parallelization of the algorithm. Comparison with the parallel implementation of the Block SOR and Block Gauss-Seidel schemes shows that parallel Block Jacobi performs significantly better with a reduction of the elapsed times, which depends on the rate of leakage between neighboring aquifers.