In this paper, we propose a wave-equation-based travel-time seismic tomography method with a detailed description of its step-by-step process. First, a linear relationship between the travel-time residual $ Δt = T^{obs}–T^{syn} $ and the relative velocity perturbation $ δc(x)/c(x) $ connected by a finite-frequency travel-time sensitivity kernel $ K(x) $ is theoretically derived using the adjoint method. To accurately calculate the travel-time residual $ Δt $, two automatic arrival-time picking techniques including the envelop energy ratio method and the combined ray and cross-correlation method are then developed to compute the arrival times Tsyn for synthetic seismograms. The arrival times $ T^{obs} $ of observed seismograms are usually determined by manual hand picking in real applications. Travel-time sensitivity kernel $ K(x) $ is constructed by convolving a~forward wavefield $ u(t,x) $ with an adjoint wavefield $ q(t,x) $. The calculations of synthetic seismograms and sensitivity kernels rely on forward modeling. To make it computationally feasible for tomographic problems involving a large number of seismic records, the forward problem is solved in the two-dimensional (2-D) vertical plane passing through the source and the receiver by a high-order central difference method. The final model is parameterized on 3-D regular grid (inversion) nodes with variable spacings, while model values on each 2-D forward modeling node are linearly interpolated by the values at its eight surrounding 3-D inversion grid nodes. Finally, the tomographic inverse problem is formulated as a regularized optimization problem, which can be iteratively solved by either the LSQR solver or a~nonlinear conjugate-gradient method. To provide some insights into future 3-D tomographic inversions, Fréchet kernels for different seismic phases are also demonstrated in this study.