The aim of this study was to investigate the formation, extents, and evolution of quasi-saturated zones during and after rainfall using numerical method. First, the TOUGH2/EOS3 simulator is modified to consider the hysteretic characteristic curves, which are verified by simulating a laboratory experiment conducted by other researchers and comparing the results with the measured data. Then, the validated model is used to investigate the features of quasi-saturated zones induced by rainfall infiltration of four soil types. As rainwater infiltrates soil, a quasi-saturated zone forms near the ground surface because infiltrating water drives pore air downwards; then, the pore-air pressure increases, the matric suction decreases, and the pore-water pressure becomes positive. Increased pore-air pressure in unsaturated zones provides excess energy to moisture above the groundwater table, allowing the moisture to overcome the matric suction and transfer the pore-water pressure when the pore-air pressure exceeds the matric suction. Then, a quasi-saturated zone forms above the groundwater table as well. The quasi-saturated zone near the surface becomes unsaturated when the rain stops, and airflow in the unsaturated zone recedes, which causes the quasi-saturated zone above the groundwater table to start dissipating. The quasi-saturated zone above the groundwater table reforms after the wetting front reaches the groundwater table. The formation of quasi-saturated zones is therefore induced by airflow pushing or water recharging, causing the pore-air pressure to exceed the matric suction; i.e., the pore-water pressure becomes positive. Differences between the quasi-saturated zones in four soil types are also explored. We conclude that the formation and dissipation of quasi-saturated zones are delayed in clay and silt because of delayed airflow.