Super-resolving sparse observations in partial differential equations: A
  physics-constrained convolutional neural network approach
We propose the physics-constrained convolutional neural network (PC-CNN) to
infer the high-resolution solution from sparse observations of spatiotemporal
and nonlinear partial differential equations. Results are shown for a chaotic
and turbulent fluid motion, whose solution is high-dimensional, and has fine
spatiotemporal scales. We show that, by constraining prior physical knowledge
in the CNN, we can infer the unresolved physical dynamics without using the
high-resolution dataset in the training. This opens opportunities for
super-resolution of experimental data and low-resolution simulations.
Daniel Kelshaw, Luca Magri
# Super-resolving sparse observations in partial differential equations:

###### Abstract

We propose the physics-constrained convolutional neural network (PC-CNN) to infer the high-resolution solution from sparse observations of spatiotemporal and nonlinear partial differential equations. Results are shown for a chaotic and turbulent fluid motion, whose solution is high-dimensional, and has fine spatiotemporal scales. We show that, by constraining prior physical knowledge in the CNN, we can infer the unresolved physical dynamics without using the high-resolution dataset in the training. This opens opportunities for super-resolution of experimental data and low-resolution simulations.

## 1 Introduction

Observations of turbulent flows and physical systems are limited in many cases, with only sparse or partial measurements being accessible. Access to limited information obscures the underlying dynamics and provides a challenge for system identification [e.g., 1]. Super-resolution methods offer the means for high-resolution state reconstruction from limited observations, which is a key objective for experimentalists and computational scientists alike. Current methods for super-resolution or image reconstruction primarily make use of convolutional neural networks due to their ability to exploit spatial correlations. In the classic, data-driven approach there is a requirement for access to pairs of low- and high-resolution samples, which are needed to produce a parametric mapping for the super-resolution task [e.g., 2; 3; 4; 5]. Considering observations of physical systems, the problems encountered in the absence of ground-truth high-resolution observations can be mitigated by employing a physics-constrained approach by imposing prior knowledge of the governing equations [e.g., 6; 7; 8].

Physics-informed neural networks (PINNs) [7] have provided a tool for physically-motivated problems, which exploits automatic-differentiation to constrain the governing equations. Although the use of PINNs for super-resolution shows promising results for simple systems [9], they remain challenging to train, and are not designed to exploit spatial correlations [10; 11]. On the other hand, convolutional neural networks are designed to exploit spatial correlations, but they cannot naturally leverage automatic differentiation to evaluate the physical loss, as PINNs do, because they provide a mapping between states, rather than a mapping from the spatiotemporal coordinates to states as in PINNs. As such, the design of physics-informed convolutional neural networks is more challenging, and rely on finite-different approximations or differentiable solvers [12; 13]. For example, the authors of [13] show results for a steady flow field (with no temporal dynamics), which produces a mapping for stationary solutions of the Navier-Stokes equations. In this work, we design a framework to tackle spatiotemporal partial differential equations. For this, we propose a super-resolution method that does not require the full set of high-resolution observations. To accomplish this, we design and propose the physics-constrained convolutional neural network (PC-CNN) with a time windowing scheme. We apply this to a two-dimensional chaotic and turbulent fluid motion. The paper is structured as follows. In SS2 we introduce the super-resolution task in the context of a mapping from low- to high- resolution. The choice of architecture used to represent
the mapping is discussed in SS3 before providing an overview of the methodology SS4. We introduce the two-dimensional, turbulent flow in SS5 and provide details on the pseudospectral method used for discretisation. Finally, results are demonstrated in SS6. Conclusions SS7 end the paper.

## 2 Super-resolution task

We consider partial differential equations (PDEs) of the form

\[\mathcal{R}(\mathbf{\tilde{u}};\lambda)\equiv\partial_{t}\mathbf{\tilde{u}}-\mathcal{N} (\mathbf{\tilde{u}};\lambda), \tag{1}\]

where \(\mathbf{x}\in\Omega\subset\mathbb{R}^{n}\) denotes the spatial location; \(t\in[0,T]\subset\mathbb{R}_{\geq 0}\) is the time, where \(n\) is the space dimension; \(\mathbf{\tilde{u}}:\Omega\times[0,T]\rightarrow\mathbb{R}^{m}\) is the state, where \(m\) is the number of components of the vector function \(\mathbf{\tilde{u}}\); \(\lambda\) are the physical parameters of the system; \(\mathcal{N}\) is a sufficiently smooth differential operator, which represents the nonlinear part of the system of \(m\) partial differential equations; \(\mathcal{R}\) is the residual; and \(\partial_{t}\) is the partial derivative with respect to time. A solution \(\mathbf{u}\) of the PDE (1) is the function that makes the residual vanish, i.e., \(\mathcal{R}(\mathbf{u};\lambda)=0\).

Given sparse observations, \(\mathbf{u}(\mathbf{\Omega}_{L},t)\), we aim to reconstruct the underlying solution to the partial differential equation on a high-resolution grid, \(\mathbf{u}(\mathbf{\Omega}_{H},t)\), where the domain \(\Omega\) is discretised on low- and high-resolution uniform grids, \(\mathbf{\Omega}_{L}\subset\mathbb{R}^{N^{n}}\) and \(\mathbf{\Omega}_{H}\subset\mathbb{R}^{M^{n}}\), respectively, such that \(\mathbf{\Omega}_{L}\cap\mathbf{\Omega}_{H}=\mathbf{\Omega}_{L}\); \(M=\kappa N\); and \(\kappa\in\mathbb{N}^{+}\) is the up-sampling factor. Our objective is to find a parameterised mapping \(\mathbf{f_{\theta}}\) such that

\[\mathbf{f_{\theta}}:\mathbf{u}(\mathbf{\Omega}_{L},t)\rightarrow\mathbf{u}(\mathbf{\Omega}_{H},t). \tag{2}\]

We consider the solution, \(\mathbf{u}\), to be discretised with \(N_{t}\) times steps, \(t_{i}\), which are contained in the set \(\mathcal{T}=\left\{t_{i}\in[0,T]\right\}_{i=0}^{N_{t}}\). We approximate the mapping \(\mathbf{f_{\theta}}\) by a convolutional neural network parameterised by \(\mathbf{\theta}\).

## 3 Convolutional neural networks

Parameterising \(\mathbf{f_{\theta}}\) as a convolutional neural network is a design choice, which allows us to exploit two key properties. First, shift invariance, which is consequence of kernel-based operations acting on structured grids [14]. Second, partial differential equations are defined by local operators, which is a property that we wish the machine to induce naturally. Therefore, we choose convolutional neural networks, which employ an architectural paradigm that leverages spatial correlations [15], through the composition of functions

\[\mathbf{f_{\theta}}=\mathbf{f_{\theta_{Q}}^{Q}}\circ\cdots\circ h(\mathbf{f_{\theta_{1}}^ {1}})\circ h(\mathbf{f_{\theta_{0}}^{0}}), \tag{3}\]

where \(\mathbf{f_{\theta_{i}}^{i}}\) denote discrete convolutional layers; \(h\) is an element-wise nonlinear activation function, which increases the expressive power of the network [e.g., 16] and yields a universal function approximator [17]; and \(Q\) is the number of layers. The convolutional layers are responsible for the discrete

Figure 1: Super-resolution task. The model \(\mathbf{f_{\theta}}\) is responsible for mapping the low-resolution field \(\mathbf{u}(\mathbf{\Omega}_{L},t)\) to the high-resolution field \(\mathbf{u}(\mathbf{\Omega}_{H},t)\). The upsampling layers (middle white layers in the leftmost set of data) perform bi-cubic upsampling to obtain the correct spatial dimensions. Convolutional layers are parameterised as Conv2D(c_in, c_out, k), where c_in, c_out denote the number of input and output channels, respectively; and k is the spatial extent of the filter. The terms \(\mathcal{L_{O}},\mathcal{L_{R}}\) denote the observation-based loss and residual loss, respectively, the combination of which forms the objective loss \(\mathcal{L_{\theta}}\), which is used to update the network’s parameters, \(\mathbf{\theta}\). We provide an in-depth explanation of these losses in §4.
operation \((\mathbf{x},\mathbf{w},\mathbf{b})\mapsto\mathbf{w}*\mathbf{x}+\mathbf{b}\), where \(\mathbf{\theta}=(\mathbf{w},\mathbf{b})\). As the kernel operates locally around each pixel, information is leveraged from the surrounding grid cells. This makes convolutional layers an excellent choice for learning and exploiting spatial correlations, as are often important in the solutions to partial differential equations (e.g., 18; 13).

We propose a physics-constrained convolutional neural network (PC-CNN), which consists of three successive convolutional layers that are prepended by an upsampling operation tasked with increasing the spatial resolution of the input. Knowledge of the boundary conditions can be imposed through the use of padding; for instance periodic padding for periodic boundary conditions. Figure 1 provides an overview of the super-resolution task with information about the architecture employed.

## 4 Methodology

Having defined the mapping \(\mathbf{f_{\theta}}\) as a convolutional neural network parameterised by \(\mathbf{\theta}\), we formalise an optimisation problem to minimise the cost function \(\mathcal{L}_{\mathbf{\theta}}\)

\[\mathbf{\theta}^{*}=\operatorname*{argmin}_{\mathbf{\theta}}\mathcal{L}_{\mathbf{\theta}} \quad\text{where}\quad\mathcal{L}_{\mathbf{\theta}}=\mathcal{L}_{\mathcal{O}}+ \alpha\mathcal{L}_{\mathcal{R}}, \tag{4}\]

where \(\alpha\) is a non-negative empirical regularisation factor, which determines the relative importance of the corresponding loss terms. Given low-resolution observations \(\mathbf{u}(\mathbf{\Omega}_{L},t)\) at arbitrary times \(t\in\mathcal{T}\), we define each of the loss terms as

\[\mathcal{L}_{\mathcal{O}} =\frac{1}{N_{t}}\sum_{t\in\mathcal{T}}\lVert\mathbf{f_{\theta}}(\mathbf{ u}(\mathbf{\Omega}_{L},t))\big{|}^{\mathbf{\Omega}_{L}}-\mathbf{u}(\mathbf{\Omega}_{L},t) \rVert^{2}_{\mathbf{\Omega}_{L}}, \tag{5}\] \[\mathcal{L}_{\mathcal{R}} =\frac{1}{N_{t}}\sum_{t\in\mathcal{T}}\lVert\mathcal{R}(\mathbf{f_{ \theta}}(\mathbf{u}(\mathbf{\Omega}_{L},t));\lambda)\rVert^{2}_{\mathbf{\Omega}_{H}},\]

where \(\mathbf{f_{\theta}}(\,\cdot\,)\big{|}^{\mathbf{\Omega}_{L}}\) denotes the corestriction of \(\mathbf{\Omega}_{H}\) on \(\mathbf{\Omega}_{L}\), and \(\lVert\,\cdot\,\rVert_{\mathbf{\Omega}}\) represents the \(\ell^{2}\)-norm over the domain \(\mathbf{\Omega}\). In order to find an optimal set of parameters \(\mathbf{\theta}^{*}\), the loss is designed to regularise predictions that do not conform to the desired output. Given sensor observations \(\mathbf{u}(\mathbf{\Omega}_{L},t)\) in the low-resolution input, the observation-based loss, \(\mathcal{L}_{\mathcal{O}}\), is defined to minimise the distance between known observations, \(\mathbf{u}(\mathbf{\Omega}_{L},t)\), and their corresponding predictions on the high-resolution grid, \(\mathbf{f_{\theta}}(\mathbf{u}(\mathbf{\Omega}_{L},t))\).

We impose the prior knowledge that we have on the dynamical system by defining the residual loss, \(\mathcal{L}_{\mathcal{R}}\), which penalises the parameters that yield predictions that violate the governing equations (1)1. Mathematically, this means that weIn absence of observations (i.e., data), the residual loss \(\mathcal{L}_{\mathcal{R}}\) alone does not provide a unique solution. By augmenting the residual loss with the data loss, \(\mathcal{L}_{\mathcal{O}}\) (see Eq. (4)), we ensure that network realisations conform to the observations (data) whilst fulfilling the governing equations, e.g., conservation laws. Crucially, as consequence of the proposed training objective, we do not need the high-resolution field as labelled dataset, which is required with conventional super-resolution methods.

Footnote 1: For example, in turbulence, the fluid mass and momentum must be in balance with mass sources, and forces, respectively. Thus, (1) represent conservation laws.

### Time-windowing of the residual loss

Computing the residual of a partial differential equation is a temporal task, as shown in Eq. (1). We employ a time-windowing approach to allow the network \(\mathbf{f_{\theta}}\) to learn the sequentiality of the data. This provides a means to compute the time-derivative \(\partial_{t}\mathbf{u}\) required for the residual loss \(\mathcal{L}_{\mathcal{R}}\). The network takes time-windowed samples as inputs, each sample consisting of \(\tau\) sequential time-steps. The time-derivative \(\partial_{t}\mathbf{u}\) is computed by applying a forward Euler approximation to the loss

\[\mathcal{L}_{\mathcal{R}} =\frac{1}{\tau N_{t}}\sum_{t\in\mathcal{T}}\sum_{n=0}^{\tau} \big{\lVert}\partial_{t}\mathbf{f_{\theta}}(\mathbf{u}(\mathbf{\Omega}_{L},t+n\Delta t)) \tag{6}\] \[\qquad-\mathcal{N}(\mathbf{f_{\theta}}(\mathbf{u}(\mathbf{\Omega}_{L},t+n \Delta t));\lambda)\big{\rVert}^{2}_{\mathbf{\Omega}_{H}}.\]

Using this approach, we are able to obtain the residual for the predictions in a temporally local sense; computing derivatives across discrete time windows rather than the entire simulation domain. The network \(\mathbf{f_{\theta}}\) is augmented to operate on these time windows, which vectorises the operations over the time-window. For a fixed quantity of training data, the choice of \(\tau\) introduces a trade-off between the number of input samples \(N_{t}\), and the size of each time-window \(\tau\). We take a value \(\tau=2\), which is the minimum window size for computing the residual loss \(\mathcal{L}_{\mathcal{R}}\). We find that this is sufficient for training the network whilst simultaneously maximising the number of independent samples used for training. To avoid duplication of the data in the training set, we ensure that all samples are at least \(\tau\) time-steps apart so that independent time-windows are guaranteed to not contain any overlaps.
## 5 Chaotic and turbulent dataset

As nonlinear partial differential equations, we consider the Navier-Stokes equations, which are the expressions of the conservation of mass and momentum of fluid motion, respectively

\[\begin{split}\nabla\cdot\mathbf{u}&=0,\\ \partial_{t}\mathbf{u}+\left(\mathbf{u}\cdot\mathbf{\nabla}\right)\mathbf{u}& =-\nabla p+\nu\Delta\mathbf{u}+\mathbf{g},\end{split} \tag{7}\]

where \(p,\nu\) denote the scalar pressure field and kinematic viscosity, respectively. The flow velocity, \(\mathbf{u}\in\mathbb{R}^{m=2}\) evolves on the domain \(\Omega\in[0,2\pi)\subset\mathbb{R}^{n=2}\) with periodic boundary conditions applied on \(\partial\Omega\), and a stationary, spatially-varying sinusoidal forcing term \(\mathbf{g}(\mathbf{x})\). In this paper, we take \(\nu=\nicefrac{{1}}{{42}}\) to ensure chaotic and turbulent dynamics, and employ the forcing term \(\mathbf{g}(\mathbf{x})=\left[\sin\left(4\mathbf{x_{2}}\right),0\right]^{\top}\)[19], where \(\mathbf{x_{2}}\) is the transverse coordinate. This flow, which is also known as the Kolmogorov flow [19], provides a nonlinear and multi-scale dateset, which allows us to evaluate the quality of predictions across the turbulent spectrum.

### Differentiable pseudospectral discretisation

To produce a solution for the Kolmogorov flow, we utilise a differentiable pseudospectral spatial discretisation \(\hat{\mathbf{u}}_{k}=\hat{\mathbf{u}}(\mathbf{k},t)\) where \(\hat{\mathbf{u}}=\mathcal{F}\circ\mathbf{u}\); \(\mathcal{F}\) is the Fourier transform; and \(\mathbf{k}\in\hat{\mathbf{\Omega}}_{k}\subset\mathbb{C}^{K^{n}}\) is the spectral discretisation of the spatial domain \(\Omega\in[0,2\pi)\)[20]. Operating in the Fourier domain eliminates the continuity term [21]. The equations for the spectral representation of the Kolmogorov flow are

\[\mathcal{R}(\hat{\mathbf{u}}_{k};\lambda)=\left(\tfrac{d}{dt}+\nu|\mathbf{k}|^{2} \right)\hat{\mathbf{u}}_{k}-\hat{\mathbf{f}}_{k}+\mathbf{k}\frac{\mathbf{k}\cdot\hat{\mathbf{f}}_{ k}}{|\mathbf{k}|^{2}}-\hat{\mathbf{g}}_{k}, \tag{8}\]

with \(\hat{\mathbf{f}}_{k}=-\left(\widehat{\mathbf{u}\cdot\nabla\mathbf{u}}\right)_{k}\), where nonlinear terms are handled pseudospectrally, employing the \(\nicefrac{{2}}{{3}}\) delaliasing rule to avoid unphysical culmination of energy at the high frequencies [21]. A solution is produced by time-integration of the dynamical system with the explicit forward-Euler scheme, choosing a time-step \(\Delta t\) that satisfies the Courant-Friedrichs-Lewy (CFL) condition. Initial conditions are generated by producing a random field scaled by the wavenumber, which retains the spatial structures of varying lengthscale in the physical domain [22]. The initial transient, which is approximately \(T_{t}=180s\), is removed from the dataset to ensure that the results are statistically stationary. (The transient time is case dependant.) For a spatial resolution \(\mathbf{\Omega}_{H}\in\mathbb{R}^{70\times 70}\), we use a spectral discretisation \(\hat{\mathbf{\Omega}}_{k}\in\mathbb{C}^{35\times 35}\) to avoid aliasing and resolve the smallest lengthscales possible.

### Residual loss in the Fourier domain

The pseudospectral discretisation provides an efficient means to compute the differential operator \(\mathcal{N}\), which allows us to evaluate the residual loss \(\mathcal{L}_{\mathcal{R}}\) in the Fourier domain

\[\begin{split}\mathcal{L}_{\mathcal{R}}=\frac{1}{\tau N_{t}}\sum_ {t\in\mathcal{T}}\sum_{n=0}^{\tau}&\|\partial_{t}\hat{\mathbf{f}}_{ \mathbf{\theta}}(\mathbf{u}(\mathbf{\Omega}_{L},t+n\Delta t))\\ -\hat{\mathcal{N}}(\hat{\mathbf{f}}_{\mathbf{\theta}}(\mathbf{u}(\mathbf{\Omega }_{L},t+n\Delta t)))\|^{2}_{\hat{\mathbf{\Omega}}_{k}},\end{split} \tag{9}\]

where \(\hat{\mathbf{f}}_{\mathbf{\theta}}=\mathcal{F}\circ\mathbf{f}_{\mathbf{\theta}}\), and \(\hat{\mathcal{N}}\) denotes the Fourier-transformed differential operator. The pseudospectral discretisation is fully differentiable, which allows us to numerically compute gradients with respect to the parameters of the network \(\mathbf{f}_{\mathbf{\theta}}\). Computing the loss \(\mathcal{L}_{\mathcal{R}}\) in the Fourier domain provides two advantages: _(i)_ periodic boundary conditions are naturally enforced, which enforces the prior knowledge in the loss calculations; and _(ii)_ gradient calculations yield spectral accuracy. In contrast, a conventional finite differences approach requires a computational stencil, the spatial extent of which places an error bound on the gradient computation. This error bound is a function of the spatial resolution of the field.

## 6 Results

First, we discuss the generation of the low-resolution data for the super-resolution task. Next, we show the ability of the PC-CNN to infer the high-resolution solution of the partial differential equation on points that are not present in the training set.

### Obtaining the low-resolution data

A high-resolution solution of the partial differential equation is generated on the grid \(\mathbf{\Omega}_{H}\) prior to extracting a low-resolution grid \(\mathbf{\Omega}_{L}\) with the downsampling factor of \(\kappa=\nicefrac{{N}}{{M}}\). Both the solver and residual loss are discretised with \(K=\nicefrac{{N}}{{2}}\) wavenumbers in the Fourier domain, which complies with the Nyquist-Shannon sampling criterion. The downsampling by \(\kappa\) is performed by extracting the value of the solution at spatial locations \(\mathbf{\Omega}_{L}\cap\mathbf{\Omega}_{H}\), i.e., a low-resolution representation of the high-resolution solution

\[\mathbf{u}(\mathbf{\Omega}_{L},t)\triangleq\mathbf{u}(\mathbf{\Omega}_{H},t)\big{|}^{\mathbf{ \Omega}_{L}}. \tag{10}\]
(In contrast, the use of a pooling method for downsampling would distort values in the low-resolution representation, which effectively modifies the high-resolution solution.)

### Comparison with standard upsampling

We showcase the results for a downsampling factor \(\kappa=7\). Results are compared with interpolating upsampling methods, i.e., bi-linear, and bi-cubic interpolation to demonstrate the ability of the method. We provide a notion of quantitative accuracy by computing the relative \(\ell^{2}\)-error between the true solution, \(\mathbf{u}(\mathbf{\Omega}_{H},t)\), and the corresponding network realisation, \(\mathbf{f}_{\mathbf{\theta}}(\mathbf{u}(\mathbf{\Omega}_{L},t))\)

\[e=\sqrt{\frac{\sum_{t\in\mathcal{T}}\lVert\mathbf{u}(\mathbf{\Omega}_{H},t)-\mathbf{f}_{ \mathbf{\theta}}(\mathbf{u}(\mathbf{\Omega}_{L},t))\rVert_{\mathbf{\Omega}_{H}}^{2}}{\sum_{t \in\mathcal{T}}\lVert\mathbf{u}(\mathbf{\Omega}_{H},t)\rVert_{\mathbf{\Omega}_{H}}^{2}}}. \tag{11}\]

Upon discarding the transient, a solution \(\mathbf{u}(\mathbf{\Omega}_{H},t)\in\mathbb{R}^{70\times 70}\) is generated by time-integration over \(12\times 10^{3}\) time-steps, with \(\Delta t=1\times 10^{-3}\). We extract the low-resolution solution \(\mathbf{u}(\mathbf{\Omega}_{L},t)\in\mathbb{R}^{10\times 10}\) as a candidate for the super-resolution task. We extract 2048 samples at random from the time-domain of the solution, each sample consisting of \(\tau=2\) consecutive time-steps. The _adam_ optimiser [23] is employed for training with a learning rate of \(3\times 10^{-4}\). We take \(\alpha=10^{3}\) as the regularisation factor for the loss \(\mathcal{L}_{\mathbf{\theta}}\) and train for a total of \(10^{3}\) epochs, which is empirically determined to provide sufficient convergence.

Figure 2 shows a snapshot of results for the streamwise component of the velocity field, comparing network realisations \(\mathbf{f}_{\mathbf{\theta}}(\mathbf{u}(\mathbf{\Omega}_{L},t))\) with the interpolated alternatives. Bi-linear and bi-cubic interpolation are denoted by \(BL(\mathbf{u}(\mathbf{\Omega}_{L},t)),BC(\mathbf{u}(\mathbf{\Omega}_{L},t))\) respectively. We observe that network realisations

Figure 2: Super-resolution results. Physics-constrained convolutional neural network (PC-CNN) compared with interpolation methods. Panel _(i)_ shows the low-resolution input, \(\mathbf{u}(\mathbf{\Omega}_{L},t)\); _(ii)_ bi-linear interpolation, \(BL(\mathbf{u}(\mathbf{\Omega}_{L},t))\); _(iii)_ bi-cubic interpolation, \(BC(\mathbf{u}(\mathbf{\Omega}_{L},t))\); _(iv)_ true high-resolution field, \(\mathbf{u}(\mathbf{\Omega}_{H},t)\); _(v)_ model prediction of the high-resolution field, \(\mathbf{f}_{\mathbf{\theta}}(\mathbf{u}(\mathbf{\Omega}_{L},t))\); and _(vi)_ energy spectra for each of the predictions.
of the high-resolution solution yield qualitatively more accurate results as compared with the interpolation. Artefacts indicative of the interpolation scheme used are present in both of the interpolated fields, whereas the network realisation captures the structures present in the high-resolution field correctly. Across the entire solution domain the model \(\mathbf{f_{\theta}}\) achieves a relative \(\ell^{2}\)-error of \(e=3.449\times 10^{-2}\) compared with \(e=2.091\times 10^{-1}\) for bi-linear interpolation and \(e=1.717\times 10^{-1}\) for bi-cubic interpolation.

Although the relative \(\ell^{2}\)-error provides a notion of predictive accuracy, it is crucial to assess the physical characteristics of the super-resolved field [24]. The energy spectrum, which is characteristic of turbulent flows, represents a multi-scale phenomenon where energy content decreases with the wavenumber. From the energy spectrum of the network's prediction, \(\mathbf{f_{\theta}}(\mathbf{u}(\mathbf{\Omega}_{L},t))\), we gain physical insight into the multi-scale nature of the solution. Results in Figure 2 show that the energy content of the low-resolution field diverges from that of the high-resolution field, which is a consequence of spectral aliasing. Network realisations \(\mathbf{f_{\theta}}(\mathbf{u}(\mathbf{\Omega}_{L},t))\) are capable of capturing finer scales of turbulence compared to both interpolation approaches, prior to diverging from the true spectrum as \(|\mathbf{k}|=18\). The residual loss, \(\mathcal{L_{R}}\), enables the network to act beyond simple interpolation. The network is capable of de-aliasing, thereby inferring unresolved physics. Parametric studies show similar results across a range of super-resolution factors \(\kappa\) (result not shown).

## 7 Conclusions

In this paper, we introduce a method for physics-constrained super-resolution of observations in partial differential equations without access to the full high-resolution samples. First, we define the super-resolution task and introduce the physics-constrained convolutional neural network (PC-CNN), which provides the framework to compute physical residuals for spatiotemporal systems. Second, we formulate an optimisation problem by leveraging knowledge of the partial differential equations and low-resolution observations to regularise the predictions from the network. Third, we showcase the PC-CNN on a turbulent flow, which is a spatiotemporally chaotic solution of the nonlinear partial differential equations of fluid motion (Navier-Stokes). Finally, we demonstrate that the proposed PC-CNN provides more accurate physical results, both qualitatively and quantitatively, as compared to interpolating upsampling methods. This work opens opportunities for the accurate reconstruction of solutions of partial differential equations from sparse observations, as is prevalent in experimental settings, without the full set of high-resolution images.

## 8 Acknowledgements

D. Kelshaw. and L. Magri. acknowledge support from the UK Engineering and Physical Sciences Research Council. L. Magri gratefully acknowledges financial support from the ERC Starting Grant PhyCo 949388.
