In this research, the variability analysis methods are developed to evaluate the effects of variation in the soil constitutive properties on theseismic ground motions and hence on the response of soil-(embedded) structure systems. The covariance matrices for the stresses, strains, and displacements in the soil and structure regions are analyzed in detail according to first- and second-order approximations. These results are verified against Monte Carlo simulations using a super computer. The algebraic complexities of the formulate depicting the higher-order statistical approximations will be handled by symbolic manipulation languages germain to artificial intelligence procedures. Automated sequential refinements of the final results of the statistical correlations are explored. Further, a set of nonlinear constitutive equations for the soil which are amenable both to finite element analysis and stochastic equivalent linearization are developed.