Modeling of seismic-wave propagation in anisotropic medium with irregular topography is beneficial to interpret seismic data acquired by active and passive source seismology conducted in areas of interest such as mountain ranges and basins. The major challenge in this context is the difficulty in tackling the irregular free-surface boundary condition in a Cartesian coordinate system. To implement surface topography, we use the boundary-conforming grid and map a rectangular grid onto a curved grid. We use a stable and explicit second-order accurate finite-difference scheme to discretize the elastic wave equations (in a curvilinear coordinate system) in a 3D heterogeneous transversely isotropic medium. The free-surface boundary conditions are accurately applied by introducing a discretization that uses boundary-modified difference operators for the mixed derivatives in the governing equations. The accuracy of the proposed method is checked by comparing the numerical results obtained by the trial algorithm with the analytical solutions of the Lamb’s problem, for an isotropic medium and a transversely isotropic medium with a vertical symmetry axis, respectively. Efficiency tests performed by different numerical experiments illustrate clearly the influence of an irregular (nonflat) free surface on seismic-wave propagation.