Frequency-domain methods, which are typically applied to 3D magnetotelluric (MT) modeling, require solving a system of linear equations for every frequency of interest. This is memory and computationally intensive. We developed a finite-difference time-domain algorithm to perform 3D MT modeling in a marine environment in which Maxwell’s equations are solved in a so-called fictitious-wave domain. Boundary conditions are efficiently treated via convolutional perfectly matched layers, for which we evaluated optimized parameter values obtained by testing over a large number of models. In comparison to the typically applied frequency-domain methods, two advantages of the finite-difference time-domain method are (1) that it is an explicit, low-memory method that entirely avoids the solution of systems of linear equations and (2) that it allows the computation of the electromagnetic field unknowns at all frequencies of interest in a single simulation. We derive a design criterion for vertical node spacing in a nonuniform grid using dispersion analysis as a starting point. Modeling results obtained using our finite-difference time-domain algorithm are compared with results obtained using an integral equation method. The agreement was found to be very good. We also discuss a real data inversion example in which MT modeling was done with our algorithm.