## 5 Application on the Lorenz-96 Model  
We apply the aforementioned methodology in the Lorenz-96 equations [25], which model the evolution of atmospheric quantities, such as vorticity or temperature at a discrete periodic lattice representing a latitude circle on earth. The model is given by the system of \(N\) equations  
\[\frac{dX_{i}}{dt}=(X_{i+1}-X_{i-2})X_{i-1}-X_{i}+F, \tag{31}\]  
where \(i=1,\ldots,N\). In this system, \(-X_{i}\) is a damping term, \(F\) is an external force, and the quadratic terms resemble the advection term and conserve kinetic energy. For more details on the solution of (31) and the treatment of \(X_{-1},X_{0}\), and \(X_{N}\), the reader is referred to References [25; 26]. This model is excellent for testing, since it is exhibits different dynamics depending on \(F\); thus, it is ergodic, and its behavior displays similarities with turbulent fluid flows [26]. In Figure 2, we plot the maximum Lyapunov exponent of the system as a function of the forcing term \(F\). Notice that the exponents are positive in the range of \(F\) considered, indicating that the Lorenz-96 model displays chaotic behavior. The exponents were computed using QR decomposition, as in Reference [27], and match well with the results from Reference [28]. The QoI considered is the time-averaged integral of the average value of \(X_{i}\), defined as  
\[J(F)=\frac{1}{NT}\sum_{i=1}^{N}\int_{0}^{T}X_{i}dt. \tag{32}\]Uncertainty is introduced through \(F=10+0.3\zeta\), where \(\zeta\sim\mathcal{N}(0,1)\) is Gaussian. To choose an appropriate chaos order \(c\), we perform a spectral convergence analysis for the first four moments of Equation (32), which is shown in Figure 3. Notice that, for \(c=4\), all four moments stabilize to the value predicted by a Monte-Carlo simulation with \(N=10,000\) samples. In Figure 4, we compare the first four statistical moments of Equation (32) obtained using the Monte-Carlo and gPC for \(c=4\). Notice that the Monte-Carlo is in good agreement with the gPC,  
Figure 3: Spectral convergence of the first four moments of the Quantity of Interest (QoI) (32) in the Lorenz-96 system. Here, \(T=100\), \(F=10+0.3\zeta\), where \(\zeta\sim\mathcal{N}(0,1)\). Figure 2: Maximum Lyapunov exponent of the Lorenz-96 model for different values of the forcing term \(F\) (blue dots). The integration time is \(T=500\) and \(N=40\). The red dots represent the exponents reported in Reference [28]. indicating that this approach can be a valuable method for UQ in the time averages of chaotic systems. Note that the skewness and kurtosis are close to 0 and 3, respectively, as the values for Gaussian distribution. A possible explanation is that \(J(F)\) is the sum of random variables that follow the same distribution, and from the central limit theorem, it follows that the PDF of \(J(F)\) is Gaussian. The same principle applies in the UQ of the sensitivity of QoI (32) with respect to \(F\). This sensitivity is obtained using the MSS algorithm described in Section 2. To that end, Equation (31) is differentiated with respect to \(F\)  
\[\frac{dV_{i}}{dt}=X_{i-1}V_{i+1}+(X_{i+1}-X_{i-2})V_{i-1}-X_{i-1}V_{i-2}-V_{i}+ 1-\eta\left((X_{i+1}-X_{i-2})X_{i-1}-X_{i}+F\right), \tag{33}\]  
where \(V_{i}=\frac{dX_{i}}{dF}\). Equation (33) is used as a constraint in the MSS algorithm (refer to Equation (5c)). Results can seen in Figure 5, where the sensitivity obtained from MSS is compared with FD. Figure 4: Uncertainty Quantification (UQ) in the Lorenz-96 model (\(N=40\) and \(T=100\)), for the QoI (32). Comparison of statistical moments computed by Monte–Carlo and generalized Polynomial Chaos (gPC) for (\(N=40\) and \(T=100\)) \(c=4\). The convergence of the GMRes algorithm applied to system (14) is shown in Figure 6. The residuals are reduced by more than 8 orders of magnitude in 13 iterations, which demonstrates the effectiveness of the preconditioner. By coupling the gPC with the MSS, the statistics of the sensitivity of the QoI can be computed efficiently.