  We propose a statistical-stochastic surrogate modeling approach to predict
the response of the mean and variance statistics under various initial
conditions and external forcing perturbations. The proposed modeling framework
extends the purely statistical modeling approach that is practically limited to
the homogeneous statistical regime for high-dimensional state variables. The
new closure system allows one to overcome several practical issues that emerge
in the non-homogeneous statistical regimes. First, the proposed ensemble
modeling that couples the mean statistics and stochastic fluctuations naturally
produces positive-definite covariance matrix estimation, which is a challenging
issue that hampers the purely statistical modeling approaches. Second, the
proposed closure model, which embeds a non-Markovian neural-network model for
the unresolved fluxes such that the variance of the dynamics is consistent,
overcomes the inherent instability of the stochastic fluctuation dynamics.
Effectively, the proposed framework extends the classical stochastic parametric
modeling paradigm for the unresolved dynamics to a semi-parametric
parameterization with a residual Long-Short-Term-Memory neural network
architecture. Third, based on empirical information metric, we provide an
efficient and effective training procedure by fitting a loss function that
measures the differences between response statistics. Supporting numerical
examples are provided with the Lorenz-96 model, a system of ODEs that admits
the characteristic of chaotic dynamics with both homogeneous and inhomogeneous
statistical regimes. In the latter case, we will see the effectiveness of the
statistical prediction even though the resolved Fourier modes corresponding to
the leading mean energy and variance spectra do not coincide.
