An upscaling algorithm to simulate unsaturated flow in horizontally heterogeneous soils at field scale under flood irrigation was developed. Based on the hydraulic functions of van Genuchten's type and the variables of dimensionless time, spatial position, negative water pressure were applied, the Richards' flow equation was transformed into the dimensionless form by using these dimensionless variables. Because there was no any parameter strongly dependent on scale in the transformed dimensionless Richards' equation, the solved problem was in nature equivalent to upscaling original problem. Solving the transformed equation by Galerkin finite element method only one time, combined with the dimensionless simulating time of each soil column, the distribution of pressure head of each column representing the soil profile of the column can be obtained. For the parameter n, slightly dependent on scale in the transformed Richards' equation and the parameter á, strongly dependent on scale in initial and boundary condition, we used the power-averaging technique to upscale these parameters. Firstly, the calculating procedures of new method were calibrated by two numerical experiments. Soil moisture in two different soil textures could be reasonably well simulated. Secondly, to verify the simulating efficiency of the developed method, we adopted the equivalent á values corresponding to the averaging exponent p=0, ±1, and ±0.5 to carry out numerical calculation. Different combinations of á and n were applied for comparing the accuracy of numerical simulation and it showed that n had little effect on the simulated result. Compared with the fine-mesh solution of horizontal plane using HYDRUS-1D software, it indicated that for example 1, whose initial profile of water negative pressure was in equilibrium, when p=0 and -0.5, the simulated results were well represented, in the wetting front dramatically fluctuation depths 10 cm and 20 cm, the simulated pressure heads were of small relative errors: -4.06% (p=0) and -10.41% (p=-0.5), 7.24% (p=0) and 2.53% (p=-0.5) respectively; for example 2, which had a constant initial profile, when p=-0.5, in general, the simulated results were best, in the wetting front dramatically fluctuation depths 50 cm and 60 cm, the simulated pressure heads were of small relative errors: -5.27% and 3.58% respectively. Obviously, the proposed method ensured better simulated results and provided a new approach to solve water movement in horizontally heterogeneous unsaturated soils.