We present deterministic ground motion simulations that account for the cyclic multiaxial response of sediments in the shallow crust. We use the Garner Valley in Southern California as a test case. The multiaxial constitutive model is based on the bounding surface plasticity theory in terms of total stress and is implemented in a high-performance computing finite-element parallel code. A major advantage of this model is the small number of free parameters that need to be calibrated given a shear modulus reduction curve and the ultimate soil strength. This, in turn, makes the model suitable for regional-scale simulations, where geotechnical data in the shallow crust are scarce. In this paper, we first describe a series of numerical experiments designed to verify the model implementation. This is followed by a series of idealized large-scale simulations in a 35 × 26 × 4.5 km 3 domain that encompasses the Garner Valley downhole array site, which is an instrumented and well-characterized site in Southern California. Material properties were extracted from the Southern California Earthquake Center Community velocity model, CVM-S4.26, considering its optional geotechnical layer, while the modulus reduction curves and soil strength were selected empirically to constrain the nonlinear soil model parameters. Our nonlinear simulations suggest that peak ground displacements within the valley increase relative to the linear case, while peak ground accelerations can increase or decrease, depending on the frequency content of the excitation. The comparisons of our simulations against hybrid three-dimensional-one-dimensional site response analyses suggest the inadequacy of the latter to capture the complexity of fully three-dimensional simulations.