Accurately modeling full-wavefield solutions at and near the seafloor is challenging for conventional single-domain elastic finite-difference (FD) methods. Because they treat the fluid layer as a solid with zero shear-wave velocity, the energy partitioning for body and surface waves at the seafloor is distorted. This results in incorrect fluid/solid boundary conditions, which has significant implications for imaging and inversion applications that use amplitude information for model building. To address these issues, here we use mimetic FD (MFD) operators to develop and test a numerical approach for accurately implementing the boundary conditions at a fluid/solid interface. Instead of employing a single “global” model domain, we partition the full grid into two subdomains that represent the acoustic and elastic (possibly anisotropic) media. A novel split-node approach based on one-sided MFD operators is introduced to distribute grid points at the fluid/solid interface and satisfy the wave equation and the boundary conditions. Numerical examples demonstrate that such MFD operators achieve stable implementation of the boundary conditions with the same (fourth) order of spatial accuracy as that inside the split-domain interiors. We compare the wavefields produced by the MFD scheme with those from a more computationally expensive spectral-element method to validate our algorithm. The modeling results help analyze the events associated with the fluid/solid (seafloor) interface and provide valuable insights into the horizontal displacement or velocity components (e.g., recorded in ocean-bottom-node data sets). The developed MFD approach can be efficiently used in elastic anisotropic imaging and inversion applications involving ocean-bottom seismic data.