{'fraction_non_alphanumeric': 0.042410037783120844, 'fraction_numerical': 0.014830380053744974, 'mean_word_length': 4.607749830047586, 'pattern_counts': {'":': 1, '<': 6, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal Closures in a Simple Model for Turbulent Flows

In this work we introduce a computational framework for determining optimal closures of the eddy-viscosity type for Large-Eddy Simulations (LES) of a broad class of PDE models, such as the Navier-Stokes equation. This problem is cast in terms of PDE-constrained optimization where an error functional representing the misfit between the target and predicted observations is minimized with respect to the functional form of the eddy viscosity in the closure relation. Since this leads to a PDE optimization problem with a nonstandard structure, the solution is obtained computationally with a flexible and efficient gradient approach relying on a combination of modified adjoint-based analysis and Sobolev gradients. By formulating this problem in the continuous setting we are able to determine the optimal closure relations in a very general form subject only to some minimal assumptions. The proposed framework is thoroughly tested on a model problem involving the LES of the 1D Kuramoto-Sivashinsky equation, where optimal forms of the eddy viscosity are obtained as generalizations of the standard Smagorinsky model. It is demonstrated that while the solution trajectories corresponding to the DNS and LES still diverge exponentially, with such optimal eddy viscosities the rate of divergence is significantly reduced as compared to the Smagorinsky model. By systematically finding optimal forms of the eddy viscosity within a certain general class of closure models, this framework can thus provide insights about the fundamental performance limitations of these models.


Introduction and Problem Statement
Turbulent flows at high Reynolds numbers continue to challenge both scientists studying their fundamental properties and engineers interested in diverse technical applications involving fluid mechanics. In particular, accurate and efficient numerical simulation of turbulent flows will for the foreseeable future remain an open problem in computational science, despite advances in both algorithms and computer architectures. This is because the solutions of the three-dimensional (3D) Navier-Stokes equation, which is the main mathematical model assumed to describe the motion of viscous incompressible fluids, are chaotic and exhibit extreme spatio-temporal complexity at Reynolds numbers characterizing developed turbulence. More specifically, such flows are known to involve complex interactions of motions at vastly different time-and length-scales [1] which cannot be simultaneously resolved even with the largest available computational resources. Under the assumptions of spatial homogeneity and isotropy, it is convenient to consider solutions to the Navier-Stokes system as energy spectra in terms of their Fourier-space representation. In this context a standard paradigm for describing turbulent flows is expressed in terms of Kolmogorov's K41 theory [2] and its various refinements [1], which stipulate that in the statistically stationary setting the energy is injected into the flow at large length scales of order L > 0 and is dissipated at some small length scales 0 < η ≪ L, where η is typically referred to as the Kolmogorov length scale. The energy is then cascaded through the range of intermediate wavenumbers L −1 k η −1 via nonlinear interactions. With the Reynolds number defined as Re = UL/ν N , where U is the characteristic velocity and ν N is the coefficient of the kinematic viscosity of the fluid, simple dimensional arguments lead to the conclusion that the number of discrete degrees of freedom necessary to resolve the flow down to the Kolmogorov length η scales as O(Re (9/4) ) [1]. This hints at fundamental limitations on the largest Reynolds numbers for which direct numerical simulation (DNS) can be performed on the 3D Navier-Stokes system.
An approach which allows one to get around the aforementioned difficulty and obtain approximate solutions of the flow problem in some practical situations relies on the socalled Large-Eddy Simulation (LES) in which one solves a suitably filtered version of the governing equations. To define this approach, we consider the 3D Navier-Stokes equation (for brevity, we omit here the specification of the flow domain Ω together with the initial and boundary conditions which are assumed generic) where v = [v 1 , v 2 , v 3 ] T and p are the velocity and pressure fields, respectively, whereas ρ denotes the density (here and below, Einstein's convention is used with repeated indices implying summation). Given the linear filtering operation v(x):= (G δ ⋆ v)(x) = R + → R + in which δ > 0 denotes the cut-off length scale and the symbol ":=" means "equal to by definition", the LES formulation of the Navier-Stokes system is obtained by applying this filtering operation to (1) and takes the form where M ij := v i v j − v i v j , i, j = 1, 2, 3, is by analogy with the dissipative term already present in the Navier-Stokes system referred to as the "subgrid-scale" (SGS) stresses [3] (the subgrid stresses are closely related to the Reynolds stresses arising in the context of the Reynolds-averaged Navier-Stokes (RANS) equations). System (2) describes the evolution of the filtered (large-scale) velocity field v and, evidently, is not closed because the SGS stresses depend on the original (unfiltered) velocity field v. Since the filtering operation defined by G δ typically involves attenuation of velocity components with length scales smaller than δ, the SGS stresses thus represent the averaged effect of these neglected motions on the evolution of the resolved flow field v. In order to close system (2) one therefore needs to represent the SGS stresses in terms of the resolved field v in some way, which constitutes the celebrated "turbulence closure problem" [4].
There is a very large body of results concerning the closure problem formulated in different flow conditions, especially in the engineering literature. Even a brief survey of these results would be outside the scope of the present study and we refer the reader to the monographs [4,5] for more information. Arguably, the most commonly used family of closure models is of the eddy-viscosity type in which the SGS stresses are expressed as [6] M ij = ν S ij , where S ij := 1 2 in which ν is the eddy viscosity (to simplify the notation used below and in contrast to the commonly employed convention, we choose to adopt a simple symbol for the eddy viscosity and put a subscript on the kinematic viscosity used in (1a)). The approaches to determining this quantity can be classified as algebraic, in which some algebraic relation is postulated between the filtered field v and the eddy viscosity ν (such as the celebrated Smagorinsky model discussed below), and differential, in which the eddy viscosity is assumed to depend on some additional quantities whose transport is described by suitable partial-differential equations (PDEs), such as, e.g., the fluctuating kinetic energy k and dissipation ǫ in the k − ǫ model. Since the SGS stresses are assumed in (3) to depend on the strain field S ij , the eddy-viscosity models have a similar structure to the dissipative term already present in the Navier-Stokes equation (1a), so its inclusion in the equation has the effect of changing the coefficient of this term from ν N to (ν N +ν). We note however that, unlike the kinematic (molecular) viscosity coefficient ν N which is constant, the eddy viscosity ν depends on the resolved field v and therefore introduces an additional nonlinearity. In addition to the classical Smagorinsky model, there exist many other approaches to approximate the eddy viscosity, including dynamic Smagorinsky models [7] relying on Germano's commutator identity [8] and the structure-function models [9], to mention just two. For a survey of recent progress in the field of turbulence modelling we refer the reader to [10], whereas some recent algorithmic developments in the solution of the LES equations (2) are discussed in [11,12,13]. Regardless of details, in all cases these closure models are postulated based on empirical grounds, albeit usually with a strong physical justification, with a small number of parameters requiring calibration from data. Since most closure models are derived in a heuristic manner, such approaches to turbulence modelling have been sometimes criticized as lacking scientific rigor and therefore difficult to generalize to flows different from the ones for which they have been calibrated. The objective of the present study is to provide insights about how well eddy-viscosity models can in principle perform. This is achieved by finding, via solution of a suitable optimization problem, a mathematically optimal form of the eddy viscosity for a given flow. More precisely, while in algebraic closure models a simple relationship is typically postulated for the dependence of the eddy viscosity ν on the resolved flow field v involving a small number of free parameters, in our proposed approach we will determine the functional form of this dependence optimally in a very general setting subject only to some minimal assumptions.
To fix attention, we will consider what is arguably the most common algebraic closure model, namely, the Smagorinsky model postulating the eddy viscosity in the form ν = in which C s is an adjustable constant known as the Smagorinsky coefficient [6]. Although the Smagorinsky model is rather simple, it is quite popular and serves as the "workhorse" for many LES computations. It is known, however, to possess certain deficiencies which include being too dissipative near walls [3]. Another major drawback of the Smagorinsky model is that it assumes the eddy viscosity to be zero if the resolved strain S vanishes and the fact that the eddy viscosity is positive otherwise, implying that the closure is strictly dissipative [14]. In our investigation we introduce a computational framework allowing one to determine an optimal Smagorinsky model in which the eddy viscosity is allowed to have a very general functional dependence on the magnitude of the resolved strain field 2 S ij S ij 1/2 which is found by matching the predictions of the LES model against a given "target" field (e.g., obtained by solving the original Navier-Stokes problem (1) via DNS or from an experiment). Since this eddy viscosity is optimal within the class of Smagorinsky-type models, its properties will provide insights about how well this class of models can in principle perform.
There have been earlier attempts to determine turbulence closure models with some optimality properties. Langford & Moser [15] and then Das & Moser [16] developed an approach for isotropic turbulence in which subgrid scales were treated as stochastic and the closure was determined by minimizing the modelling error using stochastic estimation techniques. This approach was tested on a range of models, including a stochastically forced one-dimensional (1D) Burgers equation and 3D Navier-Stokes system. An emerging family of approaches uses various machine-learning techniques such as neural networks to deduce closure models with certain optimality properties from data. In this context we mention the investigations [17,18,19,20,21,22], whereas the state-of-the-art in this field is discussed in the review papers [23,24,25]. Datadriven machine-learning methods have also been utilized for computational prediction and diagnosis of various turbulent flows [26,27,28]. Closure models obtained with some other data-driven methods are discussed in [29,30,31]. We note that, while the approach developed in the present study can also be classified as "data-based", it does not rely on machine learning, but rather on the calculus of variations and rigorous methods of PDE-constrained optimization. More specifically, recognizing that closure models are in fact forms of constitutive relations, we adapt the method developed in [32,33] to infer optimal constitutive relations from data. In the context of hydrodynamics such techniques have already been used to tackle the simpler problem of finding optimal closures for finite-dimensional reduced-order models in [34,35] and in vortex dynamics [36]. Applications of this approach in the field of electrochemistry are discussed in [37].
Since the goal of this paper is to provide a "proof of the concept" for the proposed approach by introducing and validating it in a general context, we shall focus here on a 1D model problem which will allow us to avoid the technical complexities inherent in dealing with the 3D Navier-Stokes system (1). We will thus consider the Kuramoto-Sivashinsky equation defined on the periodic domain Ω = [0, 2π] where T > 0 is the length of the time window considered, ν 4 , ν 2 ∈ R + are parameters whereas w 0 is an appropriate initial condition. The reason for choosing this particular model problem is that its solutions are known to exhibit important features characteristic of actual turbulence governed by the 3D Navier-Stokes system, namely, chaotic and multiscale dynamics with significant spatio-temporal complexity. These properties arise from an interplay between the linear and nonlinear terms in (4a): the second-order negative diffusion term injects energy at low wavenumbers which is then transferred by the nonlinear interactions to higher wavenumbers where it is eventually dissipated by the fourth-order dissipative term. Unlike the Burgers equation, in which a similar behavior may only arise from the inclusion of a somewhat artificial stochastic forcing term [38], the Kuramoto-Sivashinsky equation intrinsically exhibits a more turbulence-like behavior. Originally this equation was proposed as a model for instabilities of interfaces and flame fronts [39], and "phase turbulence" in chemical reactions [40]. Beyond its original purpose, this equation has been used as a model for hydrodynamic turbulence and is commonly employed as a testbed to study new approaches [41].
The structure of the paper is as follows: in the next section we introduce the optimal closure problem in the context of the Kuramoto-Sivashinsky equation (4); its solution based on a gradient approach is described in a general context in Section 3, whereas the set-up of the particular problem considered here is described in Section 4; details of the numerical approach are presented in Section 5; our computational results are discussed in Section 6, whereas a summary and final conclusions are deferred to Section 7; some technical material is collected in an appendix.

Eddy-Viscosity Closures for the Kuramoto-Sivashinsky Equation
In this section we formulate an LES system corresponding to the Kuramoto-Sivashinsky system (4) where the closure uses a Smagorinsky-type eddy-viscosity model of the general form (3) reduced to 1D. We must first define the filtering operation to extract the resolved scales from the solutions and for this purpose we shall use a sharp spectral filter also employed in [16]. It is defined in terms of the Fourier transform G δ (k) of its kernel (hats " · " will hereafter denote Fourier coefficients) where the cut-off length scale δ = 2π/k max and k max is the maximum resolved wavenumber. Clearly, (5) defines a low-pass filter which removes all Fourier components with wavenumbers larger than k max . Applying this filter to (4a), we obtain the filtered version of the Kuramoto-Sivashinsky equation where w := G δ ⋆ w and the last term represents the effect of the SGS stresses The reason for the sign difference with respect to the corresponding expression in (2a) is the sign of the dissipative term in (6). For clarity, we will hereafter use the symbol u to denote the solution of the LES problem for the Kuramoto-Sivashinsky equation, which should be contrasted with w obtained by filtering the solution w of the DNS problem (4). Since expression (7) depends on the original unresolved field w, it must be modelled in terms of w in order to close equation (6). For this purpose we will use a 1D analogue of the closure model (3) adapted to the Kuramoto-Sivashinsky equation, namely, in which ν(|s|) is the eddy viscosity. The ansatz in (8) is chosen such that the order of derivatives in the resulting model term will match the order (four) of the dissipative term in the Kuramoto-Sivashinsky system (4). The filtered Kuramoto-Sivashinsky system (6) equipped with such a closure model then becomes the LES system with the following form We now introduce two important definitions: where α := min x∈Ω, t∈[0,T ] |s| and β := max x∈Ω, t∈[0,T ] |s|, referred to as the "identifiability interval", is the range of values attained by the magnitude of the resolved strain |s| in the solution of the LES problem (9) with the initial data w 0 over the time interval [0, T ], , where a < α and b > β, will serve as the domain of definition of the function defining the eddy viscosity, i.e., ν : L → R; since the identifiability interval will in general depend on the initial data w 0 and the length of the time window T , i.e., I = I( w 0 , T ), it is important to choose the domain L such that it will contain all possible identifiability intervals, that is such that ∪ w 0 ,T I( w 0 , T ) ⊂ L, as this will ensure that the eddy viscosity is always defined; in practice, it is convenient to adopt a larger domain L possibly also including points outside any identifiability interval I( w 0 , T ), i.e., L \ ∪ w 0 ,T I( w 0 , T ) = ∅; with this in mind, we shall set a = inf x∈Ω, t∈[0,T ] |s| = 0 and b > sup x∈Ω, t∈[0,T ] |s|.
The counterpart of the Smagorinsky model in the present setting will then take the form Our goal now is to find an optimal form for the eddy viscosity as a function of the resolved strain s, ν = ν(|s|), generalizing the Smagorinsky model (10). This eddy viscosity will be optimal in the sense that the corresponding solutions of the LES system (9) will be as close as possible (in a suitable least-squares sense) to solutions of the original Kuramoto-Sivashinsky system (4) obtained for the same initial data w 0 . Formulation and solution of this optimization problem are presented below.

Optimization Approach to finding Eddy Viscosity
In this section we first state the optimization problem which will be used to determine the optimal form of the eddy viscosity. It is formulated here in a very general continuous setting and to solve this problem we use a gradient-descent approach in which the key element is a suitably-defined gradient representing the sensitivity of solutions to the LES system (9) to modifications of the functional form of the eddy viscosity. Finally, we ensure that these gradients are sufficiently smooth such that the resulting optimal eddy viscosity will be well defined. Starting from some initial guess ν 0 , the optimization procedure will iteratively adjust the eddy viscosity such that the corresponding solutions u of the LES problem (9) will match as closely as possible the "true" evolution governed by the original Kuramoto-Sivashinsky system (4), i.e., the DNS. To fix attention, this matching will be determined in terms of N "observations" made by applying suitably-defined observation operators . . , N, to the LES and DNS solutions, u(t, ·) and w(t, ·), continually for all t ∈ [0, T ]. The symbol H 1 (Ω) denotes the Sobolev space of continuous functions with square-integrable gradients [42]. We will use this general formulation to introduce our approach here and will define specific forms of these observation operators in Section 4 which will then be used in our computational examples in Section 6. The "target" observations will thus have the form m i (t) : We see that in order for the LES system (9) to be satisfied in the classical (strong) sense, the eddy viscosity ν(|s|) must possess certain minimum regularity as a function of s. Due to some technical reasons which will become apparent below, we must have Since optimization problems are most conveniently formulated in Hilbert spaces [43], we will assume the eddy viscosity ν(|s|) to be an element of the Sobolev space H 3 (L) of functions defined on L with square-integrable third derivatives [42] (a precise definition of the norm in this space will be provided below). The objective functional J : H 3 (L) → R will therefore have the form of the least-squares error between the target observations {m i (t)} N i=1 and the corresponding observations of solutions u of the LES problem (9) obtained for the given eddy viscosity ν, i.e., such that the optimization problem takes the form where q ν is the optimal eddy viscosity. To find a local minimizer of (12), we shall use a gradient-based optimization approach in which the optimal eddy viscosity q ν can be computed iteratively as q ν = lim n→∞ ν (n) , where in which ∇ ν J (ν) is the gradient of the cost functional (12) with respect to the eddy viscosity ν, τ (n) is the step length along the descent direction at the nth iteration, and ν 0 is the initial guess for the eddy viscosity. An optimal step size τ (n) can be found by solving the following line-minimization problem [44] We add that due to the local nature of this approach, iterations (14) can only produce local minimizers and determining whether any of them is also a global minimizer is in general not possible. A key element of the gradient-descent approach (14) is evaluation of the gradient ∇ ν J (ν (n) ) and this step is discussed below.

Adjoint-Based Gradients
While adjoint calculus has had a long history in PDE-constrained optimization starting with [45], which also includes problems in fluid mechanics [46], the optimization problem defined by (9), (12) and (13) has in fact a somewhat non-standard structure and requires therefore special techniques. The reason is that the control variable ν in (13) is a function of s, which itself is a function of the dependent variable u in (9), cf. (8), whereas standard adjoint-based methods allow one to solve PDE optimization problems in which the control is a function of the independent variables only (here, x and t). The approach used in this study was initially formulated in [32,33] and is suitably modified for our current purposes. Here we present its main elements only and the reader is referred to A for more technical details. We begin by computing the Gâteaux (directional) differential of the cost functional (12) with respect to ν where ν ′ is a perturbation of the eddy viscosity ν and u ′ (t, x; ν, ν ′ ) satisfies the corresponding linear perturbation system obtained from (9), cf. relations (35)-(36) in A. The (local) minimizer q ν of (13) requires the directional derivative (16) to vanish for all perturbations ν ′ , i.e., ∀ Away from the minimizer q ν we can use the Gâteaux differential to obtain the gradient ∇ ν J required by the descent algorithm (14). To do this, we invoke the Riesz representation theorem [47] and the fact that the directional derivative (16) is a bounded linear functional when viewed as a function of ν ′ , to obtain where ·, · X is an inner product in the Hilbert space X . As regards the choice of this space, we will consider X = L 2 (L) and X = H 3 (L) endowed with the following inner products where ℓ 1 , ℓ 2 and ℓ 3 in (18b) are "length-scale" parameters (we note that as long as 0 < ℓ 1 , ℓ 2 , ℓ 3 < ∞, the corresponding inner products are equivalent, in the precise sense of norm equivalence). While the Sobolev gradient obtained in the space H 3 (L) must be used in computations in (14), it is convenient to first derive the gradient defined with respect to the L 2 topology. We note that relation (16) is not consistent with the Riesz form (17), because the perturbation ν ′ does not appear in it explicitly as a factor, but is instead hidden in the perturbation equation (36a). As shown in A, relation (16) can be transformed to the desired Riesz form (17) using an adjoint state u * : [0, 2π] × [0, T ] → R, such that from this identity with X = L 2 (L) we can obtain the following expression for the L 2 gradient of the objective functional (12) where χ [a 1 ,a 2 ] (s) is the characteristic function of the interval [a 1 , a 2 ] ∈ I whereas u * is the solution of the adjoint system in which H * i are the adjoints of the observation operators H i , i = 1, . . . , N. We remark that the adjoint system (20) is a terminal-value problem and must be therefore integrated backwards in time whereas its coefficients are determined by the solution u(t, x) of the (forward) LES problem (9) around which linearization is performed (see A). When the adjoint system is properly defined, its solutions contain information about the sensitivity of the solutions to the LES problem (9), and hence also the error functional (12), to perturbations of the functional form of the eddy viscosity in (8). From the structure of the last term on the left-hand side (LHS) in (20a) it is also clear that in order for the adjoint system to be satisfied in the classical (strong) sense, the eddy viscosity must possess the minimum regularity specified in (11).
As defined in (19), the L 2 gradient may not be used in optimization algorithm (14) because it does not possess the required regularity, cf. (11), and is defined only on the identifiability interval I which in general is smaller than the domain of definition L of the eddy viscosity (this latter issue could in principle be remedied by extending the L 2 gradient (19) onto L\I with zeros). In order to get around these difficulties we will derive the corresponding Sobolev gradients [48,49] by setting X = H 3 (L) in the Riesz identity (17) which, upon noting (18b), leads to Sobolev gradients are determined subject to certain boundary conditions imposed on s = a and s = b [32,33] which in turn determine the behavior of the corresponding properties of the optimal eddy viscosity q ν(s) at s = a, b. There is some freedom as regards this choice and we shall use which implies that in the gradient iterations (14) the odd-degree derivatives of the eddy viscosity with respect to s will remain unchanged with respect to the initial guess ν 0 at s = a and the value of ν and its first two derivatives will remain unchanged at s = b. We remark that, in particular, the gradient iterations (14) are allowed to modify the value of q ν at s = a. It should be noted, however, that when the domain L is much larger than the identifiability interval I , i.e., a ≪ α and β ≪ b, the boundary conditions such as (22) have little effect on the optimal eddy viscosity q ν. Performing integration by parts with respect to s in (21) required number of times and noting that due to the judicious choice of the boundary conditions (22) all boundary terms vanish we obtain which because of the arbitrariness of the perturbation ν ′ ∈ H 3 (L) implies In (23) the L 2 gradient appearing on the RHS is extended from the identifiability interval I to the entire domain L with zeros. We note that extraction of Sobolev gradients by solving the boundary-value problem (22)-(23) can be interpreted as applying a low-pass filter to the L 2 gradient. Indeed, the Fourier transform of (23) yields  (23). Both axes have logarithmic scaling and it is assumed that ℓ 1 = 0 and ℓ 2 > ℓ 3 .
which shows that the cut-off for filtering is determined by the length scales ℓ 1 , ℓ 2 and ℓ 3 , cf.

Observation Operators
We start by defining the observation operators and consider two choices corresponding to observations in the physical and in the spectral (Fourier) space.

Physical-Space Observations
We will assume here that for all times t ∈ [0, T ] the solution u(t, x) of the LES system (9) is observed at a certain number of observation points {x i } N i=1 uniformly distributed over the spatial domain [0, 2π]. The observation operator H i : H 1 (0, 2π) → R associated with the ith point is then given by (we express it here in an integral form in order to facilitate obtaining its adjoint H * i ) The adjoint system (20) also involves the adjoints H * i : R → H −1 (0, 2π) of H i , i = 1, . . . , N, which can be obtained from the following duality relations where f ∈ R and ·, · R denotes the (scalar) product of two real numbers. From this we thus deduce

Fourier-Space Observations
Since the periodic Kuramoto-Sivashinsky system (4) is employed here as a "toy model" for homogeneous turbulence, another natural way to define the observation operators H i is in terms of the Fourier (e.g., cosine) transform of the state and this is the second choice we shall consider where K is the set of wavenumbers corresponding to the observed Fourier components (with cardinality |K| = N). The adjoints H * i of these observation operators are then obtained by considering the duality relations

Physical Parameters
The long-time behavior of the solutions w of the Kuramoto-Sivashinsky system (4) is determined by the parameters ν 4 and ν 2 [41]. In our study we shall use the values ν 4 = 1 and ν 2 = 100 chosen such that, after an initial transient, the solution w will on average feature 7 waves ("coherent structures") present in the domain during the evolution (this number coincides with the wavenumber k 0 of the most (linearly) unstable mode of the Kuramoto-Sivashinsky system (4) linearized about the zero state w(t, x) = 0 [41]). Given the chaotic nature of the Kuramoto-Sivashinsky system in this parameter regime and our interest in the long-time evolution, as the initial condition w 0 we will take a certain state on the turbulent attractor. When defining the LES system (9) we will take the cut-off wavenumber in filter (5) to be k max = 16, which falls in the "inertial range" not too far from the wavenumber k 0 = 7 characterizing the most unstable modes (which can be interpreted as "forcing"), cf. Figure 2. In the solution of the optimization problem (28b) The number of observations N is chosen to be smaller than the number k max of the Fourier components resolved in the LES (9), cf. Figure 2. As regards the initial guess ν 0 for the eddy viscosity in the gradient descent (14) we will take the Smagorinsky model

Numerical Approach
The gradient-descent approach (14) is formulated in the continuous ("optimize-thendiscretize") setting [46] and evaluation of the gradient expression (19) requires solutions of the LES and the adjoint systems (9) and (20). In this section we first discuss the numerical solution of these PDE problems and the computation of the Sobolev gradients via (22)- (23). Then we describe the implementation of the gradient-descent algorithm (14).

Discretization
The LES and adjoint systems (9) and (20) involve model terms with state-dependent eddy viscosity ν(|s|) and in order to represent this expression, in addition to discretizing the space and time domains [0, 2π] and [0, T ], we also need to discretize the state domain L. The former two domains are discretized using grids with equispaced points and steps sizes ∆x = 2π/N x and ∆t, where N x is the number of grid points in space, whereas the state domain L is discretized with N s Chebyshev points. The original Kuramoto-Sivashinsky system (4), its LES version (9) and the adjoint system (20) are solved using the standard Fourier pseudo-spectral method [50] where dealiasing based on the 3/2 rule is performed in the case of the Kuramoto-Sivashinsky system (4), but due to aggressive filtering, cf. (5), is unnecessary in the latter two problems. Evaluation of the model terms in (9) and (20) requires differentiation of the eddy viscosity ν(|s|) with respect to s which is performed using spectrally-accurate Chebyshev differentiation matrices [51] defined in the state domain L. The eddy viscosity ν(|s|) and its derivative dν(|s|)/ds are then interpolated from the state space to the physical space using the barycentric formulas which are also spectrally accurate [52]. The time-discretization of systems (4), (9) and (20) is performed using the exponential time-differencing fourth-order Runge-Kutta method (ETDRK4) [53], originally introduced in [54], which is fourth-order accurate. The different integrals are approximated using Gaussian quadratures, which are given by the trapezoidal rule on the periodic domain Ω (e.g., in (12)), and by the Clenshaw-Curtis formula on the bounded domain L (e.g., in (18)). The boundary-value problem (22)-(23) defining the Sobolev gradients is solved using the chebop feature of Chebfun [55], where the discretization is performed based on ultraspherical polynomials. With most computations carried out with spectral accuracy, approximation errors are dominated by time-stepping errors where the accuracy is O((∆t) 4 ). Unless mentioned otherwise, in our computations we use N x = 1024, N s = 4096 and ∆t = 3.0 × 10 −6 . In order to validate the discretization techniques discussed above, we verify the accuracy of the cost functional gradients evaluated as described in Section 3.1, cf. (19), by computing the Gâteaux differential J ′ (ν; ν ′ ) in terms of the Riesz identity (17) and comparing it with its approximation obtained with a simple forward finite-difference formula. The ratio of these two expressions is thus given by where ν ′ is an arbitrary perturbation and ǫ > 0 its magnitude. We expect κ(ǫ) to be close to unity and this is indeed evident in Figure 3 for a range of ǫ spanning several orders of magnitude. The large deviations of κ(ǫ) from unity observed for very small and very large values of ǫ are due to, respectively, the round-off and truncation errors in the finite-difference formula, both of which are well-known effects [32,33]. In Figure 3(b) we also note that, as expected, for intermediate values of ǫ, κ(ǫ) → 1 as the discretization parameter ∆t used in the numerical solution of the PDE systems (9) and (20) is refined. These results demonstrate the accuracy of the cost functional gradients evaluated as  (9) and (20). discussed in Section 3.1.

Gradient Descent
While for simplicity of presentation in Section 3 the steepest-descent (simple gradient) approach was used in (14), in actual computations we use the conjugate-gradients method [44] which is known to significantly accelerate convergence. The gradient ∇ ν J (ν (n) ) in (14) is then replaced with the descent direction g (n) defined as g (n) = ∇ ν J (ν (n) ) − γ n g (n−1) , n = 1, 2, . . . , where the "momentum" term γ n is evaluated using the Polak-Ribière formula It is a good practice [44] for the conjugate-gradients approach (30) to be periodically restarted with a gradient step after a certain number of iterations. The line-minimization problem (15) is efficiently solved using Brent's algorithm [56], which is a standard approach. Gradient iterations (14) are declared converged when the following termination condition is satisfied Algorithm 1 Implementation of the conjugate-gradients variant of descent approach (14). Input: -target observations (e.g., of the DNS, cf. (4); N is number of observations) N x , N s , ∆t -numerical discretization parameters ℓ 1 , ℓ 2 , ℓ 3 -Sobolev length scales, cf. (18b) ǫ J -tolerance in the termination criterion (32) ν 0 -initial guess for the eddy viscosity (e.g., (10)) Output: q ν -optimal eddy viscosity • determine the conjugate descent direction g (n) via (30) • determine the optimal step length τ (n) by solving (15) with Brent's algorithm • update the eddy viscosity ν (n) ← ν (n) − τ n g (n) until termination criterion (32) is satisfied where ǫ J is a prescribed tolerance (we will use ǫ J = 10 −7 ). The choice of the length-scale parameters ℓ 1 , ℓ 2 and ℓ 3 defining the Sobolev inner product (18b) will be discussed in the next section. The different steps in the solution of the optimization problem (13) are summarized as Algorithm 1.

Results
In this section we present and analyze the results obtained with our approach to determining optimal closure models, cf. Algorithm 1, and compare them to the results obtained with other closure models including the approach of Das & Moser [16], where the closure model also has some optimality properties. We consider observation operators H i , i = 1, . . . , N, defined both in the physical space, cf. Section 4.1.1, and in the Fourier space, cf. Section 4.1.2, the latter with different distributions of the observed wavenumbers (28). With most other problem parameters fixed as discussed in Section 4.2, there remains one key characteristic defining the optimization problem (12) Table 1: Summary information about the different cases considered in our computations.
namely, the length of the time window and different values of T are considered, covering from a few to several typical events in the evolution of the Kuramoto-Sivashinsky system (4) (these "events" are the merging of crests or formation on new ones). Summary information about these computations is compiled in Table 1, where we also indicate the values of the length-scale parameters ℓ 2 and ℓ 3 (ℓ 1 = 0) used to determine the Sobolev gradients, cf. (18b). The values of these parameters are chosen by trial and error to maximize the convergence of iterations in Algorithm 1. We conclude from the data compiled in Table 1 that in all cases optimization reduces the observation error (12) by a factor O(1) with its specific value depending on the length of the optimization window T and there tends to be an optimal value of T for which the largest reduction of the error functional (12) is achieved. To fix attention, we will henceforth focus on two representative configurations, namely, one with observations in the physical space and one with observations in the Fourier space, denoted respectively "Case A" and "Case B" in Table 1.
The decrease of the normalized objective functional (12) with iterations n in Cases A and B is shown in Figure 4(a). In this figure we observe a reduction of the observation error by close to one order of magnitude over O(10) iterations. We note however that the convergence rate of iterations (14) is rather nonuniform. The corresponding optimal eddy viscosities q ν(s) are shown as functions of the resolved strain s together with the Smagorinsky model (10) used as the initial guess in (14) in Figure 4(b). We see that while the optimal eddy viscosities q ν(s) are defined on a larger domain L, the deviations from the initial guess ν 0 produced by the gradient iterations (14) are essentially confined to a smaller identifiability interval I. Most importantly, in contrast to the original Smagorinsky model (10), the optimal eddy viscosities are negative for small strains such that q ν(0) < 0. It is encouraging to note that the optimal eddy viscosities obtained in Cases A and B exhibit a qualitatively similar dependence on s, despite rather different forms of observations used to define the optimization problem (12)-(13) in these two cases.
In order to obtain insights about the spatio-temporal evolution of solutions to the LES problems (9) with different closure models (no closure at all, the Smagorinsky model (10), and the optimal eddy viscosity q ν(|s|) from Cases A and B, cf. Table 1 To quantify these observations, we will now analyze the behavior of the following three diagnostic quantities which can be interpreted as, respectively, the correlation of the LES solution u with the DNS solution w and its the normalized kinetic energy and enstrophy (in (33c) the symbol | · | H 1 denotes the H 1 seminorm). The LES with the optimal closures obtained in Cases A and B are compared in terms these diagnostic quantities for t ∈ [0, 2T ] to the LES with no closure model, with the Smagorinsky model (10), and with the closure   (10), with (blue dash-dotted line) the optimal eddy viscosity q ν from case A, with (red dashed line) the optimal eddy viscosity q ν from case B, and with (yellow solid-dotted line) the optimal closure proposed by Das & Moser [16]. Bold and thin lines correspond to time in the "training window" t ∈ [0, T ] and beyond the "training window" t ∈ (T, 2T ], respectively. model of Das & Moser [16] in Figure 6. Given the chaotic nature of the Kuramoto-Sivashinsky system (4) resulting in an exponentially fast divergence of initially nearby trajectories, in all cases the correlation (33a) drops very rapidly, such that for short times t → 0 we approximately have C(t) ≈ 1 − r e pt for some r, p > 0 different in each case, cf. Figures 6(a). The growth rate p, characterizing the exponential divergence of the solutions to the DNS and LES problems, is smallest when the optimal closure models q ν(|s|) from cases A and B are used. As a result, the time t 0 when the DNS and the LES solutions become uncorrelated, i.e., when C(t 0 ) ≈ 0, is nearly twice as large for the optimal closures from Cases A and B than for the Smagorinsky model (10) and the closure model of Das & Moser [16]. As regards the behavior of the normalized energy and enstrophy (33b)-(33c), from Figures 6(c,e) we see that the optimal eddy viscosity q ν(|s|) on average tends to reduce these quantities below their levels in the DNS. This is in contrast to the approach of Das & Moser [16] in which the normalized energy and enstrophy are increased to levels higher than in the DNS, cf. Figures 6(d,f).
While the analysis above focused on "a posteriori" tests involving the results of solving the LES system (9) with different closure models, we close this section with a brief discussion of an "a priori" test where the errors in approximations of the SGS stresses (7) with different closure models are analyzed. In this context, comparing the SGS dissipation rate with the modeled SGS dissipation rate is often used as a standard diagnostic to assess the energetics in the LES flow [57]. We thus focus on the following normalized least-squares measure of this error and show its dependence on time t ∈ [0, 2T ] for different models in Figure 7. We note in this figure that the normalized SGS stress error S(t) tends to be quite large in all cases, which is a known property of Smagorinsky-type models [58], and is smallest in the case of Das' & Moser's approach [16].

Discussion and Conclusions
In this study we have introduced a computational framework for determining optimal closures for a broad class of PDE models, such as the Navier-Stokes equation (1) in hydrodynamics. This inverse problem is framed here as a PDE optimization problem where an error functional representing the misfit between the target and predicted observations is minimized with respect to the functional form of the closure relation which determines the "shape" of the nonlinearity of the model term. Because of this latter aspect, such a problem is not amenable to solution using standard adjoint-based tools for PDE optimization and recently developed extensions of these techniques [32,33] need to be used. In addition, by formulating the problem in the "optimize-then-discretize" setting, rather than with a "discretize-then-optimize" approach often used to solve PDE optimization problems [46], we are able to determine the optimal closure relations in a very general  (10), with (blue dash-dotted line) the optimal eddy viscosity q ν from case A, with (red dashed line) the optimal eddy viscosity q ν from case B, and with (yellow solid-dotted line) the optimal closure proposed by Das & Moser [16]. Bold and thin lines correspond to time in the "training window" t ∈ [0, T ] and beyond the "training window" t ∈ (T, 2T ], respectively. form subject only to some minimal assumptions on smoothness, cf. (11), and the behavior for small and large values of the state variable. This ought to be contrasted with some earlier approaches in which optimal closures were determined by fitting a small number of parameters in an assumed ansatz. Thus, our proposed approach does not suffer from the limitations of such an assumed ansatz.
To provide a proof of the concept for this approach, we have focused here on the 1D Kuramoto-Sivashinsky equation (4) as a model problem computationally more tractable than the 2D or 3D Navier-Stokes system (1) we are ultimately interested in. In being chaotic and multiscale [41], the solutions of the Kuramoto-Sivashinsky system arguably better resemble actual hydrodynamic turbulence than solutions of the Burgers equation often used as a simplified model in similar situations [16]. Such simplified setting allows us to study the properties of our proposed approach more thoroughly. We find that it is possible to determine a particular dependence of the eddy viscosity on the resolved strain, q ν = q ν(|s|), generalizing Smagorinsky's relation (10) in a LES model, such that the LES model provides a systematically improved approximation of the DNS. More precisely, while the trajectories corresponding to the DNS and LES still diverge exponentially, they do so at a much slower rate than in the case of the LES based on the standard Smagorinsky model, cf. Figure 6(a). Importantly, the functional form of the optimal eddy viscosity q ν(|s|) appears to little depend on the particular form of the observations used to set up the optimization problem, cf. Figure 4(b), and we emphasize here that our approach can be formulated based on very general measurements (changing the observation operator H i , cf. Section 4.1, will only result in a modification of the source term in the adjoint system (20a)). In particular, it is also possible to simultaneously use several different sets of measurements coming, for example, from different DNS or experiments, and a natural way to formulate such a problem is in terms of multiobjective optimization.
The LES with the optimal closure models determined with the proposed approach produces solutions which better match the DNS than the LES solutions obtained with the closure model of Das & Moser [16], cf. Figure 6(b). On the other hand, that latter model leads to a more accurate prediction of the SGS stresses, cf. Figure 7. This can be understood by recognizing that the closure model of Das & Moser is formulated to optimally reconstruct the SGS stresses rather than some other a posteriori quantities. We have also considered an optimization problem with observation operators involving SGS stresses, but the results obtained were inferior to the results presented here.
In addition to providing optimal closure models which can be quite useful in many situations, the present approach also serves another, more fundamental purpose, namely, by identifying the "best" closure models within a certain family it can offer information about the fundamental performance limitations of such class of models. More specifically, it can provide insights about how well closure models based on the eddy-viscosity ansatz (3) can perform in the best case and thus how much room there is in principle for improvement of standard approaches such as the Smagorinsky model (10).
Moving forward, our next main goal is to consider an analogous problem of finding optimal closure models for the 2D and 3D Navier-Stokes system (1) as generalizations of the Smagorinsky model (3), first in the periodic setting and then in more realistic geometries. In that latter context related questions also arise as regards wall models. We are also interested in studying closure models based on formulations other than eddy viscosity.
representations in which 0 < ǫ ≪ 1 The second term on the RHS in (35b) reflects the change of the value of the (squared) strain σ for which the eddy viscosity is evaluated as a result of perturbing its functional form [32,33]. Substituting representations (35a)-(35b) into the LES system (9) and collecting terms of O(ǫ) we obtain the following perturbation system u ′ (0, x) = 0 (36c) describing the leading-order effect u ′ of perturbing the functional form of the eddy viscosity ν(σ) on solutions of the LES system (9) [32,33]. Now we integrate (36a) against the adjoint field u * over the space-time domain [0, 2π] × [0, T ] and then perform integration by parts with respect to both space and time to obtain where all the boundary terms vanish due to periodicity. Using the definition of the adjoint system (20) together with the aforementioned change of variables σ = σ(s) we then obtain for the Gâteaux differential J ′ (ν; ν ′ ) =´T 0´2 π 0 ∂ u * ∂x ∂ 3 u ∂x 3 ν ′ dx dt, which now contains the perturbation ν ′ as a factor, but is still not in the Riesz form (17) because this form involves an inner product defined with integration with respect to the resolved strain s rather than space x and time t (we shall now return back to the original variable via the substitution σ = s 2 ). The required change of variables is introduced by the representation [32,33] ν ′ ∂ u(t, x) ∂x =ˆb a δ ∂ u(t, x) ∂x − s ν ′ (s) ds, where δ(·) is the Dirac delta distribution, such that using Fubini's theorem to swap the order of integration we finally arrive at a Riesz representation of the Gâteaux differential (16) from which after selecting X = L 2 (L) in (17) we deduce the following expression for the L 2 gradient We note that evaluation of this expression for a given value of s requires computation of an integral defined on the level set ∂ u(t,x) ∂x = s in the space-time domain [0, 2π] × [0, T ] which is rather difficult. A computationally more convenient approach is obtained using the following identity (in which the differentiation is understood in the distributional sense) such that (39) becomes and expression (19) is finally obtained. We emphasize that the need for transformation (37) and the subsequent steps is a consequence of the non-standard form of the PDE optimization problem considered here.