Localized narrowband low-frequency shear waves can be non-invasively generated within tissue, by a modulated finite-amplitude radiation force, resulting from the interference of two focused quasi-CW ultrasound beams of slightly different frequencies. Assuming a Voigt viscoelastic model, this paper describes the use of a finite-element-method model, to simulate two-dimensional shear-wave propagation in viscoelastic media, containing circular inclusions (lesions). Using this model, an inverse approach is used to extract maps of the local shear modulus and viscosity. The performance is evaluated based on three metrics: the lesion contrast, the contrast-transfer-efficiency (CTE), and the contrast-to-noise ratio (CNR). Modified definitions of these metrics are proposed and used in order to account for the time-varying nature of the shear waves and the inverse reconstruction algorithm. In the absence of any noise, it is shown that accurate reconstruction can be achieved not only with the fundamental, but also with the higher harmonics, as well as, with a low-frequency component that occurs for high viscosity and high modulation frequencies. For low-viscosity conditions, the lesion contrast, CTE, and CNR are shown to exhibit very good performance not only for the fundamental, but also, for the higher harmonics. In the case of increased viscosities and modulation frequencies, the generated low-frequency component is shown to provide superior contrast performance even when compared to that of the fundamental. The effects of noise on the reconstruction quality are examined. Depending on the lesion and background properties, it is shown that noise can seriously degrade reconstruction from the higher harmonics.