Abstract: Drip irrigation has been recognized as an efficient method to improve water and nitrogen use efficiency. Nitrate leaching is becoming an important consideration for the design, operation, and management of drip irrigation systems. Recent studies have indicated that soil properties such as soil texture and saturated hydraulic conductivity could impose a significant effect on nitrate leaching. However, most of these studies were conducted on an assumption of homogeneous soil. In the current study, the effect of the spatial variability of soil properties on nitrate leaching was evaluated by a model that was established and solved numerically by using the HYDRUS-2D package. In the simulation, the statistical distribution of the saturated hydraulic conductivity and the initial soil water content in the field were assumed to be able to be represented by a logarithmic normal distribution and a normal distribution function, respectively. A moderate variability for both saturated hydraulic conductivity and initial soil water content were considered with the coefficient of variation (CV) of 0.40 and 0.15, respectively. One hundred samples of the spatial varied saturated hydraulic conductivity and initial soil water content at a given CV value were generated randomly by the Monte Carlo method. The nitrate leaching ratios for these randomly generated parameters of saturated and initial soil water content were determined by the numerical model. After the positive correlation between the nitrate leaching ratio and the saturated hydraulic conductivity and the initial soil water content was confirmed, a joint distribution function of saturated hydraulic conductivity, initial soil water content, and nitrate leaching ratio was constructed using a Gumbel-Hougaard trivariate copula function. The maximum likelihood approach was used to determine the parameters for the Gumbel-Hougaard trivariate copula function and the Kolmogorov-Smirnov test indicated that the empirical distribution of the nitrate leaching ratios could be well represented by the trivariate copula function with a root mean square error (RMSE) value of 0.046. To quantify the risk of nitrate leaching, the conditional probability distribution function of the nitrate leaching ratio under given saturated hydraulic conductivities and initial soil water contents was estimated from the derived joint distribution function. The results indicated that the probability density function of the nitrate leaching ratio could be represented by an exponential function and the variability of saturated hydraulic conductivity and initial soil water content imposed an important effect on the nitrate leaching. For the sandy loam soil having medium variabilities of saturated hydraulic conductivity and initial soil water content, the simulated nitrate leaching ratio reached 13.2%, being more than double of the value for the homogeneous soil (6.4%). The probability under the condition that the nitrate leaching ratio exceeds a given threshold (6.4%) generally increases with an increasing saturated hydraulic conductivity and initial soil water content. It might be concluded that constructing a joint distribution function between the nitrate leaching and soil properties in the field (e.g., saturated hydraulic conductivity and initial soil water content) is a promising method to assess the risk of nitrate leaching while the spatial variability needs to be considered.