  There is growing interest in discovering interpretable, closed-form equations
for subgrid-scale (SGS) closures/parameterizations of complex processes in
Earth systems. Here, we apply a common equation-discovery technique with
expansive libraries to learn closures from filtered direct numerical
simulations of 2D turbulence and Rayleigh-B\'enard convection (RBC). Across
common filters (e.g., Gaussian, box), we robustly discover closures of the same
form for momentum and heat fluxes. These closures depend on nonlinear
combinations of gradients of filtered variables, with constants that are
independent of the fluid/flow properties and only depend on filter type/size.
We show that these closures are the nonlinear gradient model (NGM), which is
derivable analytically using Taylor-series. Indeed, we suggest that with common
(physics-free) equation-discovery algorithms, for many common systems/physics,
discovered closures are consistent with the leading term of the Taylor-series
(except when cutoff filters are used). Like previous studies, we find that
large-eddy simulations with NGM closures are unstable, despite significant
similarities between the true and NGM-predicted fluxes (correlations $> 0.95$).
We identify two shortcomings as reasons for these instabilities: in 2D, NGM
produces zero kinetic energy transfer between resolved and subgrid scales,
lacking both diffusion and backscattering. In RBC, potential energy
backscattering is poorly predicted. Moreover, we show that SGS fluxes diagnosed
from data, presumed the ''truth'' for discovery, depend on filtering procedures
and are not unique. Accordingly, to learn accurate, stable closures in future
work, we propose several ideas around using physics-informed libraries, loss
functions, and metrics. These findings are relevant to closure modeling of any
multi-scale system.
