Accurate evaluation of soil dynamic properties is essential for seismic response analyses of sites. In a number of studies, site properties have been identified using one-dimensional analyses. Such analyses uncouple the two-dimensional (horizontal) response of soil deposits, which is inherently coupled. This paper presents a system identification technique that takes into account the coupled two-directional response of soil deposits. The technique employs non-parametric estimates of the shear stresses derived from acceleration records provided by a vertical (downhole) array. A multi-yield surface plasticity approach is used to model the multi-dimensional stress-strain relation. The identification technique is first verified using finite elements computational simulations. This technique was then used to assess the coupled response of the Wildlife liquefaction research site (Imperial Valley, California). The identified shear moduli and shear wave velocities were found to be in a very good agreement with those measured in the field using crosshole seismic testing.