High-resolution simulations that fully resolve all spatiotemporal scales of geophysical and turbulent flows remain a challenge in large ocean domains. Large-eddy simulations (LES) make these computations tractable by filtering out subgrid-scale (SGS) features, but require accurate closures to remain stable and faithful; without them, solutions can drift, lose energy at the wrong rate, or develop spurious coastal artifacts. Classical analytical closures based on the eddy-viscosity hypothesis, such as the Smagorinsky and Leith models and their dynamic variants, were developed primarily for three-dimensional homogeneous turbulence. A recent benchmarking study has shown that they logically only weakly capture the SGS forcing in two-dimensional vorticity flows in the presence of coastal boundaries and interior landforms, motivating the development of data-driven closures. Among such approaches, neural-network closures have shown promise but typically return only a deterministic point estimate of the SGS term, while the mapping from resolved to unresolved scales is fundamentally non-invertible and the closure is therefore intrinsically stochastic. This non-uniqueness becomes especially pronounced in nonstationary flows, where the wake statistics themselves drift in time and a single deterministic correction can likely not represent the spread of admissible SGS responses.

In this work, we develop and evaluate sparse and deep Gaussian process (GP) closures for under-resolved, non-stationary two-dimensional classical and β-plane vorticity flows past idealized obstacles, where non-stationarity is driven by a time modulated inflow velocity U_∞(t) that produces a continuously evolving wake, with shedding frequency, wake width, and subgrid-scale statistics all drifting along the trajectory.