Soft solids, such as gelatin or soft tissue, have a shear wave speed that is smaller than the compressional wave speed. Recent observations of shear shock wave generation in the brain, which can easily reach a large Mach number regime, suggest that the cubic nonlinear behavior of soft solids could be responsible for traumatic brain injury. However, currently there are no two-dimensional (2D) models describing the propagation of linearly-polarized shear shock waves in relaxing soft solids. These models are required to model the fundamental wave propagation physics like diffraction, focusing etc., and are related to traumatic brain injuries as it can be used to model the skull/brain morphology. In this work, we present a two-dimensional system of first-order equations modeling the propagation of shear shock waves in a relaxing soft solid, and then this system is solved numerically using the piecewise parabolic method, a high-order finite volume method. The numerical solutions are validated in two parts. First, the nonlinear component, which is designed for large Mach numbers, is validated with a step-shock Riemann problem. Then relaxation mechanisms based on a generalized Maxwell body, which model the non-classical attenuation that occurs in soft tissue, are compared to analytical solutions with an error of 5%-10%. The validation of attenuation also includes dispersion due to causality which is determined by the Kramers-Kronig relations. Finally, the full numerical method, which includes nonlinearity, attenuation, and dispersion, is compared to ultrasonic measurements of shear shock wave propagation in tissue-mimicking phantoms. Two experiments were performed based on high frame-rate ultrasound imaging and tracking in a gelatin phantom 1) planar shear shock waves, and 2) focused shear shock wave propagation. The experimental and numerical waveforms closely match, e.g. the RMS amplitude error is between 12.05% and 12.27%. Moreover, the frequency-content of the temporal signal was compared for third and fifth multiples of fundamental harmonic validating the generation of odd-harmonics due to cubic nonlinearity. Furthermore, the numerical tool was able to estimate the nonlinear parameter in the phantom (beta = 4.4) using a grid-parameter-search. In context of traumatic brain injury, the current method can be used to study the shear shock formation in 2D-sections of human skull, and can also be used for nonlinear parameter estimation in brain. (C) 2019 Elsevier Inc. All rights reserved.