This paper presents a full-waveform inversion method for reconstructing the temporal and spatial distribution of unknown, incoherent dynamic traction in a heterogeneous, bounded solid domain from sparse, surficial responses. This work considers SH wave motions in a two-dimensional ( D) domain. The partial-differential-equation (PDE)-constrained optimization framework is employed to search a set of control parameters, by which a misfit between measured responses at sensors on the top surface induced by targeted traction and their computed counterparts induced by estimated traction is minimized. To mitigate the solution multiplicity of the presented inverse problem, we employ the Tikhonov (TN) regularization on the estimated traction function. We present the mathematical modeling and numerical implementation of both optimize-then-discretize (OTD) and discretize-then-optimize (DTO) approaches. The finite element method (FEM) is employed to obtain the numerical solutions of state and adjoint problems. Newton's method is utilized for estimating an optimal step length in combination with the conjugate-gradient scheme, calculating a desired search direction, throughout a minimization process.Numerical results present that the complexity of a material profile in a domain increases the error between reconstructed traction and its target. Second, the OTD and DTO approaches lead to the same inversion result. Third, when the sampling rate of the measurement is equal to the timestep for discretizing estimated traction, the ratio of the size of measurement data to the number of the control parameters can be as small as : in the presented work. Fourth, it is acceptable to tackle the presented inverse modeling of dynamic traction without the TN regularization. Fifth, the inversion performance is more compromised when the noise of a larger level is added to the measurement data, and using the TN regularization does not improve the inversion performance when noise is added to the measurement. Sixth, our minimizer suffers from solution multiplicity less when it identifies dynamic traction of lower frequency content than that of higher frequency content. The wave responses in a computational domain, induced by targeted traction and its reconstructed one, are in excellent agreement with each other. Thus, if the presented dynamic-input inversion algorithm is extended in realistic D settings, it could reconstruct seismic input motions in a truncated domain and, then, replay the wave responses in a computational domain.