A method for preserving nominally-resolved flow patterns in low-resolution ocean simulations: Constrained dynamics

Inability of low-resolution ocean models to simulate many important aspects of the large-scale general circulation is a common problem. In the view of physics, the main reason for this failure are the missed dynamical effects of the unresolved small scales of motion on the explicitly resolved large scale motions. Complimentary to this mainstream physics-based perspective, we propose to address this failure from the dynamical systems point of view, namely, as the persistent tendency of phase space trajectories representing the low-resolution solution to escape the right region of the corresponding phase space, which is occupied by the reference eddy-resolving solution. Based on this concept, we propose to use methods of constrained optimization to confine the low-resolution solution to remain within the correct phase space region, without attempting to amend the eddy physics by introducing a process-based parameterisation. This approach is advocated as a novel framework for data-driven hyper-parameterisation of mesoscale oceanic eddies in non-eddy-resolving models. We tested the idea in the context of classical, baroclinic beta-plane turbulence model and showed that non-eddy-resolving solution can be substantially improved towards the reference eddy-resolving benchmark.


Introduction
The problem of reproducing large-scale flow structures in low-resolution ocean simulations is among the most challenging ones in ocean modelling, and this is mostly due to the lack of information from the small, unresolved scales.The mainstream approach to this problem is to use parameterisations, that is, mathematically simple and physically justified approximations of the key unresolved and under-resolved small-scale processes (e.g., Gent and Mcwilliams [13], Duan and Nadiga [11], Frederiksen et al. [12], Jansen and Held [15], Mana and Zanna [19], Cooper and Zanna [5], Grooms et al. [14], Berloff [1,2], Danilov et al. [10], Ryzhov et al. [21], Juricke et al. [16,17], Cotter et al. [6], Ryzhov et al. [22], Cotter et al. [7,8,9]).Despite the decades of research effort invested in this direction, the problem remains mainly unresolved for various reasons, the majority of which is due to inaccurate description of small-scales physics and energy transfers between scales.This work takes a new angle on this long-standing problem.Namely, we propose to look at the problem from the dynamical system point of view and consider the inability of the low-resolution model to reproduce the large-scale flow structures as the persistent tendency of phase space trajectories representing the lowresolution solution to escape the right region of the corresponding phase space, which is occupied by the reference eddy-resolving solution.Thus, instead of parameterising directly the action of unresolved scales onto the resolved ones, our approach is to constrain the low-resolution flow dynamics to the right region of the phase space.Two other methods based on this approach were developed and tested earlier.The first method computes the low-resolution solution as the phase space trajectory of the image point, which is dynamically governed by the reference solution [23].The second method reconstructs a dynamical system from the reference solution, and then explicitly uses this system to predict the low-resolutions circulation [24].In this study we propose a new method: the employed low-resolution model, which is integrated explicitly, is constrained by an optimization procedure that restricts the solution to stay within the right region of phase space, as defined from the reference solution.For the proof-of-concept stage, we consider, first, the Lorenz model, and, then, the classical, baroclinic quasigeostrophic (QG) model of the beta-plane turbulence.

The method
We run the mathematical model of the studied phenomenon (here, the QG model) at low resolution and restrict its solution to the reference region of phase space defined by a set of constraints g (a ball of radius r in our case).The method is based on solving a constrained minimization problem of the form (e.g., Bertsekas [4]): min f (x), where f and g are given functions.In order to inject the constraint into the minimization problem, we use the barrier function φ = −1 g(x), which leads to the following unconstrained minimization problem: where parameter µ > 0 regulates accuracy for finding the minimum In our case, a very accurate approximation of the minimum is not required as the problem is chaotic, and the goal of constraining will be achieved anyway.Without compromising the main message, we leave further technicalities beyond the scope of our work and refer readers to Bertsekas [4].
As seen in Figure 1a, the solution to the Lorenz system remains on the attractor over the whole simulation time, and there are no constraints influencing the trajectory.If trajectory is to be restricted to some region of the phase space, then a set of constraints needs to be introduced.We chose a ball of radius r = 40 centred at (0, 0, 0) to be the constraining phase space region, and the constrained solution remains in it over the whole simulation (Figure 1b).But, the constrained solution is not the same as the original one even within the ball interior -the region where trajectory is allowed to evolve with no constraints as long as the constraints are satisfied.It happens because the constraints deform the structure of the original phase space by imposing extra restrictions on the solution.If we keep the radius decreasing, the deformation eventually becomes  so large that the constrained system completely fails to reproduce the Lorenz attractor (not shown).This indicates that the Lorenz system should be somehow modified, but discussing this is beyond the scope of our study.
The above example shows how one can use the constrained optimization approach in more sophisticated settings.For instance, one can take the primitive-equations or a QG ocean model, compute its eddy-resolving reference solution, project it onto the coarse grid, and find (approximately) a spherical region of the phase space that is occupied by this solution.If the unconstrained low-resolution model configuration cannot reproduce the (nominally resolved) reference circulation, then implementation of constraint is justified, and the questions are whether it works and how can it be tuned and optimized.
The layer-wise PV anomalies and streamfunctions are related through the pair of coupled elliptic equations: with stratification parameters s 1 = 4.22 ⋅ 10 −3 km −2 and s 2 = 1.41 ⋅ 10 −3 km −2 chosen so that the first baroclinic Rossby deformation radius is Rd 1 = 25 km.The mass and momentum constraints are imposed following McWilliams (1977).System (4)-( 6) is augmented by the periodic horizontal boundary conditions set at the eastern (Γ 2 ) and western (Γ 4 ) boundaries: and no-slip boundary conditions, u In order to apply the method to the QG equations ( 4) one has to choose the variable or variables to constrain.For this role we take the PV anomaly q that leads to the inequality g(q) ≤ 0; however, the method is not limited to this particular choice, and other constraints can be implemented.For example, one can constrain the velocity field or individual terms in the equations.
As seen in Figure 2a, in each direction the reference solution has 4 well-pronounced horizontal jets, which the (unconstrained) low-resolution model fails to reproduce (Figure 2b).This is because the low-resolution solution trajectory quickly escapes the right phase space region, even when we start it from the right initial condition.The phase space is constrained as where r = 846 (non-dimensional units) and T = 4 years.Constraint (9) corresponds to the ball of radius r centered at the 4-year time mean of the reference eddy-resolving solution ⟨q⟩; and the radius is found as the mean distance of the 4-year long reference solution from ⟨q⟩.Here, the ball is the simplest approximation of the real reference-solution attractor in the phase space.A more accurate approach would be to find an hyper-ellipsoid that takes into account the actual spread of the attractor along different coordinate axes.
When the model is constrained, the 4 jets are recovered both in the instantaneous and time-mean fields (Figure 2c).The relative blurriness of the jet edges comes from the use of an order-of-magnitude larger viscosity.Perhaps, the jet edges can be sharpened up by injecting noise into the advection operator, as in the SALT approach (e.g.Cotter et al. [6,7,8,9]), or simply by adding it as a forcing.To see how the radius of the ball influences the dynamics, we took the ball of a larger radius r = 923 (this is the maximum radius for the 10-year-long reference solution).In this case, the constrained solution still has the jets, which are, however, noticeably corrupted (Figure 2d), especially after 10 years of simulation (see the jet near the southern boundary).Moreover, the separation between the two southern jets becomes less evident (see the time-average subplot) and the standard deviation does not show as pronounced jets as those of the reference solution or the solution constrained in the smaller ball.This is explained by the fact that the solution constrained by a larger ball drifts farther away from the right phase space region.Based on this evidence, one could think that a tighter ball would yield a more accurate solution, but this is not necessarily the case, since an over-constrained QG model can lead to significantly incorrect dynamics.As an example, we found the solution for r = 796 (this is the minimum radius for the 10-year-long reference solution) and obtained the results very similar to those with r = 846 (not shown).This suggests that the ball is likely not an accurate approximation of the reference phase space when finer structures of the solution have to be modelled.As an alternative, one should focus on methods that can better approximate the reference phase space.

Conclusions and discussion
In this work we have further developed alternative hyper-parameterisation approach for parameterising effects of mesoscale oceanic eddies on the large-scale ocean circulation.Complimentary to the mainstream physics-based perspective, we propose to deal with the eddy effects from the dynamical systems point of view, and interpret the lack of them as the persistent tendency of phase space trajectories representing the low-resolution solution to escape the right region of the corresponding phase space, which is occupied by the reference eddy-resolving solution.Based on this concept, we propose to use methods of constrained optimization to confine the low-resolution solution to remain within the correct phase space region, without attempting to amend the eddy physics by introducing a process-based parameterisation.We used the Lorenz 63 system as a simple conceptual toy model for the proof of concept.Then, we considered the baroclinic, quasigeostrophi (QG) model of the beta-plane turbulence and showed that the low-resolution solution cannot reproduce such nominally-resolved flow structures, as the multiple alternating zonal jets, which are robustly present in the reference eddy-resolving solution.Implementation of the proposed hyper-parameterisation approach significantly improves the low-resolution model and recovers the jets.In the reference phase space we used a ball, with the radius found as the time-mean distance of the 4-year-long reference solution from its time mean, and centered at the 4-year time mean of the reference solution.
It is worth drawing reader's attention to the fact that only 4 years of the reference solution have been used to find the centre of the ball and its radius, and this turned out to be enough to reproduce the jets in the 10-year-long run at low resolution.This demonstrates that the proposed method gives robust results, as soon as it is supplied with sufficient data.On the other hand, for more sophisticated models or even for different setups of the QG model, longer or shorter runs might be needed to accurately estimate the centre and radius.
In addition, we have studied how the radius of the ball influences the ability of the constrained QG model to reproduce nominally-resolved structures and found that too large radius eventually lead to the solution degradation, due to the detrimental drift in the phase space, whereas too small radius does not allow for an accurate approximation of the shape of the reference phase space.Thus, there is the optimal strength of the constraint.
The utility of the proposed method is that it does not require in-depth physical knowledge of interactions between the scales of motion, and it does not require any modification of the governing equations, as in the case of traditional parameterisations.The method falls into the category of data-driven methods as it requires either detailed observations or their substitute in terms of some eddy-resolving solution data.Obviously, its advantage may turn into a possible drawback, as often information presented in data is not sufficient (for example, characteristics of the constraining ball may be inaccurately estimated).This can be mitigated by using more accurate approximations of the constraining geometry (for example, a hyperellipsoid instead of the ball to take into account the spread of the solution along different coordinate axes) and of the attractor reconstruction.Another intriguing avenue for future research extension is the so called term-wise constraining, which deals with individual dynamical terms, rather than with the whole solution, like in this study.
The proposed method can be implemented into the dynamic core of oceanic general circulation models by constraining solutions of different equations and terms or their combinations.This is, in turn, the open broad research agenda on the effects of different data-driven constraints of the governing equations (e.g., Shevchenko and Berloff [25]).
is zonally periodic flat-bottom channel with horizontal dimensions L x = 1800 km and L y = L x 2, with the depth H = H 1 + H 2 , and filled out by two stacked isopycnal fluid layers of depths H 1 = 1.0 km and H 2 = 3.0 km.
are imposed at the northern (Γ 1 ) and southern (Γ 3 ) boundaries.We use 513 × 257 uniform spatial grid, take the eddy viscosity value ν = 25 m 2 s −1 , and spin up the model from the state of rest to t = 0 over the time interval T spin = [−10, 0] years, so that the statistically equilibrated flow regime is established.For the low-resolution model, we use the same set-up, except for much coarser grid 129 × 65 and much larger eddy viscosity value ν = 250 m 2 s −1 .

Figure 2 :
Figure 2: Shown are snapshots of the top-layer PV anomaly and standard deviation of: (a) the reference solution (computed with ν = 25 m 2 s −1 on the grid 513 × 257 and then projected on the grid 129 × 65); (b) the unconstrained low-resolution solution (computed on the coarse grid 129 × 65 with ν = 250 m 2 s −1 ); (c) the constrained low-resolution solution restricted to the ball with radius r = 846 (non-dimensional units); (d) the same as (c) but for r = 923.Units of the PV anomaly fields are nondimensional; A = 30, B = 10.Note that the constrained solution (c) preserves the nominally-resolved flow structures (the 4 jets) both in the snapshots and time mean, whereas the unconstrained one fails on this (b).