The spatial variability of soil hydraulic properties has to be considered to provide realistic predictions of unsaturated flow and transport at field scale. Steady unsaturated flow through randomly heterogeneous soils is analyzed and upscaled in this study. Numerical experiments of unsaturated flow have been performed on single realizations of unsaturated soil parameters generated by using threedimensional turning band random generator. In the first part of the study a linear relation for unsaturated hydraulic conductivity as a function of soil suction (.) given by K = Ksea. is used, with constant beta and saturated hydraulic conductivity (K-s) as perfectly correlated random space functions. In the first part of the study, linear conductivity relation given by Gardner is used, with Gardner's constant (beta) and saturated conductivity (Ks) as perfectly correlated random space functions. The steady flow fields generated are analyzed to explore the effect of variability of the product (a): the pore-size distribution parameter whose reciprocal represents the characteristic capillary length (h(cap)) and the correlation length (lambda Z) of the generated soil properties. A set of one-dimensional simulations are performed, with various values of the product beta lambda(Z) ratio of correlation length to capillary dispersion length (lambda(z)/ hcap = 1, >> 1 and e 1) used to compare the moisture distribution tendency. The resulting steady pressure fields are compared and analyzed. In some cases, the soil behaves as homogeneous, and in other cases as stratified. The second part of the study seeks to identify simple upscaling laws for block-scale nonlinear constitutive relationships under high tension, with the parameters of their measurement scale counterparts represented by the van Genuchten model. Random fields are generated to represent layering of perfect and imperfectly stratified soils, and numerical simulations are performed in flows that are perpendicular and parallel to stratification. Upscaling of hydraulic conductivity curve Kd.and moisture retention curve theta(Psi) is performed for perfectly stratified and imperfectly stratified soils, and anisotropy is also studied for both types of layering considered. It is demonstrated that the upscaled nonlinear curves can be represented by simple exponential relations and/ or by linear relations obtained by piecewise linearization. The linear relations obtained can be expressed in equivalent terms that can be used in analytical solutions to directly obtain effective hydraulic conductivity values. (C) 2016 American Society of Civil Engineers.