The parabolic partial differential equation describing fluid flow in partially saturated porous media, Richards' equation, is highly nonlinear due to pressure head dependencies in the specific soil moisture capacity and relative hydraulic conductivity terms. In order to solve Richards' equation several numerical techniques have been developed, which, starting from the discretization of the partial differential equation, produced even more accurate models, leading to complex and computationally expensive simulations for large-scale systems. A three-dimensional unsaturated flow modeling developed through a simulation environment based on cellular automata (CA) is described in this paper. The proposed model represents an extension of the original computational paradigm of cellular automata, because it uses a macroscopic CA approach where local laws with a clear physical meaning govern interactions among automata. This CA structure, aimed at simulating a large-scale system, is based on functionalities capable of increasing its computational capacity, both in terms of working environment and in terms of the optimal number of processors available for parallel computing. The model has been validated with reference multidimensional solutions taken from benchmarks in literature, showing a good agreement even in the cases where nonlinearity is very marked. Furthermore, some analyses have been carried out considering quantization techniques aimed at transforming the CA model into an asynchronous structure. The use of these techniques in a three-dimensional benchmark allowed a considerable reduction in the number of local interactions among adjacent automata without changing the efficiency of the model, especially when simulations are characterized by scarce mass exchanges. Finally, from a computational point of view the higher efficiency values were achieved running the model on a parallel architecture, obtaining a high speedup very close to the optimal with the maximum number of processors available.