To accurately simulate seismic wave propagation for the purpose of developing modern land data processing tools, especially full-waveform inversion (FWI), we have developed an efficient high-order finite-difference forward-modeling algorithm with the capability of handling arbitrarily shaped free-surface topography. Unlike most existing forward-modeling algorithms using curvilinear grids to fit irregular surface topography, this finite-difference algorithm, based on an improved immersed boundary method, uses a regular Cartesian grid system without suffering from staircasing error, which is inevitable in a conventional finite-difference method. In this improved immersed boundary finite-difference (IBFD) algorithm, arbitrarily curved surface topography is accounted for by imposing the free-surface boundary conditions at exact boundary locations instead of using body-conforming grids or refined grids near the boundaries, thus greatly reducing the complexity of its preprocessing procedures and the computational cost. Furthermore, local continuity, large curvatures, and subgrid curvatures are represented precisely through the employment of the so-called dual-coordinate system — a local cylindrical and a global Cartesian coordinate. To properly describe the wave behaviors near complex free-surface boundaries (e.g., overhanging structures and thin plates, or other fine geometry features), the wavefields in a ghost zone required for the boundary condition enforcement are reconstructed accurately by introducing a special recursive interpolation technique into the algorithm, which substantially simplifies the boundary treatment procedures and further improves the numerical performance of the algorithm, as demonstrated by the numerical experiments. Numerical examples revealed the performance of the IBFD method in comparison with a conventional finite-difference method.