Truncated fluid and ocean models omit subgrid physics and introduce numerical biases that degrade forecasts. We present a Bayesian, data-driven closure for 2-D finite-volume solvers that learns the dynamical discrepancy between low-resolution (LR) and high-resolution (HR) simulations. Using sparse variational Gaussian processes (GPs) and deep GPs, we map resolved features (local velocities and gradients) to a closure source term that corrects LR tendencies toward HR dynamics while quantifying predictive uncertainty. GPs can be well-suited to closure modeling in fluids because they encode smoothness/invariance via kernels, learn nonparametric mappings from data, and return uncertainty estimates alongside the mean correction. The trained GP is embedded intrusively into a numerical finite volume framework and evaluated online each coarse time step, keeping the closure consistent with the numerics.

We assess the approach on three test beds: (i) flow past a cylinder across multiple Reynolds numbers; (ii) tidally modulated flow past a cylinder with time-varying Reynolds number; and (iii) bottom gravity currents. Models are trained on HR downsamplings–LR pairs and then tested across different regimes. We evaluate performance by using field-wise errors and wake metrics: mean velocity profiles in the near and far wake, lift C_L and drag C_D coefficients, and Strouhal number St. Relative to LR baselines without closure, GP closures reduce L2 / L∞ errors of the resolved fields and bring mean velocity, C_D/C_L, and St closer to HR references across trained Reynolds numbers. The online GP closure adds negligible wall-clock cost relative to the fluid step, preserves the conservative finite-volume structure, and provides uncertainty estimates. Overall, these results demonstrate a practical, uncertainty-aware GP closure that improves coarse-grid fidelity for 2-D fluid and ocean flows, which could potentially be extended to 3-D ocean frameworks.