{'fraction_non_alphanumeric': 0.036307789241553125, 'fraction_numerical': 0.013265225756365396, 'mean_word_length': 4.646401118099232, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Physically interpretable machine learning algorithm on multidimensional non-linear fields

In an ever-increasing interest for Machine Learning (ML) and a favorable data development context, we here propose an original methodology for data-based prediction of two-dimensional physical fields. Polynomial Chaos Expansion (PCE), widely used in the Uncertainty Quantification community (UQ), has recently shown promising prediction characteristics for one-dimensional problems, with advantages that are inherent to the method such as its explicitness and adaptability to small training sets, in addition to the associated probabilistic framework. Simultaneously, Dimensionality Reduction (DR) techniques are increasingly used for pattern recognition and data compression and have gained interest due to improved data quality. In this study, the interest of Proper Orthogonal Decomposition (POD) for the construction of a statistical predictive model is demonstrated. Both POD and PCE have widely proved their worth in their respective frameworks. The goal of the present paper was to combine them for a field-measurement-based forecasting. The described steps are also useful to analyze the data. Some challenging issues encountered when using multidimensional field measurements are addressed, for example when dealing with few data. The POD-PCE coupling methodology is presented, with particular focus on input data characteristics and training-set choice. A simple methodology for evaluating the importance of each physical parameter is proposed for the PCE model and extended to the POD-PCE coupling.


Introduction
Deep Learning techniques (DL [29,59]) and more generally Machine Learning (ML [79,95]), and their applications to physical problems (fluid mechanics [9,51,67,74] ; aerodynamics [110,115] ; plasma physics [28,83] ; astrophysics and astronomy [50,106] ; particle physics [2] ; quantum mechanics [70], geosciences [27,46,86,88,92]) have made a promising take-off in the last few years. This has been particularly the case for fields where the measurement potential has dramatically increased, with increasing spatiotemporal resolution (e.g. Geoscience Data [46] ; the new SWOT satellite mission [71,78]). In particular, multi-layer Neural Networks (NN) [91] are widely used for these applications. These techniques are of interest when data are available for a real-world problem that is difficult to model using process-based equations, because prescribing the underlying relationships is complex [95]. For such cases, the learning procedure is performed through a combination of steps: Encoding, Hidden Layers or Latent Representation, Decoding [59,95]. In particular, transformation functions, also called Activation Functions (AFs), are used jointly with weight matrices, to transform the data from one layer to another. The popularity of NN comes from this complex structure, which makes them adaptable for various applications [1,94]. Indeed, the combination of AF transformations helps capture complex interactions and non-linearities, making NN a widely used approach [95].
However, some limitations prevent the use of NN for physical applications.
(i) It is difficult to provide an explicit input-to-output formulation, due to the combinations and transformations involved. Physical interpretation of the constructed model is therefore tedious. This is why NN have been applied in physics often as a black-box and rather to optimize the accuracy of prediction than to extract physical information [41]. (ii) The calibration of the NN model involves various choices, such as the number of layers, number of neurons, the AFs, the weight matrices, etc. The number of hyper-parameters grows with the number of neurons and layers (curse of dimensionality), which limits the maximum complexity and therefore the interest of NN for highly non linear problems [79]. The choice of AF is application and data-set dependent [55]. (iii) The theoretical ability of NN to approximate functions has been proven with the fundamental Universal Approximation Theorem either with a limited number of layers [5,12,18,23,38] or with a particular type of AFs (ReLu functions in [13,21,34,64,87]). This limits the utility of NN for non-linear physical applications [1,94]. Consequently, simple networks (with few layers) [41] are often used for physical problems, and the necessary assumptions on the AFs are often overlooked in favor of efficiency [55].
To overcome these limitations, we propose an alternative ML, based on a coupling between Proper Orthogonal Decomposition (POD) [14] and Polynomial Chaos Expansion (PCE) [57,58], suitable for the prediction of spatiallydistributed physical fields. Here, POD is used for both Encoding and Decoding whereas PCE is used as a Latent Representation, as represented in Figure 1. The proposed POD-PCE addresses these drawbacks of ML.
(i) It is explicit and simple to implement, as it consists of the association of two linear decompositions. POD is a linear separation of the spatiotemporal patterns [66], shown to be accurate for both linear and non-linear problems [101], combining and relevance. PCE is a well-established method in Uncertainty Quantification (UQ) [65,100,107,112,114], widely used for the study of stochastic behavior in physics [8,45,80,102]. It is a linear polynomial expansion that allows non-linearities to be gradually added to the model by increasing the polynomial degree. The linearity and orthonormality of the POD and PCE components and the probabilistic framework of PCE make the output's statistical moments easier to study [40,49,63,99], enabling straightforward physical interpretation of the model [10,24,52]. (ii) It only has two hyper-parameters: a number of POD components, and a PCE polynomial degree. Both can be chosen according to quantitative criteria [6,14]. All other forms of parameterization (choice of the polynomial basis) can be achieved with robust physical and/or statistical arguments [65,97,107], as assessed in the present paper. Furthermore, the orthonormality of the POD and PCE bases minimizes the number of components necessary to capture essential variations in data. Additionally, the POD modes capture more energy than any other decomposition [16,73], PCE is known to exponentially converge with polynomial degree [58,113], and the cardinality of the latter can be reduced by sparse basis selection [6]. (iii) It can be considered as a universal expansion for physical field approximation: a physical field has a finite variance, which implies that it belongs to the Hilbert space of random variables with finite second order moments. There therefore exists a numerable set of orthogonal random variables, that form the basis of this Hilbert space, on which the field of interest can be expanded (strict equality, not approximation) [100]. A mathematical setting for basis construction based on input was established by Soize and Ghanem [97] for the general case of dependent variables with arbitrary density, provided that the set of inputs is finite.
The present study consisted in: i) combined use of POD and PCE in ML for point-wise prediction ; ii) application to field data with the inherent challenges not encountered with numerical data (e.g. paucity) ; iii) a focus on model explicitness as a key condition for physical understanding and iv) the influence of forcing variables study, based on a classical measure of importance (Garson weights [25]) directly computed with the POD-PCE expansion coefficients. Firstly, associating regression techniques to POD, and more generally to Reduced Order Models (ROM), is not novel [32,53,76,108]. The cited studies, however, focused on dimensionality reduction, whereas the explicit formulation and applicability to complex physical processes are emphasized in the present study. Secondly, coupling PCE to ROM was recently addressed [36,54,77] and the use of PCE as ML is consistent with the methodology proposed by Torre et al. [103], where they showed that PCE is as powerful as classical ML techniques, but neither spatiotemporal fields nor physical interpretability were addressed. The data in these studies were either obtained from numerical experiments, emulated from analytical benchmark functions such as Sobol or Ishigami, or based on one-dimensional data sets [103]. Our methodology is assessed on high-resolution two-dimensional field measured data.
The assessment of the proposed methodology was based on the study of sedimentation processes in a cooling water intake located in a coastal area, subject to tide and wave forcing. This application in geosciences is characterized by high non-linearities and is conditioned a priori with various parameters [4]. The disparity of the space and time scales and the interaction between various processes make process-based modeling difficult [4,82]. A complete overview of the use of ML in coastal sediment transport modeling is given in [27].
The paper is organized as follows. Section 2 gives a detailed explanation of the methodology. Section 3 deals with assessment of the methodology based on a physical application. Firstly, the study case and data are described in 3.1. POD and PCE performances are then demonstrated independently in 3.2 with a deep physical analysis using adequate measures. The performance of the POD-PCE predictor is discussed in 3.3. A summary of the study and perspectives of the proposed methodology are presented in Section 4.

Theoretical framework
In this section, the objective is to define the framework of the proposed POD-PCE Machine Learning methodology, along with physical influence indicators for the inputs. This is the object of Subsection 2.3, but first, a reminder of the existing POD and PCE theoretical bases is presented in 2.1 and 2.2 respectively.

Proper Orthogonal Decomposition
POD is a dimensionality reduction technique [66] that is well documented in literature [14,61,62,101]. Theoretical details and demonstrations can be found in descending chronological order in [3,16,73]. For clarity's sake, the essential elements of POD are summarized below.
The goal of POD is to extract the spatial patterns of a continuous time and space function. These patterns, when added and multiplied by appropriate temporal coefficients, explain the dynamics of the variable of interest: a real-valued physical field.
Let u : Ω × T → D be a continuous function of two variables (x, t) ∈ Ω × T. The following relationships and properties hold for any Ω × T and Hilbert space D characterized by its scalar product (. , .) D and induced norm ||.|| D . However, and as is the case for a majority of physical fields, we shall consider Ω as a set of spatial coordinates (e.g. R 2 or R 3 ), T as a set of time coordinates (e.g. R + or subset [0, T ]), and D as a set of scalar real values or vector real values (e.g. R or R 2 ). POD consists in an approximation of u(x, t) at a given order d ∈ N (Lumley [66]) as in Equation 1, where {v k (.)} d k=1 ⊂ C(T, R) and {φ k (.)} d k=1 ⊂ C(Ω, D), with C(A, B) denoting the space of continuous functions defined over A and arriving at B. The objective of POD is to identify {φ k (.)} d k=1 that minimizes the distance of the approximation from the true value u(., .), over the whole Ω × T domain, with an orthogonality constraint for {φ k (.)} d k=1 using the scalar product (. , .) D . This can be defined, in the least-squares sense, as a minimization problem.
The minimization problem is defined for all orders d ∈ N, so that the members φ k are ordered according to their importance. In particular, for order 1, φ 1 is the linear generator of the sub-vector space most representative of u(x, t) in D. For D = Im(u), the family {φ k (.)} d k=1 is called the POD basis of D of rank d. The solution to this problem has already been established in literature [66,96]. The theoretical aspects of POD and demonstrations of mathematical properties can, for example, be found in [73]: the POD basis of D of order d is the orthonormal set of eigenvectors of an operator R : D → D defined as Rφ = (u, φ) D × u T , if the eigenvectors are taken in decreasing order of the corresponding eigenvalues {λ k } d k=1 .
For this expansion, an accuracy rate, also called the Explained Variance Rate (EVR), denoted e d at rank d, can be calculated as in Equation 2 . EVR tends to 1 (perfect approximation) when d → +∞.
In practice, for D = R, when u(., .) is a discrete sample on a set of m ∈ N space coordinates X = {x 1 , . . . , x m } and for n ∈ N measurement times T = {t 1 , . . . , t n }, the available data set is arranged in a matrix U(X , T ) = [u(x i , t j )] i,j ∈ R m×n , called the snapshot matrix, so as to be able to work in a discrete space. The POD problem formulated in Equation 1 can be written in its discrete form as U(X , The problem can therefore be viewed as if working with a new function U(X , .) = [u(x i , .)] i∈{1,...,m} : T → D = R M . Then, the average over T can be defined as the statistical mean over the subset T , and the scalar product (. , .) D as the canonical product over R m . The POD operator R can be written as in Equation 3, where U(X , t j ) = [u(x i , .)] i∈{1,...,m} is the column number j of the matrix U(X , T ) (i.e the measurement over X at time t j ), and Φ(X ) = [φ(x i )] i∈{1,...,m} . As finding the POD basis is equivalent to identifying the orthonormal set of eigenvectors of the operator R, then for this discrete representation the problem becomes equivalent to solving the eigen problem of the matrix R := 1 n U(X , T )U(X , T ) T , called the covariance matrix. A number d ∈ N of eigen vectors Φ(X ) are identified and stored in the columns of the matrix Φ (d) (X ). For the eigenvalues of the covariance matrix R denoted {λ k } d k=1 , the expansion in Equation 1 can also be written as in Equation 4, where {φ k (.)} d k=1 together with {a k (.)} d k=1 are bi-orthonormal, and v k (.) = a k (.) By defining the matrix When n << m, it is more computationally efficient to solve the eigenproblem of R T instead of the eigenproblem of R as highlighted by Sirovich [96] . This is often the case when a two-dimensional physical field is measured over a domain at specific times, and is the case encountered for our application described in Section 3.
When an order d << min(m, n) corresponds to a high EVR as defined in Equation 2, we speak of dimensionality reduction, because the data are projected in a sub-space that is of much smaller dimension than R m×n . When diverse enough records are available for the variable under study, we may consider that the resulting POD basis is a generator of all possible states. Predicting the associated temporal coefficients {a k (t)} d k=1 at a given time t would therefore be enough to predict the whole state at time t. Hence, we propose to use the POD as a spatial basis extractor. This would first enable study of the spatial dynamics of the variable of interest and eventually extraction of physical information, as shown in the application Section 3. Then, the basis can be used as a generator for the prediction of future states. This implies predicting the evolution of {a k (t)} d k=1 , for which we propose to use Polynomial Chaos Expansion (PCE), as described in the following Section 2.2.

Polynomial Chaos Expansion
A reminder of the theoretical base of PCE is presented in Subsection 2.2.1. Theoretical details, demonstrations and interesting references can be found in [6,57,58,99,114]. After the theoretical introduction, a simple indicator is proposed in Subsection 2.2.2 for the analysis of the variables influence on the output value. The latter is later generalized for POD-PCE in Section 2.3.

Learning
The idea behind Polynomial Chaos Expansion (PCE) is to formulate an explicit model that links a variable of interest (output) to conditioning parameters (inputs), both in a probability space. This enables the propagation path of probabilistic information (uncertainties, occurrence frequencies) to be mapped from the input space to the output space. The variable of interest, Y, and the input parameters Θ = (θ 1 , θ 2 , ..., θ V ) are therefore considered random variables, characterized by a given Probability Density Function (PDF) denoted f Θ . It should be kept in mind that the outputs of our problem are the temporal coefficients Y = [a k (t)] k∈{1,...,d} generated by POD, and that among the inputs there can be a set of physical forcings, as described later in Section 3. The objective is to derive the evolution of the temporal coefficients as the outcome of the forcings. Let us now recall some fundamentals of the mathematical probabilistic framework, taking the example of a one dimensional real-valued variable. The definitions can be easily extended to R M .
Let (Ω, F, P) be a probability space, where Ω is the event space (space of all the possible events ω) equipped with σ-algebra F (some events of Ω) and its probability measure P (likelihood of a given event occurrence). A random variable defines an application Y (ω) : Returning to the PCE construction, inputs Θ = (θ 1 , θ 2 , ..., θ V ) are considered to live in the space of real random variables with finite second moments (and finite variances). This space is denoted by L 2 R (Ω, F, P; R) and is a Hilbert space equipped with the inner product (θ 1 , θ 2 ) L 2 R := E[θ 1 θ 2 ] = Ω θ 1 (ω)θ 2 (ω)dP(ω) and its induced norm The PCE objective is to map the output space from the input space with a model M as in Equation 5: where M 0 is the expectation of Y and M I⊆{1,...,V } represents the common contribution of the variables I ⊆ {1, ..., V } to the variation in Y . For the PCE model, these contributions have a polynomial form. We shall define, for each input variable θ i , an orthonormal univariate polynomial basis ξ where p ∈ N is a chosen polynomial degree and ξ (i) β (.) is of degree β. The orthonormality is defined with respect to the inner product (., .) L 2 R . If we introduce the multi-index notation α = (α 1 , ..., α V ) ∈ N V so that |α| = V i=1 α i , we can define a multivariate basis ζ Θ α (.), |α| ∈ [|0, p|] as ζ Θ α (θ 1 , θ 2 , ..., θ V ) : . Therefore, the model in Equation 5 can be written as: where c α ∈ R are deterministic coefficients that can be estimated thanks to different methods. It can be formulated as a minimization problem, and regularization methods can be used when dealing with small data sets. In the present study, we used the Least Angle Regression Stagewise method (LARS) in order to construct an adaptive sparse PCE. In this approach, a collection of possible PCE, ordered by sparsity, is provided and an optimum can be chosen with an accuracy estimate. It was performed in this study using corrected leave-one-out error [7]. The reader can refer to the work of Blatman [6] for further details on LARS and more generally on sparse constructions.
The choice of the basis is crucial and is directly related to the choice of input variable marginals, via the inner product (., .) L 2 R . Chaos polynomials were first introduced in [109] for input variables characterized by Gaussian distributions. The orthonormal basis with respect to this marginal is the Hermite polynomials family. Later, other Askey scheme hypergeometric polynomial families were associated to some well-known parametric distributions [113]. For example, the Legendre family is orthonormal with respect to the Uniform marginals. This is called gP C (generalized Polynomial Chaos) when variables of different PDFs are used as inputs. In practice however, the input distributions of physical variables can be different from usual parametric marginals. In such cases, the marginals can be inferred by empirical methods such as the Kernel Smoother (see [35] for theoretical elements). In this case, an orthonormal polynomial basis with respect to arbitrary marginals can be built with a Gram-Schmidt orthonormalization process as in [111] or via the Stieltjes three-term recurrence procedure as in [107].
To highlight the importance of the marginals and choice of polynomial basis for the learning process, several configurations are attempted in Section 3. Different input sets and distributions (Gaussian, Uniform, inferred by Kernel Smoothing) were tested. The influence of the polynomial basis on the ML is investigated in Section 3.2.2.

Physical importance measures
Once the PCE construction is achieved, a physical interpretation can be performed. It is notable that classical NN indicators can be used [25]. The PCE can be represented in the Feedforward NN paradigm as in Figure 2. Such networks are classically composed, in addition to the input and output layers, of successive hidden layers. Each hidden layer is composed of neurons that transform the variables of the previous layer (outputs of the previous neurons) into a new set of variables. This is done by combining a linear transformation, giving different weights to the previous neurons, and a transformation function, called Activation Function (AF). This succession of layers is called the latent representation. For a number of hidden layers L ≥ 1, the latent representation, in addition to the input and output layers, can be formally written as L|] and {f k } k∈[|1,L|] are the hidden layer weight matrices and AFs, A in is the input-to-hidden connection matrix and f out is the final hidden-output transformation [59,95].
The PCE-based NN represented in figure 2 is a single layer feedforward, composed of l ∈ N neurons, that can be written as Y ≈ f out (A 1 f 1 (A in Θ)). The first matrix A in is the input-to-hidden connection matrix of dimension V × V , that links the input layer to the PCE hidden layer containing the multivariate polynomials ζ Θ α , α ∈ {α 1 , ..., α l } , where V is the number of inputs and the multivariate indexes {α 1 , ..., α l } are conditioned by the chosen polynomial degree p such as ∀i ∈ [|1, l|] 0 ≤ |α i | ≤ p, and by the number of selected features if a sparse polynomial is constructed, as in the present case using LARS [6].
Matrix A in represents the contributions of the V variables to the multivariate polynomials The final hidden-output transformation f out is a summation. Figure 2 can also be presented differently: another hidden layer can be added to the PCE latent representation as Y ≈ f out (A 2 f 2 (A 1 f 1 (A in Θ)). The first AF f 1 would represent a transformation of each input variable to a list of monomials of degrees 1 to p (here, A in is identity). The second AF f 2 therefore represents the tensor product that transforms the different monomials to multivariate features, with A 1 appropriately filled with zeros and ones, and To capture the importance of each feature, the Garson relative Weights (GW) defined in Equation 7 are a classical measure to quantify the relative importance of each neuron of the last hidden layer, and therefore of each polynomial pattern, for the output value [25,39,104].
This measure can be used to understand the importance given by the NN algorithm to the variables and their possible interactions, especially when using feature selection algorithms as LARS: "feature interactions [...] are created at hidden units with nonlinear activation functions, and the influences of the interactions are propagated layer-by-layer to the final output" [104]. In the particular case of a polynomial expansion, the interpretation is straightforward, the importance of each variable alone corresponds to its monomials, and the importance of its interactions with other variables corresponds to the multivariate polynomials in which it is involved.
For the particular case of the orthonormal basis provided by PCE, the GW defined in 7 can be interpreted in terms of Pearson's correlations between output Y and the polynomial basis elements thanks to the orthonormality of the basis with respect to the scalar product (. , .) L 2 R that guarantees: Therefore, the weights w ζ Θ α can also be computed as This means that they measure the relative importance of the basis element in the expansion of the output, in terms of linear correlation, regardless of the sign of the latter. These "relative Pearson's correlations" can be seen as a physical contribution since the PCE model is strictly linear.
The sum of the GW w ζ Θ α for all the polynomial features equals 1. This means that they allow {ζ α } |α|≤p to be ranked in terms of relative contribution to the output Y . The contributions can be analyzed either for each polynomial pattern separately, or for a single variable θ i by adding all the polynomial shares related to this variable alone (1st Sobol indice analogy [99]) or by adding all the polynomial shares related to this variable and its interactions (total Sobol indice analogy [99]).

POD-PCE based predictor
POD and PCE were introduced separately in Subsections 2.1 and 2.2 respectively. We are now fully equipped with the adequate theoretical basis and mathematical notations, to present the POD-PCE ML methodology for a data-based model learning of a physical spatiotemporal field. In this Subsection, we will first summarize the proposed approach, then the formal details of the coupling will be given with the definition of adequate accuracy measures. Finally the previously discussed importance measures will be generalized for the POD-PCE physical study.
The proposed POD-PCE ML consists of five steps, in a learning and a prediction phase, summed up as follows: • Learning phase: * POD basis construction: given a set of measurements stored in the snapshot matrix U(X , T ) = [u(x i , t j )] i,j ∈ R m×n , construct a spatial POD basis accordingly ; * PCE learning: construct PCE models that map each POD temporal coefficient, obtained in the previous step along with the spatial basis, from time t j ∈ T to time t j+1 ∈ T , using the realizations of associated input variables denoted Θ(t j → t j+1 ) on the set T ; • Prediction phase: * POD projection: given a new measurement of the physical field U(X , t k ), obtain the values of the POD temporal coefficients using appropriate projection ; * PCE prediction: given the constructed PCE models and an estimate of the inputs from current time t k to future time t k+1 denoted Θ(t k → t k+1 ), calculate an estimation of the future POD temporal coefficients ; * POD-PCE ML prediction: reconstruct the estimate of the future state U(X , t k+1 ) by calculating the expansion on the POD basis, using the estimate obtained for the future temporal coefficients.

Machine learning methodology
Here, the formal hypothesis behind the POD-PCE ML reasoning and its mathematical formulation are discussed.
) be a vector of the inputs supposed to condition the evolution of U(X , .) over time. The dynamic model, denoted H, that gives an estimation of a future state U(X , t j+1 ) from a past state U(X , t j ) and an estimation of Θ(t j → t j+1 ) over the time interval [t j , t j+1 ], where t j < t j+1 ∈ R + , is formulated as in Equation 9 .
If the field of interest has been recorded over a set of past times T = {t 1 , . . . , t n } ⊂ R + , where t n < t j < t j+1 , a POD basis can be constructed as in Section 2.1, consisting of d ∈ N vectors of dimension m stored in a matrix as d (X )) ∈ R m×n , and can be seen as a generator of all possible states if enough records are available. If so, any future state U(X , t j ) can be expanded on this POD basis and the associated temporal coefficients are simply the weights of U(X , t j ) on the POD basis. They are therefore obtained using the canonical scalar product over R m , as in Equation 10.
Hence, the variable part of U(X , t j ) is fully expressed in the temporal coefficients a k (t j ). The field of interest U(X , t j ) can be either a field measurement, a laboratory or a numerical experiment. In any-case, it can be considered as being generated by a random process "in the sense that nature happens without consideration of what could be the best realizations for the learning algorithm" [95]. Therefore, the coefficients a k (t j ) can also be seen as the j th realization of a random variable A k . We can therefore construct a PCE approximation H k that maps A k from its input space. The latter is taken as a collection of random variables, composed from the set (A 1 , ..., A d ) at a previous time, the duration of the dynamic, and the input variables Θ(t j → t j+1 ). This is formulated as a classical dynamic model in Equation 11 .
The model H in Equation 9 is approximated as in Equation 12 .
The choice of input variables for regression models is an ongoing research question in statistics [33,81]. The chosen input variables can be extracted from a transformed version of a larger input set with the help of PCA [44] for DR. However, this approach was not studied here and will be the topic of a future study. Different input configurations will be evaluated, to investigate the influence of variable selection on the proposed learning. For example, the hypothesis of dependence between the random variables (A 1 , . . . , A d ) could be relaxed. This would imply writing the approximation in Equation 11 in a relaxed form as In that case a simpler model H, under the strong independence assumption, can be formulated as in Equation 13.
Both alternatives are tested in Section 3. To investigate the influence of input selection on learning accuracy, a quantitative evaluation of the hypothesis is needed. More generally, whether for the above-mentioned simplifications or for the approximated form of the model in general, accuracy estimators are needed. These are presented below.

Accuracy tests for the approximation
There are two determining parts in the POD-PCE learning process. Firstly, the PCE learning H k (.) of each mode A k should be as accurate as possible. Secondly, the reconstructed field k (X ) for a given rank d should be as close to the real field U(X ) as possible.
The distance between each mode and its PCE approximate can be evaluated using the generalization error, denoted δ(A k , H k ) and defined as in Equation 14.
For the model defined in Equation 13 , this error can be estimated, on a set of paired realizations (a k (t 1 ), . . . , a j (t n )) and (Θ(t 0 → t 1 ), . . . , Θ(t n−1 → t n )) , as in Equation 15 as explained by Blatman [6] . This approximated version of the generalization error is called the empirical error.
Its relative estimate denoted emp (A k , H k ) can be defined as in Equation 16 .
Once the PCE learnings can be trusted, the distance at time t j between the true state U(X , t j ) and the POD- can be defined. It might be estimated using the relative Root Mean Squared Error (relative RMSE), denoted r[U, H](t j ) and calculated as in Equation 17 , where h(x i , t j ) refers to the value of the POD-PCE approximation at coordinate x i and time t j .
A mean value of the relative RMSE is calculated over a set of realizations corresponding to a set of times T = {t 1 , . . . , t n }. It is denoted r[U, H] (T ) and estimated as in Equation 18.
Once the accuracies of the PCE learnings and the POD-PCE coupling have been evaluated, a final model, which will be the most accurate one, can be chosen. This model would, for our ML set-up, be the best representation of the dependence structure between inputs and outputs. It is used to shed light on the underlying physical relationships. Therefore the inputs are ranked in terms of physical influence, using an appropriate ranking indicator, presented in the following Subsection.

Physical influence of inputs based on the POD-PCE model
The GW influence measures presented for the PCE models in Subsection 2.2 are here extended for the POD-PCE coupling. These indicators are adequate for the analysis of each PCE model H k : i.e., for interpreting the contribution of the inputs to each random variable A k separately. However, calculating the contributions to each A k independently precludes putting them in perspective according to the importance of A k in the final reconstructed model H that approximates U(X , .). Hence, adapted indicators should be calculated.
Let U(X , .) be the random spatiotemporal field approximated by the POD-PCE ML, for prediction from time t j to time t j+1 and let H k be the PCE approximation at degree p (k) that maps the random POD temporal coefficient A k from a set of input variables, using the expansion on the multivariate polynomial basis ζ Equation 12 is written as in Equation 19: Thanks to its linearity, the POD-PCE ML can be represented as a single-layered NN, as shown in Figure 3. .

PCE 1 input layer POD -PCE Output layer
These GGW indicators show that the contribution of the polynomials {ζ Equation 2). An analogy can be drawn with the generalization of Sobol indices for a reduced order model [24,52]. The property d k=1 ..,d} to be ranked altogether in terms of contribution to output U. The contributions can be analyzed either for each polynomial pattern separately, or for a single variable θ i by adding all the polynomial shares related to this variable alone or by adding all the polynomial shares related to this variable and its interactions (analogy with first and total Sobol indices respectively [99]).

Application to a cooling water intake in a coastal zone
This section deals with the application of the POD-PCE ML described in Section 2 to a physical problem, introduced, with an industrial study case and inherent challenges, in Subsection 3.1. The physics and data are described. Subsection 3.2 deals with application of the POD-PCE learning phase to the data and assessment of accuracy and robustness with respect to the numerical choices (data set, inputs, marginals and polynomial basis). Finally, the prediction phase using POD-PCE is dealt with in Subsection 3.3, and the ability of the proposed ML to predict mean quantities and spatial details is demonstrated.

Study case
Sedimentation processes in nearshore areas can be responsible for the excessive sediment deposition commonly observed in cooling water intakes in power plants. As a result, the carrying capacity of the water intake can be drastically reduced, by decreasing its effective area of transport [98].
Cooling water intakes usually incorporate jetties, of which the angle with the shoreline and position relative to the direction of the net longshore sediment transport influence the amount of sediments diverted into the channel inlet by waves and tidal currents. Jetties also reduces littoral drift, resulting in localized sediment accretion against the shore-normal structure due to the longshore sediment transport being trapped by the jetty [19]. In addition, a return current is prone to develop, in the form of a swirling vortex at the end of the structure, and can induce sediment deposition in the vicinity of the channel entrance, consequently affecting the amount sediment delivered into the cooling water intake [15]. Therefore, effective water intake management involves frequent dredging, with high operational costs and usually hindered by a tight schedule. It is consequently necessary to assess intake sedimentation under different natural forcing and plant operation scenarios in order to optimize dredging operations to help mitigate the potentially adverse impact of the waves and tidal currents and meteorological forcing combined with plant functioning.

Site characteristics
The study site is located on the eastern English Channel coast in northern France. Tide in the study zone is classified as mega-tidal and is dominated by semi-diurnal circulation, with low-tide water depth of 10 − 15 m, and a mean tidal range of approximately 8.5 m, reaching 10 m during the spring tide [56]. Hydrodynamics are influenced by asymmetrical current velocities, with flood and ebb currents in the E-NE and W-SW directions, respectively. Current velocity at 2.2 m above seabed vary from 0.70 to 0.98 m/s, depending on flood/ebb phase, respectively [68]. Wave activity in this open exposed environment is moderate, with significant annual and decennial wave height of 3.8 m and 4.7 m, respectively, with maximum values of 4.2 − 5.8 m, averaged period of 7 − 9 s and a predominant W direction. Orbital velocities measured during the spring-tide period ranges between 0.5-1.3 m/s. An example of tidal levels, wind direction and velocity and wave height and direction in January 2016 is shown in Figure 4 . In the study area, bed sediment varies from medium to fine silted sands, with a morphology characterized by the presence of mega-rides parallel to the coast. In this zone, rock occupies less than 4% of bed surface [15].

Data
Hydrodynamic and meteorological information comprise wave height, period and direction and wind velocity and direction, provided by the VAG prediction model of the sea state [31], using retrospective 3-hourly simulations between 2009 and 2018. Tidal water levels were obtained from the SHOM-REFMAR tide gauge station located in the vicinity of the study zone, with hourly survey frequency [85].
Bathymetric measurements were available from Single-Beam Echo Sounding on 39 cross-sectional profiles of intake measured at 25 m intervals, collected fortnightly between 2005 and 2018. Mean profiles were 100 m long with 0.5 m spatial resolution of bathymetric data. Additional information such as the daily coolant flow rates, and channel dredging volumes and frequency, were provided by the plant operator. A summary of the data is shown in Appendix A.
The available measurements of the forcings did not have the same frequencies. One solution to homogenize frequencies consists in reducing the measured data to representative statistics over the sedimentation interval ∆t ≈ 15 days separating two bed elevations measurements. Hence, the following statistics were used: • Tidal level indicators: average low tide (T Lmean), minimum low tide (T Lmin), maximum tidal range (T Lrange) and standard deviation (T Lstd) ; • Wind indicators: average wind velocity (W mean) and average direction weighted by velocity (W dir) ; • Wave indicators: average wave height (W vH), standard deviation (W vstd), average wave period and average wave direction weighted by height (resp. W vper and W vdir), average wave height exceeding the 90 th percentile (arbitrary storm indicator, W v2m) and percentage occurence (W v2m%) ; • Operational indicators: average pumping flowrate (Qmean) ; time lapse since last dredging (Dp), and last dredged volume (Dv).
These statistical indicators were calculated for each sedimentation interval ; a resulting scatter plot is shown in Figure 23. There was, for example, a positive correlation between mean low tide T Lmean over the studied sedimentation periods and wave parameters W vper and W vH. Mean wave periods W vper and mean wave heights W vH were also positively correlated.
For the learning part, the data overlapped only over a limited period. A maximum of 60 measurements could therefore be used, with up to 15 forcing variables. Obviously, this "small data" configuration is a considerable handicap for the dimension of the problem, especially given that the variable of interest is a two-dimensional bathymetric field. However, permanent intake monitoring ensures that the data set will always grow and can be used to update the learning. This limitation shall not prevent testing the accuracy of the methodology on small sets such as are often encountered in physical applications, as attempted below, where learning and prediction using POD-PCE is applied to the described data. For the learning algorithm, input variables are needed, corresponding to the reduced statistical indicators described above, and denoted (θ 1 , ..., θ V ), where V is the supposed dimension of the problem.

Measurement-based learning of a physical field using POD and PCE
This section concerns learning the spatio-temporal bathymetric field using POD and PCE independently. The modes are extracted in Subsection 3.2.1 using POD and the temporal patterns are learned as a function of the forcing parameters using PCE in Subsection 3.2.2. Throughout this investigation, particular attention is given to the convergence of the learning and to its robustness with respect to the numerical choices. Trusted POD-PCE learning is immediately used for physical interpretation and the most important physical insights resulting from it are summarized in Subsection 3.2.3.

Physical analysis and data reduction using POD
First, POD was applied on the bathymetry measurements. The aim was to identify morphodynamic patterns so as to better understand the sediment deposition inside the channel, and to characterize variations in depositions with the external forcing variables. After setting aside poor-quality measurements (e.g. incomplete bathymetries), a total n = 156 realizations were used. The bathymetry points were sonar boat measurements on m p = 39 crosssections inside the intake. Linear interpolation was performed on m i = 100 fixed points for each profile, in order to express all measurements on the same grid, giving a total m = m i × m p = 3, 900 spatial points. The interpolated realizations were then stored in a snapshot matrix Z(X , T ) = [z(x i , t j )] i,j ∈ R m×n and POD-processed as explained in Section 2.1. The EVR defined in Equation 2 and the mean relative RMSE between the POD approximation and the complete measurement (averaged over the realization set as in Equation 18 ) were calculated for each POD approximation rank and are plotted in Figure 5. The first pattern represents over 94% of the variance, and explains most of the variation in dynamics. The variance percentage reached 99% at rank 14, where the mean error was slightly over 10%, decreasing to 8% at rank 20. Dimensionality reduction is therefore a realistic option for this specific dynamic problem. This encouraged the learning and prediction attempts undertaken in Subsections 3.2.2 and 3.3 respectively. The spatial and temporal components of the first four POD modes corresponding to an EVR higher than 97% are respectively plotted in Figures 6 and 7. Figure 6: The first four spatial patterns of the POD applied to intake bathymetries.
(a) Coefficients a 1 (t) and a 2 (t) (b) Coefficients a 3 (t) and a 3 (t) Figure 7: The first four temporal coefficients of the POD applied to intake bathymetries.
The first spatial pattern (Figure 6-a) represents the channel's slope. Its temporal coefficient (Figure 7-a) shows regularity in time that is almost periodicity. When it increased, overall sediment deposition in the channel increased, because the difference between the upstream and the downstream bed elevations, and therefore the slope, diminished. The sediment deposition in the channel might be related to the increasing sediment supply caused by the external forcing influence. Decrease always corresponded to a dredging episode. The apparent periodicity is therefore not natural or seasonal but due to periodicity of operational intervention: sediment deposition in the channel is tolerated up to a certain level and then dredging is always undertaken at a certain point, which corresponds to the maximum of the temporal coefficient. The second pattern (Figure 6-b) acts as a geometric distribution function of the sediment deposition. In general, when the first temporal coefficient was maximal, the second coefficient was positive, meaning that the sedimentation mainly occured in the middle of the first portion of the channel (upstream), on the right bank of the bend and on the left bank in front of the pumps. This spatial distribution can be associated to the internal flow characterized by a velocity distribution inside the channel. In fact, the sediments settle where velocity is the lowest, which is probably the case where the banks appear. The computed sediment deposition and erosional patterns are analogous to those commonly observed in meandering rivers [42].
The third pattern (Figure 6-c) shows sediment deposition concentrated in the first portion of the intake, and the fourth pattern ( Figure 6-c) emphasizes sediment dynamics, particularly in the downstream part of the channel. This behavior is statistical proof and quantification of finer sediment supply. The finer sediment fraction was transported in suspension and deposited at the end of the intake channel. The temporal coefficients associated with the third and fourth mode (Figure 7-b) were less regular than those of the first and second mode, and seemed to follow a more stochastic dynamic. The peaks may represent unusual sediment supply, probably linked to extreme events (extreme tides, storms, etc.).
To check the robustness of the statistical conclusions deduced from POD, convergence analysis is necessary. This was performed on the EVR values associated with the first four patterns, using bootstrap analysis [20]. The results are shown in Figure 8. The convergence of the mean values and the tightening of the confidence intervals around the mean whith increasing matrix size are clear for these first four modes. However, whereas the confidence intervals represent at most an error of ±0.6% around the mean for the first mode, they reached respectively ±12.5%, ±25% and ±12.5% for the second, third and fourth modes. The analysis proved that the POD results could be used to pursue the learning. Firstly, a high EVR and low RMSE were associated with a small number of modes, guaranteeing optimal data reduction (d << min(m, n) as explained in Section 2.1). The number of POD modes to accurately represent the bathymetry can be chosen accordingly. In the present study, the configuration was d = 11 modes (discussed in Step 3 ), guaranteeing EVR ≥ 98% and information loss ≤ 12% (mean relative RMSE). Secondly, the EVRs were guaranteed to converge statistically at least for the first four modes, with error of ±0.6% around the mean for the most important mode, representing over 94% of the variance. Thirdly, the deduced patterns were physically coherent. Lastly, more than a decade of evolution was used to extract the POD basis, under variable operational and environmental conditions. As long as the operating conditions of the intake remained unchanged, it can be assumed that a wide range of evolutions has been covered, except for extreme events that rarely occur and that are not specifically treated in this study [26]. Hence, the POD basis can be considered as a physically trustworthy and mathematically complete basis to understand past evolutions and to predict future ones. The learning of temporal coefficients is therefore attempted in Subsection 3.2.2.

Learning of the POD patterns using PCE
The temporal coefficients calculated with the POD in Section 3.2.1 (Figure 7) were learned using PCE (theory in Section 2.2). The aim was to learn the way these coefficients evolve over time, as a function of the forcing parameters presented in Table 2 , with the ultimate objective of field prediction as explained in Section 2.3 and applied in Section 3.3. The present section focuses strictly on the learning phase and the physical analysis of the learned model, highlighting quality of learning (robustness, convergence, etc.).
The investigation of learning is organized in four steps.
• Step 1 -Sensitivity of learning to inputs and marginals: different configurations were tested to practically demonstrate the implications of these choices on the accuracy of fit. • Step 2 -Convergence and Robustness of fit: The best model resulting from Step 1 was studied more deeply.
Its convergence and robustness with respect to the choice of training members are were analyzed. • Step 3 -Physical interpretation of the best learned model: the best model was chosen, and the most influential forcings were ranked using the Evolution Weights (EW) and Generalized Evolution Weights (GEW) presented in Sections 2.2 and 2.3 respectively. • Step 4 -Robustness of the physical interpretation with respect to the learning-set members: the physical conclusions of the model were shown to be statistically meaningful.
These steps, in the aboved order, follow the logic of statistical model construction to build a trustworthy prediction algorithm, used in Section 3.3.
Step 1 -Sensitivity of learning to inputs and marginals Input variable selection is capital, and marginals must be chosen wisely. Below, we demonstrate the influence of these choices on the performance of the learning. Different configurations were tested.
Before introducing the tested configurations, the training steps that are common to all configuration need to be defined. A learning data-set is classically separated into different sub-sets, corresponding to different steps of the learning algorithm. This is commonly referred to as the "Train-Validation-Test Split" [95]: a training set is used for the learning, a validation set is used to check the learning and for further calibration, and a test set is used to assess the prediction capability of the statistical model. However, the data-set used in this study was small: the bathymetry and forcings measurements shown in Subsection 3.1 overlap for the 2012-2017 period only, leaving 64 sedimentation periods to study. Therefore, only a "Train-Predict" split was performed, where the prediction set played the role of both the test set and validation set. Hence, the numerical choices were calibrated on the training set, and validated on the prediction set for both statistical accuracy and physical prediction. The learning was then performed with an arbitrary choice of training-set size at 50, which left a prediction set of 14. The training data were chosen in chronological order (first 50 records), to mimic the learning process in an industrial context. This arbitrary training-set choice had consequences for learning ; the sensitivity of learning to training set choice was investigated (Step 2 ). All the model choices presented below (choice of inputs and marginals) were assessed on this training configuration. To assure that comparison is made between models at their best performances, the PCE polynomial degree was optimized for each separately. Degrees from 1 to 7 were tested, and the associated relative empirical errors on the training and prediction sets, respectively T and P , were calculated as in Equation 16 . The PCE degree that minimized the training and prediction errors for each model was chosen ; the corresponding result is referred to as "optimal" learning.
Three different input configurations were used for the learning of each temporal coefficient a i as generally formulated in Equation 11 .
• H i -model: a first simple configuration where all the statistical indicators described in Subsection 3.1 were used and an independence hypothesis between the POD temporal coefficients is considered. The model is written as in Equation 13 : . This is a model of dimension 17. • H F i -model: a more complex configuration where a "Full" 15-mode POD approximation is considered with possible dependencies between the temporal coefficients. Of course, the choice of the basis size and the dependency structure can be optimized, but the objective here was to make a first step toward a more optimal configuration. The model can be written as: . This is a model of dimension 31.
• H P i -model: a smaller set of inputs, used by the operators to qualitatively evaluate sediment deposition risk, was used. It corresponds to the six variables T Lmean, W vH, W vper, W vdir, W v2m and W v2m%. This mimics the physical expertise that may be engaged when building a statistical model. It is written as , where Θ P stands for the "physical". This is a model of dimension 8.
To these variable choices were associated three choices of marginals, conditioning the choice of the PCE orthonormal polynomial (Section 2.2).
• Lgd: all the variables follow a Uniform PDF. The bounds of the marginal were set to the minimum and maximum chronological values ±1% as in [103]. The associated orthonormal polynomial basis is the Legendre family. • Hrm: all the variables have Gaussian marginals characterized by the empirical mean and variance deduced from the data. The associated orthonormal polynomial basis is the Hermite family. • Stlj: the marginals were inferred from the data using Gaussian Kernel density estimates. The orthonormal polynomial basis was constructed from the knowledge of the marginal using a Stieltjes orthogonalization.
The three marginal choices (Lgd, Hrm, Stlj) were trained with the three dimension choices (H P i dim=8, H i dim=17, H F i dim=31). The empirical errors of the "optimal" learnings are compared in Figure 9 .  For Modes 1 and 2, training and prediction errors were almost identical for all configurations, although with a slight advantage with the smallest dimensions for all the marginal types in the learning of Mode 2. Starting from Mode 3, bigger differences emerged. At the learning step of Mode 3, models of dimension 17 and 31 were poorly fitted for the Stlj and Lgd configurations compared to others. At the prediction step of Mode 3, the errors of models with dimension 31 were much greater than smaller dimensions for all marginal choices. There seemed to be an overfitting of the model by selecting a larger number of inputs. The best models for Mode 3 were those of the smallest dimension, 8, with either the Stlj or the Hrm model. Lastly, for Mode 4, two orderings were observed for the prediction error. Firstly, for each marginal choice, prediction error increased with dimension, which confirmed the overfitting hypothesis. Secondly, error was the smallest with the Stlj model (Kernel density), followed by the Hrm model (Gaussian) and lastly by the Lgd model (Uniform). Here, Uniform marginals performed worst ; they were probably too different from the real data marginals and did not account for particularities in the inputs. In the parametric family, Gaussian marginals probably fitted real density better.
To conclude this comparison, the best marginal choice was the Kernel density estimate. The smallest dimensions performed the best, with the Stlj choice for the polynomial basis. The H P i ; Stlj model of dimension 8 was therefore selected. However, the training was performed with an arbitrary split of the available statistical set. The sensitivity of the model to the learning set size and members is performed in the following step.
Step 2 -Convergence and robustness of the fit Up to this point, an arbitrary number of 50 measurements was used for the training phase, leaving 14 prediction points for testing purposes. In the following, the influence of training set size on the learning and prediction error is assessed. The objective is to check the robustness of the previous best model H P i ; Stlj with respect to the data-set size and members. The evolution of the training and prediction empirical errors according to training set size is shown in Figure 10 and 11 respectively. For each training set size, members were chosen randomly among the full data-set, and the remaining members were used for the prediction phase. For the estimation of the confidence intervals, bootstrap analysis was performed [20]. For comparison, the convergence of the H P i ; Hrm model is plotted for Mode 3 and can be found in Appendix B for Modes 1 and 2.
For the first two modes, the training errors in Figure 10 show a convergence of the median value and a tightening of the confidence intervals. The trainings can be considered as converging from around training size 40. The associated median prediction errors in Figure 11 globally decreased with increasing training set size. However, although the final median values were lower for the Hrm model, the confidence intervals were much larger than for the Stlj model. The latter seems much more robust with respect to changes in training scenario.  Figure 12 and Figure 13 for training and prediction, respectively. Only the middle 80% portion of the residuals range is plotted, in order to analyze the center of the distribution ; the full residuals distribution was long tailed, because the confidence intervals associated with small training set sizes were too large and produced extreme behaviors of the model. The training residuals were generally centered around zero: i.e., the models are unbiased. A slight asymmetry was, however, observed for Mode 1, which means that a 1 (.) was more often overestimated by H P i ; Stlj. Consequently, the mean elevation in the channel and the mean global sedimentation may be slightly exaggerated. These exaggerations, however, remained within a reasonable range, as most of the residuals fell within the ±2% interval. The training residuals of Modes 2 and 3 were perfectly centered, but percentage error dramatically increased. Most of the residuals fell within the ±10% interval for Mode 2, whereas they reached ±50% for Mode 3. However, this error concerns modes that represent at most 4% of the total bathymetry variance, as more than 96% of the total variance was already captured by the addition of the first two modes. The residuals shapes (i.e. slight overestimation for Mode 1 and perfect centering for Modes 2 and 3) were maintained through the prediction phase. Furthermore, the residuals mostly fell within the ranges identified in the training phase. H P i ; Stlj model behavior was stable. The prediction uncertainty could therefore be measured and trusted and the physical interpretation was consequently robust, as discussed below in Step 3.
Step 3 -Physical interpretation of the best learned model The calibrated H P i ; Stlj PCE models were considered optimal, as they showed good fit, convergence and robustness with respect to the training choices. Here, they are analyzed to deduce physical information. Firstly, the optimal polynomial degrees selected for each mode and the associated training and prediction empirical errors are shown in Figure 14. Linear models were optimal for Modes 1 and 2 (degree 1), and the associated errors were low for both the training and the prediction sets. For Modes 3 and 4, the optimal polynomial degrees increased, which implies higher-order contributions and/or higher-order interactions for the input variables. For modes of higher ranks, the models were either linear (degree 1) or approximated by a simple average value (degree 0). This means that LARS rejects polynomial terms of higher degrees because they do not significantly improve the learning [6]. Prediction relative empirical errors increased from Mode 3, but remained under 50% up to mode 11. This percentage error seems surprisingly high, but must be interpreted according to the meaning of the relative empirical error, which, as calculated in Equation 16, is a measure of the missing variations (distance between the model and reality) relative to the variance of the data. It is also often called "the fraction of unexplained variance" [17], and represents the amount of variance that was not captured by the PCE model. For example, for mode 12, estimated with a degree 0 PCE model, 100% of the variance is not captured, which is natural because only the average value is accounted for with degree 0. For Mode 13, LARS did not include polynomial terms of degree higher than 1 (no improvement of the learning), but this appears to be a bad choice since the prediction error was much higher than the training error, meaning that the training set or inputs were not enough to guarantee accurate learning in Mode 13. For Modes 3 to 11 with error up 50%, this means that either the training set or the used inputs made it impossible to predict more than 50% of the variance. However, as presented in Step 2, this 50% error concerned at most 4% of the total bathymetry variance. Hence, the errors starting from Mode 3 represented at most 2% of missing variance. Beyond Mode 11, prediction with PE would not be optimal, as the prediction error dramatically increases.
Secondly, the explicit formulations of the PCE models (reported in Appendix B) were used to analyze the contribution of each forcing variable to the dynamics. for this, the Garson Weights (GW) defined in Equation 7 were used to estimate the influence of the forcings on each temporal coefficient. The global influence on the whole bathymetry field was quantified using the Generalized Garson Weights (GGW), as in Equation 20.
The ranking of the modes and the impact of the inputs is represented in Figure 15. The inner circle represents the contribution of each mode to the whole field. Mode 1 corresponds to a major contribution and the following modes are ranked in terms of contribution, which is the essence of POD. The outer circle represents the contribution of each polynomial term in the PCE representation. The share of each polynomial term corresponds to the GGW in relation to the global contribution (full circle). When this share is compared to the importance of the corresponding mode, it corresponds to the GW. Lastly, the polynomial terms corresponding to more than 0.5% GGW are indicated. ζ α=(.) (.) corresponds to the notation introduced in Subsection 2.1, with the multi-index notation for α that represents the polynomial degree of each monomial. For example, ζ α=(α1,α2) (θ 1 , θ 2 ) corresponds to a polynomial of degree α 1 + α 2 , where θ 1 contributes as a monomial of degree α 1 and θ 2 as a monomial of degree α 2 . The meaning of the variables that appear in Figure 15 can be found in Subsection 3.1 Figure 15: Piechart of the most influential parameters, using GW and GGW on the POD-PCE. The inner circle represents the share of each mode. The outer circle represents the share of each polynomial term. The polynomial terms corresponding to GGW higher than 0.5% are shown.
For all the temporal coefficients a i (.), the most influential contributor by far was the value of the previous state a i (t j−1 ), in the form of a monomial of degree 1. It is followed by contributions involving the mean wave height during the sedimentation period period W vH for all the modes, which makes W vH the most important external forcing, figuring in the third position among all the forcings, with a contribution of 9.6% through the first mode, a total of 12.6% if only W vH monomials are considered, and 13.3% if interactions with other variables are taken into account.
The other forcing contributions also appeared, but with much less importance: e.g., the influence of mean low tide level T Lmean, which took an interaction form with the previous bathymetry shape for Mode 3. Firstly, this interaction makes sense in terms of physics, as sediment deposition is conditioned by the value of bed shear stress [105], which depends on velocity and water depth. The water depth value is exactly the tidal level minus the bed elevation value, which here appears as a multiplicative interaction between T Lmean and Mode 3. Second, the value of this contribution was only 0.7% GGW, which is negligible when compared to the first contribution of W vH. The learned model gave much more importance to waves than to tides. This does not necessarily mean that tides have no influence on sediment deposition, but may simply suggest that, in the present configuration, sediment mobilization by the tide is always more or less the same, and that the forcing that makes a considerable difference is the variation in wave heights. Waves are a determining factor for sediment mobilization in coastal configurations, through the influence they have on bed shear stress [105]. Further more, a noticeable correlation between W vH and T Lmean can be seen in Figure 23 in Appendix A, which means that the information of low-tide levels is to a certain extent contained in the mean wave height. There is therefore a probable dependency between these variables. In case of dependencies, the iterative process used by LARS may drop a variable that is physically important because the important information is already contained in another variable, due to their dependency.
Lastly, it is important to note that the contribution of less frequent wave events was also present but to a much smaller extent. It is represented by a polynomial term in the form ζ α=(1,1)(W v2m,W v2m%) , where W v2m and W v2m% are respectively mean wave height exceeding 2m and the associated frequency of occurrence (arbitrary storm indicator chosen in Subsection 3.1). This term appears in Modes 3 and 4 for a maximum total influence of 0.3%. Higher-order interactions and less frequent events are therefore represented by modes of higher rank, associated with smaller variance percentages.
Step 4 -Robustness of the physical interpretation with respect to the learning-set members As a last proof of the robustness of the proposed learning algorithm, specifically concerning physical interpretation, a sensitivity analysis with respect to the training set members was performed. The robustness of the calculated Garson Weights (GW) with respect to the choice the training members was studied. This is equivalent to studying the robustness of the polynomial basis term selection as produced by LARS, and their associated multiplicative coefficients.
For this, a Bootstrap analysis was again used to construct different learning sets of size 50, instead of choosing the first 50 measurements. This produces a distribution of the GWs rather than a single value, for each polynomial term. The result is shown in Figure 16 for the weights of the a k (t j−1 ) and W vH monomials.
(a) Polynomial term a k (t j−1 ) (b) Polynomial term W vH Figure 16: Probability density functions of the GWs associated with the degree 1 monomials of variables a k (t j−1 ) and W vH. The training size is 50 using 20 different random picks for each size.
For modes 1 to 4, the median weights P50 (Percentile 50) of the a k (t j−1 ) monomials, represented in Figure  16-a, were always over 0.6, but the variation range was strictly less than 1, with density functions centered around the median and a small standard deviation for modes 1, 2 and 4. This means that whatever the training set, the previous state a k (t j−1 ) value was always predominant but never enough to estimate the evolution of the first four modes. A tendency (in particular linear) using the last state was not sufficient, and additional information was always needed (forcing). In parallel, Figure 16-b shows that this information is certainly the waves, as the median values of the GW for the first four modes were between 10 and 25%, corresponding to the information gap left by the previous value variable a k (t j−1 ) in the fitted PCE model.
Starting from mode 5, the median values of a k (t j−1 )'s GW had greater chance of falling around 1, which means that the associated polynomial models only rely on the last recorded value of the mode for the future guess. In other terms, the constructed model consists of a linearization around the previous value (tendency capturing) and does not incorporate the correlations between the future-state and the forcing variables (causality model). This can be explained by the small variances of the higher-rank modes and the difficulty of learning the PCE models from statistics averaged over the sedimentation periods. Additionally, the P25-P75 confidence interval moved to the upper bound of the density functions.

Summary of the physical insights from the learning
The spatial patterns as deduced by POD express the spatial correlation in the sediment deposition from the upstream to the downstream part of the channel. The EVR reached 99% with d = 20 modes only, where the mean relative RMSE between the approximation and reality was slightly over 10%. This is a statistical proof that the spatial correlations expressed in the POD patterns are explanatory of the physical dynamics over their whole range of variation (at least that observed from 2010 to 2018), with a low approximation rank. In conclusion, the dynamic problem exhibits fairly strong spatial correlations, and the solution to the problem can be expressed on a finite orthonormal basis.
The temporal patterns express the evolution of the sedimentation, as they multiply the spatial patterns. They were learned using PCE as a function of the previously cited inputs (previous states and forcings). The statistical model configuration (dimension and marginals) was chosen after an investigation of different options. The associated training and prediction error converged for the first three modes, and are characterized by tight confidence intervals. The residuals of the selected model were either negligible or centered around zero, demonstrating the unbiased character of the learning and prediction. The fitted models are of lower degree for the low-rank modes 1 and 2 and of higher degrees for modes 3 and 4, which are higher-rank, due to the emergence of interactions between the forcings, namely variables related to extreme behavior (storm events). The model mainly relies on the last state information, showing a strong correlation/continuity in time of the studied physics. Using GW, which measures the forcing influence for the first five modes, the action of waves was highlighted by the PCE model as a determining phenomenon. The first mode influenced the dynamic with a rate of 64.9%, the previous value of the second mode with a rate of 10.2% and, in third position, the mean wave height with a rate of 9.6%. The remaining 15.3% is essentially associated with previous values of higher order modes (10.3%), interactions with tides and contributions of other wave indicators. The GWs show robustness with respect to the choice of the training set members, which makes them trustworthy, at least for temporal correlation and analysis of wave influence. The main physical conclusions are that the dynamic problem is characterized by strong temporal correlations, representing more than 85% of the evolution, with an external sediment source, mainly represented by the waves, representing not more than 15%.

Prediction of a physical field using POD-PCE coupling
After performing both POD and PCE independently, the accuracy of a Machine Learning process using a POD-PCE coupling was assessed as in Section 2.3. In the continuity with Section 3.2, the first 50 historical bathymetries were used for training and the other 14 for forecasting. First, the impact of the size of the POD basis on the prediction process is assessed in Subsection 3.3.1. Then, the best size was determined and the average prediction behavior is analyzed in Subsection 3.3.2. The accuracy of the POD-PCE ML in predicting spatial details is assessed on cross-section examples in Subsection 3.3.3 and a summary is given in Subsection 3.3.4.

Influence of POD basis size
In order to track the errors generated by the various steps of the algorithm (POD, PCE and coupling), the mean relative RMSE (averaged over the prediction set, as in Equation 18 ) was calculated for each step and for each approximation rank d, as follows: • Reduction error: distance between the POD approximation d k=1 a k (t)Φ k (x) and the corresponding measured bathymetry field z(x, t) ; • Learning error : distance between the POD approximation d k=1 a k (t)Φ k (x) and the prediction using the POD-PCE coupling formulated as d k=1 H P i a i (t j ), t j+1 − t j , Θ P (t j → t j+1 ) Φ k (x) ; • Prediction error: the resulting final error between the prediction using POD-PCE coupling and the corresponding measured bathymetry z(x, t).
The results are shown in Figure 17 . Reduction error decreased from 16% to 3%, with increasing approximation rank. The error followed a logarithmic trend, with a significant slowdown from rank 9. These errors are coherent with the errors averaged over the full set (rather than the prediction set only) in the POD results Section 3.2.1 (around 8%). The learning error increased from 1% to 5% with increasing approximation rank, which is natural because the complexity of the model is increased. In fact, a prediction of rank d + 1 has an additional temporal coefficient that is predicted as compared to rank d. It is therefore natural that the distance between the approximation increased with increasing rank. The learning error order of magnitude was consistent with the empirical prediction error of 4% for mode 1 (as calculated in Section 3.2.2), associated with an EVR of over 94%. Lastly, the prediction error decrease is the balance of, on the one hand, the increase in accuracy by adding POD modes and, on the other hand, the increase in forecasting error with increasing number of temporal coefficients to be predicted. Consequently, the prediction error decreased from 16 to 6% up to rank 11, following almost the same decreasing trend as the reduction error. However, the decrease rate became slower and increasingly subdued, being overtaken by the learning errors, which dramatically increased starting from mode 12, as seen in Figure 14. Hence, a POD-PCE model of size 11 was selected for prediction.

Average performance of the chosen model
Average sediment deposition was predicted using the POD-PCE model of rank 11, for each of the 14 prediction dates. The average sedimentation rate, denoted S r , was calculated for time t j representing the sedimentation over [t j−1 , t j ], as in Equation 21. For operational estimation of sediment deposition, only the positive evolutions are of interest ; therefore, the erosion points were discarded in calculating rate S r by cancelling negative evolutions. Indeed, and therefore a null contribution to the sedimentation rate S r . Furthermore, only regions of considerable depth are of interest. Therefore only n p bathymetry points under −1 m (x i , i ∈ N p ) were taken into account. The results are shown in Figure 18 S r = 1 2n P i∈Np Figure 18: A comparison between the real sedimentation rates and the POD-PCE prediction of the sedimentation.
The POD-PCE prediction globally followed the real sedimentation trend, for example from Dates 1 to 8. When it was not equal to the real sedimentation rate, it was generally an overestimation, which is coherent with the asymmetry observed in the distribution of Mode 1 training and prediction residuals (Figure 12 in Subsection 3.2.2).
For the particular Date 11 however, half of the sedimentation was missing. Investigation of the data for this particular measurement showed that the previous record, taken as input, had been made 29 days previously, which is far from the average ∆t ≈ 15 days ; it is twice the mean interval, thus underestimating sedimentation by half. As measurement intervals were in general around the average, sedimentation time interval ∆t was not selected as a key parameter by LARS for the dynamic model, although it was given as an input and is physically significant. For the particular case of the time variable, multiplicative enhancement can be intended as a correction. However, this shows one of the limitations of statistical modeling: statistical significance can be confused with physical importance. Indeed, for the statistical conclusions to be physically significant, the measurements should be diverse enough to account for the variations in the inputs and the impact of these variations on the output. This was unfortunately not guaranteed for sedimentation measurement intervals, as they were often equal to 2 weeks. Additionally, for Dates 12 and 13, a large part of the wave measurements were missing in the sedimentation time interval. Consequently, mean wave height W vH was estimated over only a small portion of the time interval. This may lead to a good prediction (Date 13) if the interval used is representative enough of the full interval, and to bad prediction (Date 12) when not, and highlights the limitations of statistical averaging.

Spatial details of prediction by the chosen model
The spatial details of the prediction were analyzed on cross-sections for specific prediction dates. First, sediment deposition was observed on a cross-section at the entrance of the intake ( Figure 19). In accordance with the previous conclusions from Figure 18, the POD-PCE prediction captured various sedimentation ranges, as shown with Dates 2 and 7. However, a slight artificial sedimentation was predicted whereas there was no dynamics in reality, for example for Date 1 (Figure 19-a), due to the fact that the model is continuous, whereas threshold phenomena can occur in reality. Next, sediment deposition was observed on a cross-section at the middle of the first portion of the intake ( Figure  20). Mean sedimentation was well captured, but some details of the bathymetry were missing, such as formation of a new feature for Date 5 (distance 20 to 25 m) and Date 7 (distance 30 to 35 m). For Date 6, sediment deposition was slightly underestimated in the right bank and overestimated in the left bank. However, although the details of sediment deposition were not perfectly captured, the value of the sedimentation area corresponds well enough to reality. It can also be concluded that the way the RMSE and relative errors are averaged in space, for example in Figure 17, actually penalizes the accuracy of the algorithm because it does not take account of the oscillation of the prediction around an accurate mean. Then sediment deposition was observed on a cross-section at the bending part of the intake ( Figure 21). It shows that the prediction algorithm understood that the sediment deposition mainly occured in the right bank of the channel, for example for Date 7, even though it was overestimated. Furthermore, considering modes of higher rank from the previous measurement as an input, the algorithm captured the swing in the profile throughout its history, which is here observed from Date 4 to 14. Lastly, a cross-section of the last portion of the channel, in front of the downstream pumping station, is shown in Figure 22. Once again, the prediction algorithm understood where the sediment deposition occurs, this time in the left bank of the channel, coherent with the pattern represented by Mode 2 (Figure 6-b). However, it can be seen that unusual sediment deposition occurred for Date 11, which was not captured by the model, and may correspond to the arrival of less frequent fine sediment that appears in Mode 4 ( Figure 6-d). This also explains the sediment deposition error observed for Date 11 in Figure 18.

Summary of POD-PCE algorithm performance
Overall, the proposed learning algorithm showed interesting prediction characteristics. Firstly, model complexity can be increased gradually, by increasing the number of POD modes when accurate. Plotting error against the number of modes shows a convergence that helps in selecting the optimal number of modes. The average RMSE of the predicted field remains reasonably low. It was around 6% with the 11 modes selected in the present case.
Secondly, trends are well captured for spatially averaged quantities (here sedimentation rate) and for detailed spatial representations of the field. Good spatial distribution of evolution is guaranteed by the POD modes, even when evolution amplitude is over-or under-estimated.
Lastly, some disadvantages should be noted. For example, less frequent events that are represented by modes of higher ranks can be overlooked. Furthermore, sudden changes in features were not sufficiently captured, due to the high temporal correlation between last state and future state that was incorporated in the learning.

Summary and discussion
In this study, POD-PCE coupling for field-measurement based Machine Learning was proposed and assessed, in an industrial context in the field of geosciences, for complex physical phenomena involving non-linear dynamics and various forcings.
POD showed excellent performance both for reducing problem dimensionality and for physical interpretation of spatial and temporal modes. This is an important property of POD [47,48]. Even though the spatial patterns are mathematical modes that do not have explicit physical meaning, an interpretation could be given, at least for the patterns that were statistically stable, because they are statistical representations of the dynamic behavior. Our investigation showed that temporal signals are also relevant: they can show more or less regularity, which can be related either to the representation of different space or time scale physics, or to less frequent events. In our particular problem, the mean of the bathymetry field was not deduced from the field before applying POD, as can be the case for classical applications [14]. It was used as valuable information for decomposition, which allowed correlations to be studied between the first and second patterns, as well as the construction of the learning. The potential of POD for detecting biased and missing data was also assessed. POD was first applied to the whole set of measurements, but discontinuities emerged in the temporal signals of the decomposition. Such a procedure is important because, in most of cases, the data need to be filtered, which is a time-consuming task. The POD enabled fast recognition of elements that react differently from the average. However, many points of improvement are worth mentioning. Firstly, the choice of POD as a decomposition technique was here motivated by its simplicity and the possibility of interpretation when coupled to a linear learning formulation such as PCE. Other decomposition techniques exist, and many authors attempted comparisons, for example with Fourier [84], extensions of POD [30,37] or other classes of decomposition [90]. For the present application, other decomposition techniques such as Kernel Principal Component Analysis (KPCA) [69] and Sparse PCA [43] were analyzed, without significant improvements. Secondly, data filtering using POD consisted only in deleting the poor-quality measurements and extracting the spatial zones where data were always measured. POD can however be used to reconstruct missing data, by inverse projection on POD basis elements deduced from qualitative data [89]. This could help to extend the statistical set for the learning. Lastly, a linear interpolation of the bathymetry was used to project all the measurements onto the same grid for POD application. The uncertainty that emerged from this interpolation process was not treated. This uncertainty, added to the measurement errors, can shed light on model behavior. For example, comparison of mono-beam cross-sections with multi-beam measurements and uncertainty propagation of bathymetry errors through the learning could be attempted, especially because uncertainties in the bathymetric information may impact the flow field computation [11,60].
PCE was used to learn the temporal POD modes as 1D data. We showed the importance of the polynomial basis and therefore of the choice of marginals for the learning phase. In this comparison, we showed that choosing uniform distributions associated with Legendre family might not be appropriate, even though it is widely used when no input information is available [103]. Moreover, the number of inputs involved can alter the learning. When using LARS, the presence of numerous variables can mislead the algorithm in selecting useless variables that seem to decrease training error but increase prediction error. This is an overfitting phenomenon. A good combination between polynomial basis and dimension choice could significantly improve convergence speed, centering of residuals and mean training and prediction errors. The proposed contribution analysis using the PCE coefficients showed that the last-state information is often the most influential input. A robustness test was conducted by varying the training set, and the observation was stable. For the modes of small ranks associated with the largest variances (modes 1 to 4 in our example), wave height was the most influential forcing, whatever the chosen learning set. This is consistent with physical knowledge of sediment mobilization in coastal configurations, where waves are known to be determining through the influence they have on bed shear stress [105]. For modes of higher rank, however, the only selected variable by LARS is the last-state information. Firstly, the forcings that were used as PCE inputs were simple statistical estimators deduced from the data (means, percentiles, etc.). This reduction was used instead of giving all the time series as an input, because the problem would become ultrahigh-dimensional. This wastes the richness of the available data for inputs such as tidal information that are measured on an hourly basis. A more accurate statistical reduction of the inputs could be used. For example, Lataniotis et al. [54] used PCA and KPCA for surrogate modeling with PCE and Gaussian Processes on ultrahigh-dimensional problems. Secondly, dependencies were not specifically modeled. These can be incorporated using the mathematical setting for the construction of the polynomial basis established by Soize and Ghanem [97] . The dependencies, however, indirectly influenced the construction of the model via selection of basis elements by LARS, which avoids redundancy. Thirdly, the choice of tested input configurations for PCE was arbitrary. A more objective variable selection technique is necessary [22,81]. Lastly, PCE was chosen for the interpretation possibilities that it allows when combined to POD, thanks to computation of importance measures from the expansion coefficients. Other methods can be used: i.e., for high rank POD modes that were poorly learned using PCE. Indeed, PCE seems to work better for modes associated with high than low variances. Although it may be tempting to conclude that modeling of high rank modes is not necessary since they are associated with low variances, it should be noted that accurate prediction of modes of higher rank can make the difference between average forecasting and forecasting that captures less frequent events and smaller scale features of the 2D field. Therefore, the present ML could be enhanced by improving the learning of high-rank modes. For example, the construction of marginals and the use of random draw with confidence intervals, or extreme statistics models [26], instead of causal models like PCE, could be attempted.
Finally, the robustness and convergence properties added to the physical interpretability supported the choice of POD-PCE coupling as a ML prediction algorithm. It respects the PDR (Predictive, Descriptive, Relevant) framework defined in [75]. It is characterized by both predictive and descriptive accuracy (simplicity) and is stable with respect to data disturbance. It is interpretable, as the sparsity, simulatability and modularity defined in [75] are respected by construction. Finally, it is both interpretable at features level (POD components and their PCE) and at multidimensional output level (GW compared to the proposed GGW indicators). The POD-PCE ML was therefore implemented using the first 11 modes, after sensitivity test to number of modes. Mean information (e.g. sediment deposition rate) was in general well reproduced. Profile-by-profile investigation also showed that POD-PCE coupling was promising, as the spatial distribution of the sediment deposition patch locations and amplitudes were well represented. Some general limitations should be highlighted and could be good perspectives for improving the process. The small data-set was a clear handicap in our problem. Some events, such as sediment downstream the intake or variation in measurement intervals, were poorly represented. It would be interesting to test the methodology on an enriched data set in order to assess the real potential of POD-PCE Machine Learning. Due to the lack of such data, input distributions were certainly not well approximated. One way of improving POD-PCE coupling would be the development of hybrid measurement-based/process-based data learning [72,93]. This could be used to enrich the data set, not only by increasing its size (emulated realistic scenarios) but also by adding new input parameters that are not measured but obtained from process-based modeling. Nearshore grid point - Table 1: Summary of measurements, frequencies, periods, sources and spatial coverage. Descriptive statistics are given (P: Percentile).
The forcings did not have the same frequencies. A solution for the homogenization of frequencies is to reduce the measured data to representative statistics over the sedimentation interval ∆t ≈ 15 days that separates two bed elevation measurements. Hence, the statistics described in Table 2   Scatter plots of some input variables are presented in Figure 23.

Explicit fitted PCE model
For the first three modes, fitting with the model H P i ; Stlj is illustrated in Figure 26, and shows good pointwise performance on the fitting and learning intervals. The corresponding explicit PCE models for each POD mode are shown in Equation 22.