Numerical simulation models of water flow in saturated-unsaturated soil systems are important tools for quantitatively understanding hydrological processes. Boundary conditions required in the simulation models play important roles for characterizing flow processes and hydrological fluxes. Whether a numerical scheme is feasible and reliable to describe soil-water problems lies in an appropriate procedure for the boundary conditions (van Dam & Feddes, 2000). Notably, for vadose zone characterization, many methods categorized as upscaling and downscaling techniques (see reviews by Vereecken et al. [2007] and Anderson et al. [2003]), forward and inverse modeling approaches (see reviews by Hopmans and Šimůnek [1999] and Vrugt et al. [2008]), practically depend on the nature of boundary conditions and are affected by the uncertainty associated with the boundary conditions, and thus their applicability and validation require appropriate determination of boundary conditions (Vereecken et al., 2008). This study focuses specially on soil surface flux-type boundary conditions, which are prevailing at moderate weather and soil wetness conditions (van Dam & Feddes, 2000) and depend on soil-plant-atmosphere interaction and are determined by several boundary fluxes, including evaporation, precipitation, irrigation, drainage, etc. Traditionally, these fluxes are derived by experimental measurements, analytical or numerical models based on water and energy balance. However, due to the variabilities in rainfall, wind, and moisture flux across a wide range of scales, the variability in agricultural management practices within and between individual fields, and the variability in the ecology of the land surface (Vereecken et al., 2007), these fluxes often vary spatiotemporally, making their measurements and determinations difficult, expensive and time