Accurate reduced models of turbulence are desirable to facilitate the optimisation of magnetic-confinement fusion reactor designs. As a first step towards higher-dimensional turbulence applications, we use reservoir computing, a machine-learning (ML) architecture, to develop a closure model for a limiting case of electrostatic gyrokinetics. We implement a pseudo-spectral Eulerian code to solve the one-dimensional Vlasov–Poisson system on a basis of Fourier modes in configuration space and Hermite polynomials in velocity space. When cast onto the Hermite basis, the Vlasov equation becomes an infinitely coupled hierarchy of fluid moments, presenting a closure problem. We exploit the locality of interactions in the Hermite representation to introduce an ML closure model of the small-scale dynamics in velocity space. In the linear limit, when the kinetic Fourier–Hermite solver is augmented with the reservoir closure, the closure permits a reduction of the velocity resolution, with a relative error within 2 % for the Hermite moment where the reservoir closes the hierarchy. In the strongly nonlinear regime, the ML closure model more accurately resolves the low-order Fourier and Hermite spectra when compared with a naive closure by truncation and reduces the required velocity resolution by a factor of 16. One of the most important ongoing efforts in the pursuit of magnetic fusion energy is the study of microturbulence (Yoshida et al. 2025). While magnetic-confinement fusion plasmas are designed to be stable to macroscopic magnetohydrodynamic perturbations, the steep gradients in temperature and density that are present in all current designs also drive small-scale instabilities. These instabilities drive sustained turbulence, which then dominates the transport of particles, energy and momentum in the plasma. As the fusion industry turns its focus toward designing first-of-a-kind pilot plants, optimising these designs to mitigate turbulent heat losses is an urgent objective. Significant progress in turbulence optimisation for tokamaks (Highcock et al. 2018) and stellarators (Roberg-Clark et al. 2023; Kim et al. 2024) has been achieved. However, the space of candidate magnetic equilibria remains vast, and high-resolution turbulence simulations are computationally expensive. A unifying feature of fluid (Kraichnan 1959; Kolmogorov 1991), magnetohydrodynamic (Goldreich & Sridhar 1995) and gyrokinetic (Bañón Navarro et al. 2011) turbulence theory is that free energy tends to cascade from large scales to small scales through local interactions in phase space, and these cascades are well described by power laws. Therefore, models that can accurately calculate large-scale quantities, including spectra and heat fluxes, without needing to resolve small scales, are both desirable and plausible. In the field of fluid turbulence modelling, two successful approaches to constructing reduced, lower-resolution models that capture large-scale behaviour are Reynolds-averaged Navier–Stokes (Reynolds 1895) and large-eddy simulation (Deardorff 1970). Large-eddy simulation also has been used to build reduced models of gyrokinetic turbulence (Bañón Navarro et al. 2014). These formalisms accelerate turbulence simulations by averaging over time and filtering out the small-scale dynamics. Closure models provide another pathway to filter out the small scales while retaining accurate large-scale statistics. In kinetic plasma physics, Hammett & Perkins (1990) developed a landmark analytic linear closure for the low-frequency response to small Langmuir oscillations. An alternative method for building reduced models is to use data-driven, statistical techniques like machine learning (ML) and Bayesian optimisation. In recent years, the proliferation of ML models has led to advancements in pattern identification and model construction for magnetic fusion. Some example applications include improved prediction and avoidance of tokamak disruptions using the DECAF framework (Piccione et al. 2020) and on the DIII-D tokamak (Rea et al. 2019; Fu et al. 2020), improved active-feedback plasma control systems with reinforcement learning (Degrave et al. 2022; Seo et al. 2024; Kerboua-Benlarbi et al. 2024), optimised execution of gyrokinetic simulations with the PORTALS framework (Rodriguez-Fernandez, Howard & Candy 2022) and accelerated, low-resolution turbulence simulations (Greif, Jenko & Thuerey 2023; Castagna et al. 2024; Clavier et al. 2025). Data-driven, configuration space (Ma et al. 2020; Qin et al. 2023) and heat flux (Ingelsten et al. 2024; Huang, Dong & Wang 2025) closure models of Landau damping have also been developed. In this paper, we choose to use the one-dimensional, collisionless Vlasov–Poisson system as a test problem for implementing an ML closure in velocity space. This system is the one-dimensional, electrostatic limit of the full five-dimensional gyrokinetic system (Parker & Dellar 2015). It models dynamics parallel to the magnetic field, and it contains a quadratic nonlinearity present in the full gyrokinetic system. These features make the Vlasov–Poisson system a suitable foundation to assess the future potential for a closure model to be developed for gyrokinetic turbulence. In parallel velocity space, we use properties of the Hermite basis to cast the objective of reducing resolution requirements as a closure problem. The Hermite basis has been used extensively to solve problems in kinetic plasma dynamics (Armstrong 1967; Grant & Feix 1967; Joyce, Knorr & Meier 1971; Gibelli & Shizgal 2006; Zocco & Schekochihin 2011; Loureiro, Schekochihin & Zocco 2013; Kanekar et al. 2015; Parker & Dellar 2015; Jorge, Ricci & Loureiro 2017; Adkins & Schekochihin 2018; Mandell, Dorland & Landreman 2018; Frei et al. 2023). In this work, we opt to focus on velocity space because the Hermite basis naturally expresses the interaction between scales in velocity space as local interactions in a coupled hierarchy of moments. This property of the moment formalism directly presents a single-term target for a closure. Specifically, the interaction between resolved scales and unresolved scales in velocity space appears as a single unknown term in an equation in this representation. In contrast, the spectral representation of configuration space introduces some non-local interactions across scales, complicating the development of a closure model. Additionally, the Hermite basis explicitly expresses conservation laws for particle number, momentum and energy as evolution equations for the lowest-order moments. These properties allow us to develop a closure that preserves the low-order conservation laws. More details on the closure problem are presented in § 3.5. We use reservoir computing, an ML architecture, to implement a spectral, proof-of-concept, velocity-space closure for the moment hierarchy that maintains accurate spectra and reduces resolution requirements. The reservoir computing model maintains an internal representation of the system state, which it evolves through a directed graph with fixed weights. It then predicts the system state at time t+1 from the state at time t using a set of output weights trained by linear regression. Since its concurrent development by Jaeger (2001) and Maass, Natschläger & Makram (2002), reservoir computing has successfully been applied to the task of forecasting the dynamics of low-dimensional chaotic systems (Lu et al. 2017; Pathak et al. 2018). This ML paradigm is well suited to the task of modelling time-series data of physical processes. In climate physics, Arcomano et al. (2020) built a surrogate global atmospheric forecast model with reservoir computing, and in plasma physics, Jalalvand et al. (2022) used reservoir computing to analyse and classify Alfvén eigenmodes using DIII-D diagnostics data. In comparison with several other recurrent neural network architectures, reservoir computing has shown comparable ability to capture the statistics of chaotic systems, while requiring significantly less training time (Vlachas et al. 2020). These results motivate our decision to use reservoir computing for the closure model. Some ML methods, including reservoir computing, are more challenging to interpret than traditional physics models. Artificial neural network architectures are often ‘black boxes’ that generate layered, analytically intractable networks of nonlinear activation functions. Gaining insightful intuition from these networks can be daunting. We mitigate this problem by constraining the ML portion of our simulations to represent only the unresolved small scales, preserving analytic evolution equations for large scales. Our methods are analogous to a large-eddy simulation approach for phase space. The paper has the following structure. In § 2, we present our test problem, the normalised one-dimensional Vlasov–Poisson system. We then derive the projection of the system onto the pseudo-spectral Fourier–Hermite basis in § 3. We present the closure problem in § 3.5 and introduce our ML closure model in § 4 as a solution to that problem. In § 5, we present our ML model’s predictions, and we conclude with a summary and discussion in § 6. As our test problem, we study the dynamics of a one-dimensional, collisionless, electrostatic hydrogen plasma. We focus on the electron dynamics, treating the ions as a cold, neutralising background. The equations that describe this plasma are the Vlasov–Poisson system (Vlasov 1968). The calculations in this paper use Gaussian units. The Vlasov equation for electrons in one dimension is: (2.1) ∂f/∂t + v ∂f/∂z - (eE/m_e) ∂f/∂v = 0, where f(z,v,t) is the electron distribution function, e is the elementary charge and m_e is the mass of an electron. Term E(z,t) is the electric field, which we calculate using Poisson’s equation: (2.2) -∂²Φ/∂z² = 4πρ, E = -∂Φ/∂z, where ρ(z,t) is the total electric charge density, including both ions and electrons, and Φ(z,t) is the electric potential. For this quasi-neutral plasma with cold ions, we can evaluate the charge density as: (2.3) ρ(z,t) = e ( n_0 - ∫ dv f ), where n_0 is the mean number density of both species. Plasmas with near-Maxwellian velocity distributions are of particular interest for magnetic fusion applications. Therefore, we separate the distribution function into a Maxwellian mean and time-dependent perturbations: (2.4) f(z,v,t) = F_0(v) + g(z,v,t), where (2.5) F_0(v) = (n_0 / v_te) (1 / √(2π)) e^(-v² / 2v_te²) is the one-dimensional Maxwellian velocity distribution for the electrons, and g(z,v,t) is the fluctuating part of the distribution function. We have defined the electron thermal speed as v_te = √(T_e / m_e), where T_e is the temperature of the electrons in units of energy. The Vlasov equation supports longitudinal Langmuir waves, whose dynamics occurs on time scales of the order of the plasma frequency: (2.6) ω_pe ≡ √(4π n_0 e² / m_e). The characteristic length scale for those oscillations is the Debye length: (2.7) λ_D ≡ √(T_e / 4π n_0 e²). To support numerical integration of this system of equations, we normalise the variables into dimensionless forms, denoted with primes: (2.8) t ≡ (1/ω_pe) t', v ≡ v_te v', z ≡ λ_D z', E ≡ (T_e / eλ_D) E', f ≡ F_0 + g ≡ (n_0 / v_t)(F'_M + g'), ρ ≡ n_0 e ρ', Φ ≡ (T_e / e) Φ', F_M' ≡ (1/√(2π)) e^(-v'²/2), ∫ dv' F_M' = 1. After dropping the primes for legibility, we can now write the model equations in dimensionless form: (2.9) ∂f/∂t + v ∂f/∂z - E ∂f/∂v = 0, (2.10) ∂E/∂z = 1 - ∫ dv f. Finally, after splitting the distribution function into its mean and fluctuating components, we arrive at the normalised evolution equation for the fluctuating part of the electron distribution function: (2.11) ∂g/∂t + v ∂g/∂z + v E F_M - E ∂g/∂v = 0, (2.12) ∂E/∂z = -∫ dv g. We solve the Vlasov–Poisson system (2.11)–(2.12) using an orthonormal spectral basis by projecting the velocity-space dependence onto a Hermite polynomial basis and the configuration-space dependence onto Fourier harmonics. We choose this approach because pseudo-spectral methods are often more efficient than pure finite-difference algorithms (Boyd 2000), and as mentioned in the introduction, there is an extensive history of the application of the Fourier–Hermite basis to solve problems in kinetic plasma physics. 3.1. Hermite moment formalism To discretise velocity space, we impose a Hermite basis representation, using notation from Mandell et al. (2018). Including both the Maxwellian weight and the normalisation factors, the Hermite basis expansion (φ^m(v)) and projection (φ_m(v)) functions are: (3.1) φ^m(v) = (He_m(v) / √(2π m!)) e^(-v²/2), φ_m(v) = He_m(v) / √(m!), where (3.2) He_m(v) = (-1)^m e^(v²/2) (d^m/dv^m) e^(-v²/2) are the probabilists’ Hermite polynomials (Abramowitz & Stegun 1972). Note that these functions are orthonormal: (3.3) ∫_{-∞}^{∞} dv φ^m(v) φ_{m'}(v) = δ_{m,m'}, and that the background Maxwellian is identical to the m=0 basis element: (3.4) φ^0 = (1/√(2π)) e^(-v²/2) = F_M. This property makes the Hermite moment representation a natural choice for near-Maxwellian velocity distributions. We project the fluctuating portion of the distribution function onto the Hermite basis using: (3.5) g(z,v,t) = ∑_{m=0}^∞ φ^m G_m(z,t), G_m(z,t) ≡ ∫_{-∞}^{∞} dv φ_m g(z,v,t). 3.2. Vlasov equation in the Hermite basis We seek an evolution equation for the amplitudes, G_m, of the Hermite expansion functions. To derive this moment hierarchy, we project (2.11) onto the Hermite basis term by term: (3.6) ∫_{-∞}^{∞} dv φ_m ∂g/∂t = ∂G_m/∂t. We then project the linear ballistic term, v(∂g/∂z), onto the basis: (3.7) v ∂g/∂z = (∂/∂z) ∑_{m=0}^∞ v φ^m G_m = (∂/∂z) ∑_{m=0}^∞ (√(m+1) φ^{m+1} + √m φ^{m-1}) G_m, where we used a recurrence relation for the Hermite basis functions to remove the explicit dependence on v (Abramowitz & Stegun 1972): (3.8) √(m+1) φ^{m+1} = v φ^m - √m φ^{m-1}. Applying the projection function to (3.7) and using (3.3) yields: (3.9) ∫_{-∞}^{∞} dv φ_m v ∂g/∂z = (∂/∂z) (√m G_{m-1} + √(m+1) G_{m+1}). The linear portion of the electric field drive term contains the m=1 expansion function, so applying the projection function to this term yields: (3.10) ∫_{-∞}^{∞} dv φ_m v E F_M = E ∫_{-∞}^{∞} dv φ_m φ^1 = E δ_{m,1}. The final term to consider is the nonlinear term, -E(∂g/∂v). Combining a second recurrence relation (Abramowitz & Stegun 1972): (3.11) (d/dv) He_m(v) = v He_m(v) - He_{m+1}(v), with the definition of the Hermite basis functions φ^m (3.1), presents a useful new identity: (3.12) ∂φ^m/∂v = -√(m+1) φ^{m+1}. The expansion of the nonlinear term in the Hermite basis is then: (3.13) -E ∂g/∂v = -E ∑_m (∂φ^m/∂v) G_m = E ∑_m √(m+1) φ^{m+1} G_m. After projecting this term onto the basis: (3.14) ∫_{-∞}^{∞} dv φ_m (E ∑_{m'} √(m'+1) φ^{m'+1} G_{m'}) = ∑_{m'} E √(m'+1) δ_{m,m'+1} G_{m'} = E √m G_{m-1}, we can collect all of the terms to recover the general form of the equation for G_m: (3.15) ∂G_m/∂t + (∂/∂z)(√m G_{m-1} + √(m+1) G_{m+1}) + E √m G_{m-1} + E δ_{m,1} = 0. 3.3. Pseudo-spectral method in Fourier space Next, we discretise space with a Fourier basis. The Fourier transform and inverse Fourier transform operators for a spatial grid of N points are: (3.16) F[g] = (1/N) ∑_z g e^(-ikz), F⁻¹[G_k] = ∑_k G_k e^(ikz), where the Fourier wavenumbers are given by k = 2πj / N, with -N/2 ≤ j < N/2 and j ∈ ℤ. The Fourier convolution theorem (3.17) transforms the nonlinear term into a coupled sum over all wavenumbers. This coupling introduces additional computational overhead. We mitigate that problem by using a pseudo-spectral method (Boyd 2000), evaluating the nonlinear term on the original configuration space and projecting the result back onto the space of Fourier wavenumbers. Applying the Fourier transform to (3.15) and using the Fourier convolution theorem yields: (3.18) ∂G_{m,k}/∂t + ik(√m G_{m-1,k} + √(m+1) G_{m+1,k}) + √m F[ F⁻¹[E_k] F⁻¹[G_{m-1,k}] ] = -E_k δ_{m,1}. Simulating this collisionless system with finite resolution without using any regularisation is numerically infeasible, as the dynamics can generate structures in phase space at arbitrarily small scales (Parker & Dellar 2015; Loureiro et al. 2016; Mandell et al. 2024). To mitigate that effect, we introduce numerical regularisation through hyperdiffusion (-ν_k k⁴) and hypercollision (-ν_m m⁴) operators: (3.19) H_{mk} = 0 if m ≤ 2; -ν_k k⁴ - ν_m m⁴ if m > 2, where ν_k and ν_m are parameters which control the strengths of the hyperdiffusion and hypercollisions, respectively. We solve the system for a finite number of moments M. Defining compact notation for the nonlinear operator: (3.20) N[E_k, G_{m-1,k}] = √m F[ F⁻¹[E_k] · F⁻¹[G_{m-1,k}] ], we derive the final system equations (3.21)–(3.24): (3.21) ∂G_{0,k}/∂t + ik G_{1,k} = 0, (3.22) ∂G_{1,k}/∂t + ik (G_{0,k} + √2 G_{2,k}) + N[E_k, G_{0,k}] + E_k = 0, (3.23) ∂G_{m,k}/∂t + ik (√m G_{m-1,k} + √(m+1) G_{m+1,k}) + N[E_k, G_{m-1,k}] = H_{m,k} G_{m,k}, for 2 ≤ m < M-1, (3.24) ∂G_{M-1,k}/∂t + ik √(M-1) G_{M-2,k} + N[E_k, G_{M-2,k}] = H_{M-1,k} G_{M-1,k}. 3.4. Poisson’s equation in the Fourier–Hermite basis To complete the system, we project Poisson’s equation onto the Fourier–Hermite basis and solve for the electric field, E_k: (3.25) ρ = -∫_v dv g = -∫_v dv φ_0 ∑_m φ^m G_m = -G_0, (3.26) E_k = i G_{0,k} / k. We require that the k=0 mode is zero and stationary, corresponding to zero mean electric field. 3.5. The closure problem A continuous representation formally requires an infinite number of Hermite moments, but we must discretise the system with a finite number of moments, M. Observe that in the Fourier–Hermite basis, the parallel streaming term in the convective derivative couples nearest neighbours in Hermite space and forms a hierarchy of equations that must be closed. The simplest model of closure is by truncation, which we use in the direct numerical solution (DNS) by omitting the G_{m+1,k} term in the evolution equation for the final moment, G_{M-1,k}. Our objective is to introduce a dynamical closure model to close the hierarchy at a moment number m = m_c, m_c < M-1, while still achieving accurate predictions of the Fourier and Hermite spectra. Inspecting the evolution equation for G_{m,k} (3.15) reveals that the coupling between moment m and moment m+1 only explicitly appears as a linear term in the equation, acting on each wavenumber k independently. As a result, we introduce a unique and independent Hermite closure model for each k in the system, for a total of N_k reservoirs. We use methods similar to those of Pathak et al. (2018) to design the reservoirs. Each reservoir is a network of D_r nodes defined by a latent state vector r_k(t) ∈ ℝ^(D_r). The connections between nodes are defined by a D_r × D_r sparse adjacency matrix, A, where exactly κ non-zero elements per row are initialised to random values drawn from the uniform distribution on [0, 1] and scaled to set the largest eigenvalue of A equal to a specified spectral radius ρ_sp. The runtime is separated into two phases: training and prediction. In the training phase, the inputs are derived from a high-resolution DNS. In the subsequent prediction phase, the inputs are derived from a low-resolution simulation with the output of our reservoir used to compute the closure. The input to a given reservoir at each time step is a vector u_k(t) of dimension 2w, where the elements are the real and imaginary parts of the final w Hermite moments before m_c: (4.1) u(t) = [ G_{m_c-w,k}(t-wΔt), G_{m_c-(w-1),k}(t-(w-1)Δt), ..., G_{m_c-1,k}(t-Δt) ]. The input layer is a D_r × 2w matrix W_in, with elements drawn uniformly from [-σ, σ]. Matrices A and W_in remain fixed after initialisation. Only the output layer, W_out, is trained. The reservoir state vector evolves according to: (4.2) r_k(t + Δt) = tanh[ A r_k(t) + W_in u(t) ]. To account for quadratic nonlinearity which breaks odd-parity symmetry, we apply a symmetry-breaking transformation r_k -> r_k^* to the reservoir latent state before readout: (4.3) r_{i,k}^*(t) = r_{i,k}(t) for even i; r_{i,k}^2(t) for odd i. We train W_out using Tikhonov-regularised linear regression to minimise: (4.4) ∑_{j=b}^{T/Δt} ||G_{m_c,k}(jΔt) - W_out r_k^*(jΔt)||² + β ||W_out||², which yields the closed-form solution: (4.5) W_out = G_{m_c,k} R_k^{*T} ( R_k^* R_k^{*T} + β I )⁻¹. 4.1. Closure integration During the prediction phase, the system does not evolve moments beyond m = m_c. For each wavenumber k, resolved moments m < m_c are integrated using operator D: (4.6) G_{m < m_c, k}(t + Δt) = D[ G_{m ≤ m_c}(t) ]_k, and the reservoir predicts the closure moment: (4.7) G_{m_c, k}(t + Δt) = W_out r_k^*(t + Δt). 4.2. Computational complexity For the DNS, the per-step cost scaling for evolving the higher moments is: (4.8) Cost_DNS = 3 M' O(N_k log N_k) + 20 M' N_k + N_k, where M' = M - m_c. For the ML closure, the cost per time step is: (4.9) Cost_ML = N_k D_r [8w + 4κ + 13/2]. For sufficiently large (M' log N_k) / (D_r w), the reservoir closure is computationally faster than the DNS. 5.1. Linearised Vlasov–Poisson For our first test, we verify the model on the linearised Vlasov–Poisson system with a single wavenumber (k=0.4) and external driving: (5.2) ∂G_{m,k}/∂t + ik(√m G_{m-1,k} + √(m+1) G_{m+1,k}) = χ(t) δ_{m,2} - E_k δ_{m,1}, where χ(t) is drawn from U[0, 1). Energy cascades down to finer scales through Landau damping. In the absence of regularisation, numerical reflection occurs at the high-m boundary in finite-resolution simulations (e.g. M=1025). When closing the system at m_c = 3 using a reservoir with D_r = 6 nodes, the ML closure achieves excellent agreement with the high-resolution baseline, with a mean-squared error of 2.25 × 10⁻⁵ in the Hermite spectrum, significantly outperforming naive truncation and hypercollisional damping. 5.2. Weakly nonlinear Vlasov–Poisson We reintroduce the nonlinear term and simulate an initial-value problem with an initial cosine density perturbation: (5.3) g(z,v,t=0) = ε cos(k_0 z) F_M(v), (5.4) G_{m,k,t=0} = (ε/2) δ_{m,0} (δ_{k,k_0} + δ_{k,-k_0}), with ε = 0.001 (0.1% perturbation), k_0 = 0.4, N_k = 17, and M = 17 for the baseline. Closing at m_c = 3 (reducing velocity resolution to M = 4) using independent 12-node reservoirs for each k accurately reproduces the linear Landau damping rate and wave frequency, yielding mean-squared errors of 1.26 × 10⁻⁴² (Hermite spectrum) and 8.01 × 10⁻³⁸ (Fourier spectrum). 5.3. Strongly nonlinear Vlasov–Poisson For the strongly nonlinear regime, we increase the perturbation amplitude to ε = 0.18. At this high amplitude, the damping attenuates and the mode saturates. The baseline requires M = 65 moments to resolve the low-order spectrum accurately. Applying the ML closure at m_c = 3 (M = 4) using 120-node reservoirs achieves mean-squared errors of 3.66 × 10⁻¹⁸ for the Hermite spectrum and 1.19 × 10⁻¹⁴ for the Fourier spectrum, maintaining accurate wave frequency and capturing both spectra while reducing the required velocity-space resolution by a factor of 16. In this paper, we have used reservoir computing to present an ML, velocity-space closure model for the hierarchy of Hermite moments in the one-dimensional Vlasov–Poisson system. The closure model serves as a proof-of-principle that reservoir computing can be used to reduce the simulation domain requirements for studying plasma dynamics. A logical objective to pursue is to extend the closure to reduce the required resolution in the Fourier representation of position space. In practice, gyrokinetic turbulence simulations commonly use significantly more resolution for configuration space than they do for velocity space, in part because of the near-Maxwellian nature of the velocity distributions. Therefore, a Fourier-space closure model would potentially be more desirable than a velocity-space closure. However, the structure of the nonlinear term poses a challenge to achieving that goal due to non-local convolutions across spatial scales. The eventual goal is to build closure models that can reduce the domain requirements for higher-dimensional turbulence simulations, including the full gyrokinetic equation. Future work should prioritize testing this velocity-space closure in spectral gyrokinetic codes such as Gyacomo and GX, as well as exploring methods to extract compact, interpretable symbolic representations of the closures. In § 5.3, a convergence study was conducted to determine the high-resolution (M=65) baseline for the strongly nonlinear case. The low-order moments represent important physical quantities (density, momentum, and energy). The lowest Hermite resolution that captures the m=0 and m=1 moments within a factor of 2 is M=65. Comparison with an M=513 reference confirms exponential error convergence for M ≥ 65. Time traces of the m=1 and m=2 moments for both weakly and strongly nonlinear initial-value problems confirm that the ML closure accurately captures the oscillation frequencies and peak amplitudes. Calculated damping rates show exact agreement in the weakly nonlinear regime and small deviations in the strongly nonlinear regime. We express our gratitude to R. Gaur, N.R. Mandell and M. Nastac for helpful discussions about turbulence, M. Almanza, E.P. Alves, D. Carvalho, J. Chou, F. Fiuza, G. Guttormsen, B. Jang, N.F. Loureiro, M. McGrae-Menge, J. Pierce and A. Velberg for their perspectives on subgrid modelling and E. Ott, D. Patel and A. Wikner for sharing their insights on reservoir computing. Editor Nuno Loureiro thanks the referees for their advice in evaluating this article. This work was made possible through the support of the United States Department of Energy under contract numbers DE-FG02-93ER54197 and DE-SC0022039. M.L. is a consultant for Type One Energy.