Abstract: Non-point source phosphorus pollution generated from irrigated farmlands is one of the main causes of local and regional eutrophication. However, current phosphorus pollution models either do not include the water movement in irrigation and drainage process or do not consider the phosphorus transformation under the exchanging aerobic and anaerobic conditions. Therefore, we developed a physically based phosphorus pollution model to quantitatively describe the water movement and phosphorus fate and transport processes in irrigated paddy fields in plain areas. The simulation of the runoff yield in an irrigation area was based on the water balance equations describing the water input and output of the paddy fields and the motion wave equations describing the water movement in the drainage channel networks. The simulation of the excess phosphorus yield was based on the convection diffusion equations and a phosphorus transformation model considering the soil sub-stratification-the cultivated horizon was sub-divided into aerobic and anaerobic layers. In this way, the changes in dissolved oxygen and the processes of phosphorus transformation in different soil layers caused by the alternating wet and dry conditions could be quantitatively described in details. The phosphorus flux of diffusion, particle mixing and infiltration between the water layer, the aerobic soil layer and the anaerobic soil layer were also quantified. The model was calibrated and verified with the observed ponding water depth, drainage discharge, and phosphorus concentrations in the runoff and soil water in one experimental paddy field and two typical drainage ditches in Heping Irrigation District, Heilongjiang, China in 2018. The simulated drainage discharge and phosphorus concentrations of the experimental paddy field and the drain ditches agreed well with the observations. The Nash-Sutcliffe efficiency coefficient (NSE) and coefficient of determination (R2) of the simulated drainage discharge were greater than 0.820 and 0.815, respectively. And the NSE and R2 of simulated total phosphorus concentration were greater than 0.811 and 0.821, respectively. The simulated vertical distribution of the soil soluble phosphorus obtained by considering the aerobic and anaerobic layers of the cultivated horizon were closer to the in situ observation than the results obtained with the same model but do not consider the soil sub-stratification. Then, the verified phosphorus pollution model was used to estimate the non-point source phosphorus pollution in the whole Heping Irrigation District. The phosphorus loss through drainage and leakage during the growth stages of rice was 1.88 kg/hm2, which was about 5.7% of the phosphorus input from fertilization and irrigation. Among the 1.88 kg/hm2 phosphorus loss, the phosphorus output load of runoff at the tillering stage (0.85 kg/hm2), and the jointing and booting stage (0.60 kg/hm2) was the first and second largest loss, due to rainfall washout of soil phosphorus. The loss by leakage output load was the second and first largest at the soaking stage (0.11 kg/hm2) and the tillering stage (0.16 kg/hm2), due to the basic fertilizer and the early booting stage fertilizer. For the whole Heping Irrigation District, the total excess phosphorus exported from the first ditch (1.40 t) and the fourth ditch (1.39 t) were the first and second largest, due to their larger control area of the irrigation district. Overall, the physically based phosphorus pollution model developed in this study included the water movement in irrigation and drainage process, considered the phosphorus transformation under the exchanging aerobic and anaerobic conditions caused by the alternating wet and dry conditions, and provided more accurate estimation of phosphorus fate and transport in irrigated paddy fields in plain areas.