Flow parameters of environmental processes within subsurface media are represented by means of random fields (RFs) [l] due to the natural heterogeneity and the incomplete characterization of the media. In the case of saturated single-phase flow the main model parameter is the hydraulic conductivity RF, the statistical properties of which are inferred from experimental data. Samples (realizations) of the hydraulic conductivity RF that honor the statistics can be generated on a nwnerical grid by means of various simulation algorithms [2]. Then, the partial differential equation of flow can be solved using standard numerical techniques [3]. This approach is inefficient for obtaining stochastic averages of flow properties, since a large number of realizations is usually required [4]. It is well-known that effective parameter. which involve stochastic averages over small-scale fluctuations, can provide efficient estimators of large scale behavior [5]. However, the evaluation of the non-linear stochastic averages involved in effective hydraulic conductivity (EHC) estimation is very complicated. Until recently, the EHC analysis has mainly been limited to asymptotic estimates valid at the limit of vanishing correlation lengths, and low-order perturbation methods that are inaccurate for large heterogeneity. More detailed models involve non-local average constitutive relations and effective kernel functions instead of uniform effective parameters, e.g. [6]. In addition, experimental measurements in the field exhibit a strong scaling behavior of the EHC [7]. Scaling has been examined within certain general theories of measurement, e.g. [8.9], but a flexible and practical model for the prediction of the scaling behavior is still missing. In this paper we investigate a non-perturbative replica-variational approach (IO] that leads to EHC kernel functions through a coupled system of integral equations. This approach can also capture scaling effects due to the cut-off of the stochastic fluctuations by the support scale.