Abstract: Abstract: The performance evaluation of regional scale border irrigation plays an important role in the improvement of surface irrigation management level, but the existing models present the shortcomings such as less accurately capturing the spatial variability of regional scale irrigation parameters. So the exiting models cannot be effectively applied to analyze the performance of regional scale border irrigation system. To solve this problem, the border surface-subsurface flow model was applied to describe the surface and subsurface water flow. The conservative complete hydrodynamic equation and the Richards equation were applied to describe the surface water flow and subsurface water flow in border irrigation, respectively. The finite-volume approach was applied to spatially discretize these governing equations to obtain good mass conservation ability. The Picard iteration approach was introduced to obtain the linearization of this nonlinear algebraic system. The mass conservation component of surface flow model and subsurface flow model were iteratively coupled at the same time step to obtain the convergence value of surface flow depth, and then the momentum conservation component of surface flow model was externally coupled based on the convergence value of both the surface flow depth and infiltration rate to update the surface flow velocity. Solutions were numerically computed using an improved hybrid numerical method for surface flow model and a proposed numerical solution method with high-order accuracy for subsurface flow model. Monte-Carlo sampling method was used to accurately capture the spatial variability of regional scale irrigation parameters and generate a large number of border irrigation parameters samples, which were input to border surface-subsurface model, respectively. Consequently, the border surface-subsurface water flow processes of regional scale could be accurately simulated. Three times border irrigation experiments at regional scale were performed to validate the proposed model in Mawan Irrigation District, located in Dongying City, Shandong Province. Soil samples were collected at 4 depths from the top, middle, and bottom of the border field to analyze soil bulk density, soil particle size distribution, and soil moisture. The soil hydraulic parameters were transformed from the abovementioned soil properties using the Rosetta model. The relative elevation values of the border bottom were observed using a water level gauge. The surface flow depth was measured using water depth measuring devices, which were placed at every observation point before the irrigation was initiated. The surface flow depth was used to estimate Manning's roughness coefficient. And unit discharge, border length, and border width were observed in March 2012, November 2012, and March 2013. The validation results showed that the proposed model could well simulate regional scale border irrigation processes, and presented very good modeling accuracy. Specifically, the relative errors between the measured and simulated values by the proposed stochastic parameter irrigation model were 9.95%-12.23% and 8.39%-10.21% for irrigation quota and field water utilization coefficient, respectively. By contrast, the relative errors of irrigation quota and field water utilization coefficient based on the existing deterministic parameter irrigation models were 14.15%-16.78% and 13.87%-15.88%, respectively. Additionally, the cumulative probability distribution trends of average soil moisture after irrigation between the measured and simulated values present the satisfactory performance. Thus, the proposed model overcomes the shortcomings of existing models and provides a useful numerical analysis tool for the management and design of regional scale border irrigation system.