Oceans are known to have complex and nuanced effects on seismic wave propagation; changing the behaviour of phases which propagate in the solid Earth, supporting the propagation of new classes of hydroacoustic waves, and acting as the source of the dominant ambient noise, the ocean microseism. Simulations with realistic ocean layers have historically been challenging, due to the need to explicitly mesh the low-sound-speed ocean and the challenge in implementing bathymetry with a conformal mesh. This has led to a number of simplifications becoming commonplace, including modelling the ocean as a ‘weight load’, and restriction to purely axisymmetric or local-scale domains. Such approximations are not appropriate in a high-frequency global seismology context, motivating us to develop a code capable of supporting realistic ocean layers with accurate bathymetry. Here, we present such an implementation of AxiSEM3D which should provide useful modelling output for comparison with new and exciting datasets from OBSs and Mermaids. We use this method to perform a high-frequency evaluation of the oft-used ‘ocean loading’ formulation, and find that it breaks down well above the 5-10 seconds period routinely suggested. Bathymetry is found to have a particularly substantial impact on surface-reflected phases (e.g. PP), and on the generation of ocean microseismic noise. We also consider the use of this method for interpreting observed water depth phases, and intend to work toward direct comparison to data.