Title: DEEP LEARNING FOR PHYSICAL PROCESSES: INCORPORATING PRIOR SCIENTIFIC KNOWLEDGE

Abstract: We consider the use of Deep Learning methods for modeling complex phenomena like those occurring in natural physical processes. With the large amount of data gathered on these phenomena the data intensive paradigm could begin to challenge more traditional approaches elaborated over the years in fields like maths or physics. However, despite considerable successes in a variety of application domains, the machine learning field is not yet ready to handle the level of complexity required by such problems. Using an example application, namely Sea Surface Temperature Prediction, we show how general background knowledge gained from the physics could be used as a guideline for designing efficient Deep Learning models. In order to motivate the approach and to assess its generality we demonstrate a formal link between the solution of a class of differential equations underlying a large family of physical phenomena and the proposed model. Experiments and comparison with series of baselines including a state of the art numerical approach is then provided.

1. INTRODUCTION:

A physical process is a sustained phenomenon marked by gradual changes through a series of states occurring in the physical world. Physicists and environmental scientists attempt to model these processes in a principled way through analytic descriptions of the scientist's prior knowledge of the underlying processes. Conservation laws, physical principles or phenomenological behaviors are generally formalized using differential equations. This physical paradigm has been, and still is the main framework for modeling complex natural phenomena like e.g. those involved in climate. With the availability of very large datasets captured via different types of sensors, this physical modeling paradigm is being challenged by the statistical Machine Learning (ML) paradigm, which offers a prior-agnostic approach. However, despite impressive successes in a variety of domains as demonstrated by the deployment of Deep Learning methods in fields such as vision, language, speech, etc, the statistical approach is not yet ready to challenge the physical paradigm for modeling complex natural phenomena, or at least it has not demonstrated how to. This is a new challenge for this field and an emerging research direction in the ML community. We believe that knowledge and techniques accumulated for modeling physical processes in well developed fields such as maths or physics could be useful as a guideline to design efficient learning systems and conversely, that the ML paradigm could open new directions for modeling such complex phenomena. In this paper we then raise two issues: 1) are modern ML techniques ready to be used to model complex physical phenomena, and 2) how general knowledge gained from the physical modeling paradigm could help designing efficient ML models.

In this work, we tackle these questions by considering a specific physical modeling problem: forecasting sea surface temperature (SST). SST plays a significant role in analyzing and assessing the dynamics of weather and other biological systems. Accurately modeling and predicting such dynamics is critical in various applications such as weather forecasting, or planning of coastal activities. Since 1982, weather satellites have made huge quantities of very high resolution SST data available Bernstein (1982) . Standard physical methods for forecasting SST use coupled ocean-atmosphere prediction systems, based on the Navier Stokes equations. These models rely on multiple physical hypotheses and do not optimally exploit the information available in the data. On the other hand, despite the availability of large amounts of data, direct applications of ML methods do not lead to competitive state of the art results, as will be seen in section 4. We use SST as a typical and representative problem of intermediate complexity. Our goal is not to offer one more solution to this problem, but to use it as an illustration for advancing on the challenges mentioned above. The way we handle this problem is general enough to be transfered to a more general class of transport problems.

We propose a Deep Neural Network (NN) model, inspired from general physical motivations which offers a new approach for solving this family of problems. We first motivate our approach by introducing in section 2 the solution of a general class of partial differential equations (PDE) which is a core component of a large family of transport and propagation phenomena in physics. This general solution is used as a guideline for introducing a Deep Learning architecture for SST prediction which is described in section 3. Experiments and comparison with a series of baselines is introduced in section 4. A review of related work is finally presented in section 5.

The main contributions of this work are: 1) an example showing how to incorporate general physical background for designing a NN aimed at modeling a relatively complex prediction task. We believe the approach to be general enough to be used for a family of transport problems obeying general advection-diffusion principles. 2) formal links between our model's prediction and the solution of a general advection diffusion PDE 3) an unsupervised model for estimating motion fields, given a sequence of images. 4) a proof, on a relatively complex physical modeling problem, that full data intensive approaches based on deep architectures can be competitive with state of the art dedicated numerical methods.

2. PHYSICAL MOTIVATION:

Forecasting consists in predicting future temperature maps using past records. Temperatures are acquired via satellite imagery. If we focus on a specific area, we can formulate the problem as prediction of future temperature images of this area using past images. The classical approach to forecasting SST consists in using numerical models representing prior knowledge on the conservation laws and physical principles, which take the form of PDEs. These models are then coupled with SST data using assimilation techniques in order to adjust to initial conditions. It is then integrated forward in time to predict SST evolution. For the sea surface, temperature variation is related to a fluid transport problem. In fluids, transport occurs through the combination of two principles: advection and diffusion. During advection, a fluid transports some conserved quantity I (the temperature for SST) or material via bulk motion, i.e.for small variations ∆x, ∆t, conservation is expressed as:

EQUATION

applying a first order approximation of the right hand side and moving the resulting terms to the left hand side of equation 1, we obtain the advection equation, known also as the Brightness Constancy Constraint Equation (BCCE):

∂I ∂t + (w.∇)I = 0 (2)

where ∇ denotes the gradient operator, and w the motion vector ∆x ∆t . This equation describes the temporal evolution of quantity I for displacement w. Note that this equation is also the basis for many variational methods for Optical Flow. To retrieve the motion, numerical schemes are applied, and the resulting system of equations, along with a an additional constraint on w is solved for w. This motion can then be used to forecast the future value of I. ∂I ∂t

+ (w.∇)I = D∇ 2 I (3)

∇ 2 denotes the Laplacian operator and D the diffusion coefficient. Note that when D → 0, we recover the advection equation 2. This equation describes a large family of physical processes (e.g. fluid dynamics, heat conduction, wind dynamics, etc). Let us now state a result, characterizing the general solutions of equation 3.

Theorem 1.1 For any initial condition I 0 ∈ L 1 (R 2 ) with I 0 (±∞) = 0, there exists a unique global solution I(x, t) to the advection-diffusion equation 3:

I(x, t) = R 2 k( x -w, y) I 0 (y) dy (4)

where k(u, v) = 1 4πDt e -1 4Dt u-v 2 is a radial basis function kernel, or alternatively, a 2 dimensional Gaussian probability density with mean u and variance 2Dt.

Equation 4 provides a principled way to calculate I(x, t) for any time t using the initial condition I 0 , provided the motion w and the diffusion coefficient D are known. It states that quantity I(x, t) can be computed from the initial condition I 0 via a convolution with a Gaussian probability density function. In other words, if I was used as a model for the evolution of the SST and the surface's underlying advecting mechanisms were known, future surface temperatures could be predicted from previous ones. Unfortunately neither the initial conditions, the motion vector nor the diffusion coefficient are known. They have to be estimated from the data. Inspired from the general form of solution 4, we propose a ML method, expressed as a Deep Learning architecture for predicting SST. This model will learn to predict a motion field analog to the w in equation 4, which will be used to predict future images. ) with a convolutional neural network (top left CDNN component). A warping scheme then displaces the last input image along this motion estimate to produce the future image. The error signal is calculated using the target future image I t+1 , and is backprogated through the warping scheme to correct the CDNN. To produce multiple time-step forecasts, the predicted image is fed back in the CDNN in an autoregressive manner.

3. MODEL:

The model consists of two main components, as illustrated in Figure 1 . One predicts the motion field from a sequence of past input images, this is convolutional-deconvolutional (CDNN) module on the top of figure 1 , and the other warps the last input image using the motion field from the first component, in order to produce an image forecast. The entire system is trained in an end-to-end fashion, using only the supervision from the target SST image. By doing so, we are able to produce an interpretable latent state which corresponds in our problem to the velocity field advecting the temperatures.

Let us first introduce some notations. Each SST image I t is acquired on a bounded rectangle of R 2 , named Ω. We denote I t (x) and w t (x) the sea surface temperature and the two-dimensional motion vector at time t ∈ R at position x ∈ Ω. I t : Ω → R and w t : Ω → R 2 represent the temperatures and the motion vector field at time t defined on Ω. When time t and position x are available from the context, we will drop the subscript t from w t (x) and I t (x), along with x for clarity. Given a sequence of k consecutive SST images {I t-k-1 , ..., I t } (also denoted as I t-k-1:t ), our goal is to predict the next image I t+1 . Generally, and this is the case for our problem, we do not have a direct supervision on the motion vector field, since the target motion is usually not available. Using the warping scheme introduced below, we will nonetheless be able to (weakly) supervise w, based on the discrepancy of the warped version of the I t image and the target image I t+1 . To calculate the pixel value for time t+1 at position x, we first compute its previous position at time t, i.e. xw. We then center a Gaussian in that position in order to obtain a weight value for each pixel in I t based on its distance with xw, and compute a weighted average of the pixel values of I t . This weighted average will correspond to the new pixel value at x in I t+1 .

3.2. WARPING SCHEME:

Discretizing the solution of the advection-diffusion equation in section 2 by replacing the integral with a sum, and setting image I t as the initial condition, we obtain a method to calculate the future image, based on the motion field estimate ŵ. The latter is used as a warping scheme:

EQUATION

where

k(x -ŵ, y) = 1 4πD∆t e -1 4D∆t x-ŵ-y 2

is a radial basis function kernel, as in equation 4, parameterized by the diffusion coefficient D and the time step value ∆t between t and t + 1 and ŵ is the estimated value of the vector flow w. To calculate the temperature for time t + 1 at position x, we compute the scalar product between k(xŵ, .), a Gaussian centered in xŵ, and the previous image I t . Simply put, it is a weighted average of the temperatures I t , where the weight values are larger when the pixel's positions that are closer to xŵ. Informally, xŵ corresponds to the pixel's previous position at time t. See figure 3 .

As seen by the relation with the solution of the advection-diffusion equation, the proposed warping mechanism is then clearly adapted to the modeling of phenomena governed by the advectiondiffusion equation. SST forecasting is a particular case, but the proposed scheme can be used for any problems in which advection and diffusion are occurring. Moreover, this warping scheme is entirely differentiable, allowing backpropagation of the error signal to the motion fireld estimating module.

This warping mechanism has been inspired by the Spatial Transformer Network (STN) Jaderberg et al. (2015) , originally designed to be incorporated as a layer in a convolutional neural network architecture in order to gain invariance under geometric transformations. Using the notations in Jaderberg et al. (2015) , when the inverse geometric transformation T θ of the grid generator step is set to T θ (x) = xŵ(x), and the kernels k( . ; Φ x ) and k( . ; Φ y ) in the sampling step are radial basis function kernels, we recover our warping scheme. The latter can be seen as a specific case of the STN, without the localization step. This result theoretically grounds the use of the STN for Optical Flow in many recent articles Zhu et al. (2017) For training, supervision is provided at the output of the warping module. It consists in minimizing the discrepancy between the warped image Ît+1 and the target image I t+1 . The loss is measured via a differentiable function and the gradient is back propagated through the warping function in order to adjust the parameters of the convolutional-deconvolutional module generating the vector field. This is detailed in the next section.

3.3. LOSS FUNCTION:

At each iteration, the model aims at forecasting the next observation, given the previous ones. We evaluate the discrepancy between the warped image Ît+1 and the target image I t+1 using the Charbonnier penalty function ρ(x) = (x + ) 1 α , where and α are parameters to be set. Note that with = 0 and α = 1 2 , we recover the 2 loss. The Charbonnier penalty function is known to reduce the influence of outliers compared to an l 2 norm. We have also tested the Laplacian pyramid loss Ling & Okada (2006) , where we enforce convolutions of all deconvolutional layers to be close to down-sampled versions of the target image in the Charbonnier penalty sense, but we have observed an overall decrease in generalization performance.

The proposed NN model has been designed according to the intuition gained from general background knowledge of a physical phenomenon, here advection-diffusion equations. Additional prior knowledge -expressed as partial differential equations, or through constraints -can be easily incorporated in our model, by adding penalty terms in the loss function. As the displacement w is explicitly part of our model, one strength of our model is its capacity to apply some regularization term directly on the motion field. In our experiments, we tested the influence of different terms: divergence ∇. w t (x)2 , magnitude w t (x)

2 and smoothness ∇w t (x) 2 .

L t = x∈Ω ρ( Ît+1 (x) -I t+1 (x)) + λ div (∇. w t (x)) 2 + λ magn w t (x) 2 + λ grad ∇w t (x) 2 (6) 4 EXPERIMENTS 4.1 DATASET DESCRIPTION

Since 1982, high resolution SST data has been made available by the NOAA6 weather satellite, Bernstein (1982) . Dealing directly with these data requires a lot of preprocessing (e.g. some regions are not available due to clouds, hindering temperature acquisition). In order to avoid such complications which are beyond the scope of this work, we used synthetic but realistic SST data of the Atlantic ocean generated by a sophisticated simulation engine: NEMO (Nucleus for European Modeling of the Ocean) engine 2 , Madec (2008) . NEMO is a state-of-the-art modelling framework of ocean related engines. It is a primitive equation model adapted to the regional and global ocean circulation problems. Historical data is accumulated in the model to generate a synthesized estimate of the states of the system using data analysis, a specific data assimilation scheme, which means that the data does follow the true temperatures. The resulting dataset is constituted of daily temperature acquisitions of 481 by 781 pixels, from 2006-12-28 to 2017-04-05 (3734 acquisitions).

We We make the simplifying hypothesis that the data in a single sub-region contains enough information to forecast the future of the sub-region. As the forecast is for a small temporal horizon we can assume that the influence from outside the region is small enough. We normalize the daily SST acquisitions of each sub region using the mean and the standard deviation of all the SST data of the sub-region acquired on the same day of the year for all the years in the training set, i.e. the SST acquisition of sub-region 2 on date September 8th 2017 is standardized using the data of all the September 8th available in the dataset. This removes the seasonal component from SST data.

4.2. BASELINE COMPARISON:

We compare our model with several baselines. Each model is evaluated with a mean square error metric, forecasting images on a horizon of 6 (we forecast from I t+1 to I t+6 and then average the MSE). The hyperparameters are tuned using the validation set. Neural network based models are run on a Titan Xp GPU, and runtime is given for comparison purpose.

Concerning the constraints on the vector field w (equation 6. the regularization coefficients selected via validation are λ div = 1, λ magn = -0.03 and λ grad = 0.4. The coefficient diffusion D was set to 0.45 by cross validation. We also compare the results with the model without any regularization.

Our reference model for forecasting is Béréziat & Herlin (2015) , a numerical assimilation model which relies on data assimilation. In Béréziat & Herlin (2015) , the ocean's dynamics are modeled using shallow water equations Vallis (2017) and the initial conditions, along with other terms, are estimated using assimilation techniques Trémolet (2006) . This is a state of the art assimilation model for predicting ocean dynamics, here SST.

The other baselines are 1) an autoregressive convolutional-deconvolutional NN (ACNN), with an architecture similar to our CDNN module, but trained to predict the future image directly, without explicitly representing the motion vector field. Each past observation is used as an input channel (the 4 input images used in the experiments are concatenated), and the output is used as new input for multi step forecasting, 2) a ConvLSTM model Shi et al. (2015) , which uses convolutional transitions in the inner LSTM module, and 3) the model in Mathieu et al. (2015) which is a multi-scale ACNN trained as a Generative Adversial Network (GAN). We have used a non-official code for Mathieu et al. (2015) , which is made available at https://github.com/dyelax/Adversarial_ Video_Generation. For Béréziat & Herlin (2015) , the code has been provided by the authors of the paper. We have implemented the ACNN and ConvLSTM models ourselves. The code for our models, along with these baselines will be made available. The regularization coefficients for our model have been set using a validation set with λ div = 1, λ magn = -0.03 and λ grad = 0.4.

4.3. QUANTITATIVE RESULTS:



Model:

Quantitatively, our model performs well. The MSE score is better than any of the baselines.

The closest NN baseline is Mathieu et al. (2015) which regularizes a regression convolutiondeconvolution model with a GAN. The performance is however clearly below the proposed model and it does not allow to easily incorporate prior constraints inspired from the physics of the phenomenon. ACNN is a direct predictor of the image sequence, implemented via a CDNN module identical to the one used in our model. Its performance is poor. Clearly, a straightforward use of prediction models is not adapted to the complexity of the phenomenon. ConvLSTM performs better: as opposed to the ACNN, it seems to be able to capture a dynamic, although not very accurately. Overall, direct prediction models are not able to capture the complex underlying dynamics and they produce blurry sequences of images. The GAN explicitly forces the network output to eliminate the blurring effect and then makes it able to capture short term dynamics. The state of the art numerical model Béréziat & Herlin (2015) , performs well but has a slighthly lower performance than our regularized model, although it incorporates more prior constraints. This shows that pure ML models, when conceived adequately and when trained with enough data, can be competitive with state of the art dedicated models. Regularizing the motion vector w notably increases the performance w.r.t. to the unregularized model. The choice of the constraints (divergence, magnitude and smoothness) inspired here by physical background correspond to relevant priors on the dynamics of the model.

I t I t+1 I t+3 I t+6

As for the running time, the proposed model is extremely fast, being just above the ConvLSTM model of Shi et al. (2015) . The running time of Béréziat & Herlin (2015) 's model is not comparable to the others. It was run on a CPU (no GPU code) when all the others were run on Titan Xp GPU. However, an optimization procedure is required to estimate the motion field, and it is clearly slower than the straightforward NN predictions. Moreover, in order to prevent the numerical scheme from diverging, multiple intermediate forecasts are required.

Besides MSE, we need to analyze the prediction samples qualitatively. Figure 4 .3 shows the predictions obtained by the different models. On the top row, the ground truth for a sequence of 4 temperature images corresponding to time t, t + 1, t + 3 and t + 6. The second row corresponds to our regularized model prediction at times t + 1, t + 3 and t + 6 (time t corresponds to the last input image, it is repeated on each row). The model seems to conserve temperatures. The prediction is close to the target for t + 1, t + 3 and starts to move away at time t + 6. The third row shows the motion flow estimated by the model. Each color in the flow images corresponds to a motion vector. There is clearly a strong evolving dynamic captured for this sequence. Row 4 is the numerical assimilation model of Béréziat & Herlin (2015) . It also clearly captures some dynamics and shows interesting patterns, but it tends to diverge when the prediction horizon increases. The ACNN model (row 5) rapidly produces blurry images; it does not preserve the temperatures and does not seem to capture any dynamics. On row 6 are plotted the predictions of the ConvLSTM model. Temperature is not preserved and although a dynamic is captured, it does not correspond to the target. Overall, the proposed model seems to forecast SST quite accurately, while retrieving a coherent motion vector field.

5. RELATED WORK:

ML for Physical modeling Close to this work is the field of spatio-temporal statistics. In their reference book Cressie & Wikle (2015) also advocate the use of physical background knowledge to build statistical models. They show how the design of statistical models can be inspired from partial differential equations linked to an observed physical phenomenon. They mainly consider auto-regressive models within a hierarchical Bayesian framework. In Raissi et al. (2017) , Archambeau et al. (2007) and Alvarez et al. (2011) the author use PDE-inspired gaussian process to model physical process.

Even if the methods and the application are different, the motivation and arguments are similar to the ones developed here.

Another interesting research direction is the use of NNs for reducing the complexity of numerical simulation for physical processes. Generally, in these approaches statistical models are used in place of a computational demanding component of the numerical simulation process. For example in the domain of fluid dynamics, Tompson et al. (2017) and Ladický et al. (2015) propose to use regressors for simulating fluid and smoke animation. Ladický et al. (2015) use a random forest to compute particle location and Tompson et al. (2017) use a CNN to approximate part of a numerical PDE scheme. In these approaches, ML is only a component of a numerical simulation scheme whereas we aim at modeling the whole physical process via a Deep Learning approach. Farther to our objective, Rudy et al. ( 2017) make use of a sparse regression method for discovering the governing partial differential equation(s) of a given system by time series measurements in the spatial domain.

Other works have suggested using neural networks for physical process forecast, such as Brajard et al. (2017) .

Our work is also related to recent developments in computer vision, in the related but distinct fields of video prediction and motion estimation in videos. Our goal and the domain of application are clearly different from video modeling, but since our solution involves predicting a motion field and the next SST image, the solutions share some similarities. Motion estimation and video predictions by deep architectures have motivated a series of work over the last two years. We briefly review them below and outline the differences.

Optical Flow Optical flow consists in retrieving the apparent motion of objects, surfaces, or particles between two consecutive frames of a video. The extracted motion can be used in many areas such as object detection, object tracking, movement detection, robot navigation and visual odometry. In the vision community, this is considered as a problem by itself and several papers are dedicated to this topic. Classical methods rely on the brightness constancy constrain equation (BCCE) (equation 2), derived from the observation that surfaces usually persist over time and hence the intensity value of a small region remains the same despite its position change Sun et al. (2008) . Since using BCCE directly leads to complicated optimizing issues, classic approaches -namely differential methodsapproximate BCCE using a first order Taylor expansion and develop variational methods.

As an alternative to these methods, Deep Learning models have been recently proposed for estimating the optical flow between 2 images. Fischer et al. (2015) formulate optical flow as a supervised regression problem, using a CNN to predict motion. Ilg et al. (2016) build on this approach and propose to use an ensemble of these CNN architectures. They assess results on par with state of the art methods for optical flow, while maintaining a small computational overhead. The difficulty here is that these methods require a notable quantity of target data, i.e. optical flow images, while because of the complexity of manually annotating flow images, there are only a few small annotated collections available. Fischer et al. (2015) and Ilg et al. (2016) chose to pretrain their model on a synthetic dataset made of computer animations and their associated motion and show that this information transfers well to real videos. Yu et al. (2016) demonstrate that it is possible to predict the optical flow between two input images in an unsupervised way using a CNN and a Spatial Transformer Network. This is however not extensible for prediction as is done in our setting since this requires the two images I t and I t+1 as input while I t+1 is not available at inference time for prediction.

Video prediction:

It is only very recently that video prediction emerged as a task in the Deep Learning community. For this task, people are generally interested at predicting accurately the displacement/ emergence/ disappearing of objects in the video. In our application, the goal is clearly different since we are interested into modeling the whole dynamics behind image changes and not at following moving objects. Let us first introduce some methods that perform prediction by computing optical flow or a similar transformation. Both Patraucean et al. (2015) and Finn et al. (2016) use some form of motion flow estimation. For next frame prediction Patraucean et al. (2015) introduce a STN module at the hidden layer of a LSTM in order do estimate a motion field in this latent space. The resulting image is then decoded in the original image space for prediction. This approach clearly does not allow introducing prior knowledge on the field vector as this has been done in our work. Finn et al. (2016) learn affine transformations on image parts in order to predict object displacement and Van Amersfoort et al. ( 2017) proposed a similar model.

Let us now consider models that directly attempt to predict the next frame without estimating a motion field. As shown in the experimental section, plain application of autoregressive models produces blurred images. Mathieu et al. (2015) , one of our baseline proposed to use different loss functions and a GAN regularization of a CDNN predictor which led to sharper and higher quality predictions. Significant improvements have been obtained with the Video Pixel Network of Kalchbrenner et al. (2016) , which is a sophisticated architecture composed of resolution preserving CNN encoders, LSTM and PixelCNN decoders which form a conditional Spatio-temporal video autoencoder with differentiable memory. This model is probably state of the art today for video prediction, They reach a high accuracy on moving MNIST and good performance on a robot video dataset. A drawback is the complexity of the model and the number of parameters: they are using respectively 20 M and 1 M frames on these two datasets. We did not test this model since up to our knowledge no code was available.

6. CONCLUSION AND FUTURE WORK:

The evolution in time of the proposed model is deterministic. Predicting future observations should also deal with the inherent ambiguity and lack of information for the prediction task. A natural future direction would be to incorporate uncertainty in the model's evolution in the proposed framework. We can extend the proposed model by incorporating a stochastic latent variable in the flow field generation process. A promising direction is the development of generative models which has become popular in Deep Learning, leading to different families of innovative models. For example, the Stochastic Gradient Variational Bayes algorithm (SGVB) Kingma & Welling (2014) provides a framework for learning stochastic latent variables with deep neural networks, and has recently been used by some authors to model time series Karl et al. (2016) ; Chung et al. (2015) ; Krishnan et al. (2015) . A recent work where both spatial and temporal information are considered is Walker et al. (2016) who model pixel trajectories in a video. As a follow up of our work, we plan to consider such extensions in the future.

The data intensive paradigm offers alternative directions to the classical physical approaches for modeling complex natural processes. Our belief is that cross fertilization of both paradigms is essential for pushing further the frontier of complex data modeling. By using as an example application a relatively complex problem concerning ocean dynamics, we proposed a principled way to design Deep Learning models using inspiration from the physics. The proposed approach can be easily generalized to a class of problems for which the underlying dynamics follow advection-diffusion principles. We have compared the proposed approach to a series of baselines. It is able to reach performance comparable to a state of the art numerical model and clearly outperform alternative NN models used as baselines.

A PROOF OF THE THEOREM IN SECTION 2 1

Proof. In the following, bold x and y will denote vectors of R 2 , while x and y will correspond to the first and second components of x, respectively. Analogously, u and v will correspond to the components of w. The 2D Fourier Transformation F of f : R 2 → R is defined as

EQUATION

We apply the Fourier Transform F to both sides of 3. As consequence of the linearity of the Fourier transform, we can calculate decompose the Fourier transform of the left hand side in the sum of the transforms of each term. We have three terms: ∂I ∂t , (w.∇)I and -D∇ 2 I.

EQUATION

Regrouping all three previously calculated terms, we obtain

∂F(I) ∂t + (i < ξ, w > +D ξ 2 )F(I) = 0 (11)

This is a first order ordinary differential equation of the form f (t) + af (t) = 0, which admits a known solution f (t) = f (0)e -at . Thus, the solution of 11 is

F(I) = F(I) 0 e -(i<ξ,w>+D ξ 2 )t = F(I) 0 e -i<ξ,w>t e -Dt ξ 2 (12)

where F(I) 0 denotes the initial condition of the advection diffusion equation in the frequency domain. In order to obtain a solution of 3 in the spatial domain, we calculate the inverse Fourier Transform F -1 of 12. The multiplication of two functions in the frequency domain is equivalent to their convolution in the spatial domain, i.e. F(f * g) = F(f )F(g). Hence, the inverse of both terms F(I) 0 e -i<ξ,w>t and e -Dt ξ 2 can be calculated separately:

Multiplication by a complex exponential in the frequency domain is equivalent to a shift in the spatial domain : e -i<ξ,w> F(f (x)) = F(f (xw)), for v ∈ R 2 . Thus, for the first term,

EQUATION

For the second term, we use the fact that the Fourier Transform of a Gaussian function also is a Gaussian function, i.e. F( 1 2πσ 2 e -1 2σ 2 x 2

) = e -1 2 σ 2 ξ 2

. Identifying σ 2 with 2Dt, we have:

EQUATION

As has been stated above, the solution is a convolution of both previously calculated terms:

I(x, t) = R 2 1 4πDt

e -1 4Dt y 2 I 0 (xwy)dy

= R 2 1 4πDt

e -1 4Dt x-w-y 2 I 0 (y)dy

B ON THE GENERALIZATION IN SPACE AND TIME:

The ability of the model to adapt to other conditions should be evaluated on other regions. This, however, requires a complete study by itself and is beyond the scope of this paper. We, however, present below complementary experiments aimed at assessing the potential of the proposed model for forecasting SST on sequences distant in time and space from the ones used for training. We describe below some results providing indications on how the model performs on regions different from the training ones. For these experiments, the model has been trained on regions 17 and 18 in Figure 4 .1 and tested on two other regions (regions 8 and 9), and vice versa (trained on 8 and 9 and tested on 17 and 18). The two couples of regions have been selected so as to have different latitude and longitude. The underlying physical processes generating the data are known to be different in these regions: the overall motion in regions 17 and 18 is greater, and the difference between extreme 

:

A proof of the theorem is provided in appendix A

NEMO data are available at http://marine.copernicus.eu/services-portfolio/ access-to-products/?option=com_csw&view=details&product_id=GLOBAL_ ANALYSIS_FORECAST_PHY_001_024

non numbered regions correspond to land and not sea on the figure

Figures and Tables:

Figure1: Motion is estimated from the input images (I t-k-1:t ) with a convolutional neural network (top left CDNN component). A warping scheme then displaces the last input image along this motion estimate to produce the future image. The error signal is calculated using the target future image I t+1 , and is backprogated through the warping scheme to correct the CDNN. To produce multiple time-step forecasts, the predicted image is fed back in the CDNN in an autoregressive manner.

Figure 2: Architecture of the CDNN motion estimation component. For the estimated motion flow ŵt , colours correspond to the flow orientation and colour intensity to the flow intensity

Figure3: Warping scheme. To calculate the pixel value for time t+1 at position x, we first compute its previous position at time t, i.e. xw. We then center a Gaussian in that position in order to obtain a weight value for each pixel in I t based on its distance with xw, and compute a weighted average of the pixel values of I t . This weighted average will correspond to the new pixel value at x in I t+1 .

, Yu et al. (2016), Patraucean et al. (2015), Finn et al. (2016): in equation 3, when D → 0, we recover the brightness constancy constraint equation, used in the latter.

extract 64 by 64 pixel sized sub-regions as indicated in figure 4.1. We use data from years 2006 to 2015 for training and validation (94743 training examples), and years 2016 to 2017 for testing. We withhold 20% of the training data for validation, selected uniformly at random at the beginning of each experiment. For the tests we used sub-regions enumerated 17 to 20 in figure 4.1, where the interactions between hot and cold waters make the dynamics interesting to study. All the regions numbered in figure 4.1, from 2006 to 2015 where used for training 3 . Each sequence of images used for training or for evaluation corresponds to a specific numbered sub-region.

Figure 4: Sub regions extracted for the dataset. Test regions are regions 17 to 20.

Figure 5: From top to bottom: target, our model prediction, our model flow, numerical assimilation model , ACNN, ConvLSTM. Data correspond to daily temperatures from January 17 to January 23, 2017

TEMPORAL DIMENSION In section 4, training has been performed on data from 2006 -2015 and testing on the period 2016-2017. In order to provide some indication of the model behavior on more distant time intervals between train and test data, we have performed experiments using the same regions (17 to 20) as in section 4, but using the period 2011 to 2017 for training and period 2006 to 2010 for testing. Figure B.1 shows the MSE curve on this test set, each point corresponding to the mean MSE on predictions performed on 6 days ahead the current date. The most important conclusion is probably that the MSE error remains in the same range for all these years. All the yearly error curve show a clear seasonal phenomenon with a higher prediction error during summer. A similar behavior has been observed when exchanging train and test data.

Figure 6: Evaluation of our model's accuracy in time on data from 2006 to 2010 using data from 2011 to 2017 for training. Regions 17 to 20 were used for both periods. Each day, we produce daily forecasts for 6 days ahead and calculate the associated mean square error.

Figure 7: Output for the 6 of May to the 9 of May 2016, Output , From top to bottom: target, our model prediction, our model flow

Average score and average time on test data. Average score is calculated using the mean square error metric (MSE), time is in seconds.

Table content: <table><tr><td>Average Score (MSE) Average Time</td></tr></table>

Evaluation of our model's spatial generalization ability. We train our model on two distinct regions and calculate the MSE on both regions for each trained model.

Table content: <table><tr><td/><td colspan="2">Test Regions 17 &amp; 18 Test Regions 8 &amp; 9</td></tr><tr><td>Model trained on Regions 17 &amp; 18</td><td>1.43</td><td>1.22</td></tr><tr><td>Model trained on Regions 8 &amp; 9</td><td>1.90</td><td>1.19</td></tr></table>

Peer review (iclr2018):
 The authors use deep learning to learn a surrogate model for the motion vector in the advection-diffusion equation that they use to forecast sea surface temperature. In particular, they use a CNN encoder-decoder to learn a motion field, and a warping function from the last component to provide forecasting. I like the idea of using deep learning for physical equations. I would like to see a description of the algorithm with the pseudo-code in order to understand the flow of the method. I got confused at several points because it was not clear what was exactly being estimated with the CNN. Having an algorithmic environment would make the description easier. I know that authors are going to publish the code, but this is not enough at this point of the revision. Physical processes in Machine learning have been studied from the perspective of Gaussian processes. Just to mention a couple of references “Linear latent force models using Gaussian processes” and "Numerical Gaussian Processes for Time-dependent and Non-linear Partial Differential Equations" In Theorem 2, do you need to care about boundary conditions for your equation? I didn’t see any mention to those in the definition for I(x,t). You only mention initial conditions. How do you estimate the diffusion parameter D? Are you assuming isotropic diffusion? Is that realistic? Can you provide more details about how you run the data assimilation model in the experiments? Did you use your own code?