The computational burden of frequency-domain full-waveform inversion (FWI) of wide-aperture fixed-spread data is conventionally reduced by limiting the inversion to a few discrete frequencies. In this framework, frequency-domain seismic modeling is performed efficiently for multiple sources by solving the linear system resulting from the discretization of the time-harmonic wave equation with the massively parallel sparse direct solver. Frequency-domain seismic modeling based on the sparse direct solver (DSFDM) requires specific design finite-difference stencils of compact support to minimize the computational cost of the lower-upper decomposition of the impedance matrix in terms of memory demand and floating-point operations. A straightforward adaptation of such finite-difference stencil, originally developed for the (isotropic) acoustic-wave equation, is proposed to introduce vertical transverse isotropy (VTI) in the modeling without any extra computational cost. The method relies on a fourth-order wave equation, which is decomposed as the sum of a second-order elliptic-wave equation plus an anellipticity correction term. The stiffness matrix of the elliptic-wave equation is easily built from the isotropic stiffness matrix by multiplying its coefficients with factors that depend on Thomsen’s parameters, whereas the anelliptic term is discretized with a parsimonious second-order staggered-grid stencil. Validation of DSFDM against finite-difference time-domain modeling performed in various synthetic models shows that a discretization rule of four grid points per minimum wavelength provides accurate DSFDM solutions. Moreover, comparison between real data from the Valhall field and DSFDM solutions computed in a smooth VTI subsurface model supports that the method can be used as a fast and accurate modeling engine to perform multiparameter VTI FWI of fixed-spread data in the viscoacoustic approximation.