Marine vibrators are increasingly being recognized as a viable alternative to seismic airguns for ocean-bottom acquisition due to their ability to generate more low-frequency content and their less adverse impact on marine wildlife. However, their use introduces processing challenges, such as the Doppler effect and time-dependent source-receiver offsets, which are negligible in conventional airgun acquisition. In addition, the time-varying nature of the sea surface during the multi-second acquisition time introduces further challenges for processing and inversion. To accurately account for source motion and time-varying sea surface effects in seismic data processing, we develop a reliable and robust numerical modeling tool. We use a mimetic finite-difference method in a generalized coordinate system to model the full acoustic wavefield triggered by a moving source in the presence of a time-varying sea surface. Our approach uses a time- and space-dependent coordinate transformation, which tracks the source movement and conforms to the irregular time-varying sea surface, to map an irregular physical domain in Cartesian coordinates to a regular computational domain in generalized coordinates. We formulate this coordinate transformation such that both coordinate systems conformally match below the ocean-bottom level. Numerical examples demonstrate that this approach is accurate and stable, even for an unrealistically exaggerated sea state. This computational tool is not limited to modeling, but could also be used to develop advanced processing techniques for marine vibrator data, such as imaging and inversion.