The main purpose of this paper is to reduce the dimensionality of unknown coefficient vectors in finite element (FE) solutions of two-grid FE Crank-Nicolson (CN) (TGFECN) format for the unsaturated soil water flow (USWF) problem with two strong nonlinear terms by using proper orthogonal decomposition (POD). For this purpose, our first step involves designing a time semi-discrete CN (TSDCN) scheme for the USWF problem and demonstrate the existence, boundedness, and error estimations of TSDCN solutions. Thereafter, we discretize the TSDCN scheme using the two-grid FE method to create a new TGFECN format and prove the existence, boundedness, and error estimations of TGFECN solutions. The primary focus should be on reducing the dimension of unknown coefficient vectors of TGFECN solutions through the utilization of the POD technique in creating a novel format, referred to as dimension reduction iterative TGFECN (DRITGFECN) format, while establishing the existence, boundedness, and error estimations for DRITGFECN solutions. Lastly, we use two sets of numerical tests to exhibit the advantage of the DRITGFECN format. Due to the presence of two highly nonlinear terms in the unsaturated soil flow problem, the development and analysis of DRITGFECN format pose greater challenges and necessitates more advanced technical skills compared to previous studies. However, the significance and broad applications of this research make it a valuable subject for investigation.