Wave-equation reflection waveform inversion (RWI) is a promising method to reconstruct the background velocity model with reflection data. But it is difficult to precondition this highly nonlinear inverse problem for efficient convergence and reliable model updating. In the context of full-waveform inversion (FWI), the second-order derivative of the objective function with respect to the model parameters (i.e., the Hessian operator) can be used with the second-order optimization to overcome similar difficulties. To our knowledge, however, there are few such studies for RWI. This motivates us to derive and investigate the Hessian operator using the Born approximation in the context of RWI. For the same experimental setup with a toy model, we observe that the Hessian matrix of RWI exhibits a more compact structure than its counterpart in the context of FWI because they are constructed with the sensitivity kernels of very different spatial distribution features. Accordingly, we have developed a truncated Gauss-Newton RWI method consisting of two nested loops: An outer loop updates the background model in the manner of nonlinear waveform inversion, and an inner loop aims to iteratively estimate the background model update in a matrix-free manner by solving the Gauss-Newton equation. The workflow consists of alternating estimation of the perturbation and the background velocity model components, of which the former is done with a linear least-squares reverse time migration using the conjugate-gradient method. A synthetic example on the overthrust model and a field data example demonstrate that the Hessian-based RWI can significantly improve the estimation of intermediate-wavenumber components and promote velocity model building and seismic imaging with precritical reflection data.