1 Introduction
Many problems in computational mechanics consist of known physics described by partial differential equations (PDEs) for some fields of interest (e.g. fluid velocity, solid displacement) but include other physical fields whose truth values are not known exactly. For instance these uncertain fields can be (1) initial conditions of the fields of interest, (2) constant physical properties such as mass (density field) or stress distribution, or (3) physical fields related to the fields of interest through unknown relationships, i.e. the problem is unclosed. Here we refer to such fields as
latent physical fields, as they are embedded in the PDEs but not solved for directly. In inverse problems these are treated as model parameters that need to be inferred based on observations of the field of interest. For fields of interest and latent fields such problems can be written as(1) 
where represents the governing partial differential equations (PDEs) including geometry and boundary conditions on . The fields of interest can be considered a function of the latent fields as
(2) 
An example of such a problem is the Reynoldsaveraged Navier–Stokes (RANS) equations, which describe the mean velocity and pressure fields of a fluid but include the unclosed Reynolds stress field. The unclosed Reynolds stress field is approximated through imprecise turbulence models. Similarly, in many solid mechanics problems the field of interest is displacement and the latent field is the stress in the material, which can be a highly nonlinear function of the strain and strain history. In the rest of this work the RANS equations are used as the working example. We are interested in the inverse problem: inferring the true latent fields using an initial estimate of the latent fields and sparse experimental observations of the output fields.
1.1 Inverse Problems in Computational Mechanics
There are many frameworks to improve the predicted fields by using sparse observations (e.g., kennedy_bayesian_2001, ; oliver_validating_2015, ) (duraisamy_turbulence_2019, ; xiao_quantification_2018, , for comprehensive reviews for RANS simulations), and the choice is between different representations of the uncertainty in the problem. If the governing PDEs are considered exact and therefore the only uncertainty is in the latent field, then the embedded uncertainty formulation by Oliver et al. oliver_validating_2015 is most appropriate. In this formulation the problem is written as as in Equation (2) and the latent fields are considered uncertain and need to be inferred. This is in contrast to black box approaches where a correction to the model output (e.g. in Equation (2) for some approximation of ) is inferred. The embedded discrepancy formulation is used in this work since the Reynolds stress field is considered the major source of uncertainty in the RANS equations.
In Bayesian inversion approaches the uncertain field is considered a random field, and a statistical model describing this field must be specified. This is referred to as the prior distribution, which is later updated based on the observation data to obtain a posterior distribution. The chosen statistical model, particularly its covariance kernel, can have a large impact in the Bayesian inference. For simplicity in writing the equations, in the rest of the paper it will be assumed that there is a single scalar latent field which needs to be modeled as a random field. The initial approximation (baseline solution) of the latent field is usually close to the truth, and, in the absence of additional information, the latent field is taken to be a Gaussian process with mean equal to this baseline solution (e.g. xiao_quantifying_2016, ; dow_quantification_2011, ; singh_using_2016, ; cheung_bayesian_2011, ). The prior distribution is then
(3) 
where is the specified covariance. Using the embedded discrepancy formulation described by Equations (2) and (3), observations of the output fields can be used to infer a posterior distribution for the latent field. The inferred latent field can then be propagated through the model in Equation (2) to obtain the inferred output fields of interest. There are many approaches to solve the Bayesian inversion problem directly or approximately, including ensemble based methods (e.g. xiao_quantifying_2016, ; iglesias_ensemble_2013, ) and optimization/adjoint based methods (e.g. dow_quantification_2011, ; singh_using_2016, ).
The Bayesian embedded discrepancy framework has been used in several works to address the RANS equations closure problem. Edeling et al. edeling2014bayesian ; edeling2014predictive and Ray et al. ray2016bayesian ; ray_robust_2018 ; ray2018learning
inferred model coefficients in turbulence models by using Markov Chain Monte Carlo methods. Dow and Wang
dow_quantification_2011inferred the eddy viscosity by observing the velocity field, which was among the first efforts to infer a full field rather than model coefficients. Their work used the linear eddy viscosity approximation in the RANS equations, which allows the Reynolds stress tensor field to be represented by a single scalar eddy viscosity field. They obtained the maximum likelihood estimate using an adjoint method to compute the derivatives. Similarly, Singh and Duraisamy
singh_using_2016 also inferred the eddy viscosity by observing velocities but do so by inferring an embedded multiplicative discrepancy into the production term of the Spalart–Allmaras model for eddy viscosity. They obtained the maximum a posteriori (MAP) estimate using an adjoint method to compute the derivatives. Xiao et al. xiao_quantifying_2016 inferred the components of the Reynolds stress tensor field by observing velocities at sparse locations. They used an iterative ensemble Kalman method to infer the Reynolds stress field. Meldi et al. meldi2017reduced ; meldi2018augmented integrated the Kalman method into CFD solvers in an intrusive manner for data assimilation of turbulent flows. Their approach has the advantage of involving only a single simulation rather than an ensemble of simulations as in the ensemble Kalman methods. More recently, Edeling et al. edeling2018data proposed transport equations to describe the deviation from the eddy viscosity model and use Bayesian inference to calibrate the constants in the transport equations. In related efforts, researchers have also used datadriven methods to address closure problems in large eddy simulations vollant2017subgrid and multiphase flows ma2015using ; ma2016using .1.2 IllPosedness and Enforcing Physical Constraints
Inverse problems are inherently illposed in that numerous possible latent fields can result in output fields that are close to the observation values. One example is estimating the Reynolds stress from inplane velocity measurements in a square duct case xiao_quantifying_2016 . In this flow the inplane velocities are driven by the imbalance in the Reynolds stress, and any Reynolds stress field with the correct imbalance is a possible solution. Illposed problems are usually regularized by enforcing additional constraints, such as smoothness of the inferred field. For example, smoothness can be enforced by adding the magnitude of the gradient into the cost function.
Because of the illposed nature of the problem, most previous research fails at inferring the true latent field while still obtaining global improvements in the fields of interest. Since the latent field is a physical field with known physical constraints, the search space for possible latent fields can be reduced by enforcing these constraints before any additional regularization. One type of physical constraints are global physical constraints that the solution must conform to in the entire domain. These are constraints such as the divergencefree requirement for an incompressible velocity field or the positivity requirement for a pollutant concentration field. Some of these global constraints can be enforced by using a different choice for the statistical model of the random fields in Equation (3), i.e. using a more complex representation than a simple Gaussian process. The representation of the random latent field is chosen such that the physical constraints are automatically met by any realization of the random field. For example, a positivity constraint can be enforced by modeling the latent field as a lognormal process as
(4a)  
(4b) 
where the baseline solution is now the median of the distribution, and there is an implicit normalization of , , and . Xiao et al. xiao_quantifying_2016 use this statistical model to enforce positivity on the turbulent kinetic energy while Dow and Wang dow_quantification_2011 use it to enforce positivity on the eddy viscosity.
A second type of physical constraints is boundary conditions (B.C.) on the latent field. Boundary conditions are physical constraints on fields at the boundaries of the domain. Typical examples include fixed values (Dirichlet B.C.), fixed gradient (Neumann B.C.), fixed value and fixed gradient (mixed B.C.), and periodic conditions. For example, in fluid mechanics applications, (i) solid walls enforce a zerovelocity condition referred to as the noslip condition, (ii) symmetry planes (used as a computational simplification for symmetric problems) enforce a zerogradient condition on all fields, and (iii) far field conditions enforce zero gradient and fixed value conditions on all fields. Many works have simply ignored some of these boundary conditions and still obtained improvements in the output fields thanks to the illposedness of the problem. For instance Xiao et al. xiao_quantifying_2016 do not enforce the periodic boundary condition on the Reynolds stress field. Similarly Dow and Wang dow_quantification_2011 do not enforce the symmetry condition on the latent eddy viscosity field. On the other hand, simpler boundary conditions have been often enforced in a problemspecific manner. Particularly, the two works described above xiao_quantifying_2016 ; dow_quantification_2011 both enforce the zerovalue boundary condition at the wall on the latent fields thanks to the multiplicative nature of lognormal process used as the statistical model. Similarly, Cheung et al. cheung_bayesian_2011 enforce the noslip (zerovalue Dirichlet) condition on the inferred velocity by inferring a multiplicative discrepancy on a baseline velocity that already satisfies the noslip condition. However, there clearly has been no uniform approach to enforce general boundary conditions on physical latent fields in Bayesian inversion problems. In this work we demonstrate how to enforce arbitrary boundary conditions on the latent fields by choice of statistical model.
1.3 Contribution of Current Work
In this paper we demonstrate how to enforce general boundary conditions on physical latent fields by choice of the statistical model used to represent them in the Bayesian inversion. We start with the assumption that the latent field is represented by a Gaussian process or a function thereof (e.g. a lognormal process), and the boundary constraints are then enforced by choice of covariance matrix. An appropriate covariance matrix is obtained by first defining an initial Gaussian process and then modifying the covariance matrix in a way that enforces the boundary conditions. We also present how to enforce internal observations of the latent field, which is done in a mathematically similar way by modifying not only the covariance but also the mean of the initial Gaussian process. This can be a realistic scenario since some observations of the latent field can be obtained in practice. For example, Reynolds stress can be obtained from instantaneous measurements of the velocity components. Enforcing both of these types of constraints further reduces the search space in the optimization for the latent field by ensuring it satisfies more of the known physics, driving the inference closer to the goal of obtaining the true latent field. The method is tested for fluid mechanics problems using the RANS equations, but we emphasize that the method is general with applications to many PDEdescribed physics problems.
1.4 Related Work
In what follows we differentiate our contribution from other works using Gaussian processes to represent physical fields. Gaussian processes are commonly used as the Bayesian framework for field inferrence from observations gibbs1998bayesian ; rasmussen1999evaluation (e.g. shahriari2015taking, ; hensman2013gaussian, , forcontemporaryapplications). In these cases the field being inferred is the same field being observed. The class of problems considered here are inverse problems where the inferred field is related to the observed field through some model (set of PDEs, e.g. RANS equations). We focus on the statistical model used to represent the random field, and not on the particular choice of Bayesian framework used in the inversion. Specifically we seek a choice of statistical model that enforces the boundary conditions for all realizations of the random field. However, the methodology showcased here is directly adopted from these works that use the Gaussian process as the Bayesian framework. The theory of incorporating observed values and derivatives of the random field into the Gaussian process is well developed solak_derivative_2003 ; rasmussen_gaussian_2006 , but, to the author’s knowledge, it has not been used in the context of enforcing boundary conditions on the inferred fields in inverse modeling. In many such problems the boundary conditions have simply been ignored as illustrated earlier.
Similarly, there are many works using Gaussian processes to represent solutions to stochastic PDEs (e.g. sarkka2011linear, ; graepel2003solving, ; raissi2017inferring, ). In these works the form of the PDE is embedded into the covariance kernel and realizations of this random fields satisfy the PDE. These are powerful techniques for solving PDEs in a stochastic sense and can be naturally used to solve the inverse problem for model parameters in the PDE (e.g. raissi2017machine, ; dondelinger2013ode, )). However this is not directly related to the class of problems considered here. Here we consider PDEs with an uncertain latent field which needs to be inferred, but the PDEs are solved deterministically for each proposed latent field during the chosen Bayesian optimization process. As mentioned earlier our focus is seeking a statistical model for the latent field that incorporates known physical constraints. Recently Wu et al. wu2019physics presented a method that incorporates some aspects of Gaussian processes for stochastic PDEs into the problem of specifying a statistical model that incorporates known physical constraints. Specifically, for problems where the latent field is itself governed by PDEs with unclosed fields the structure of these PDEs are incorporated into the covariance kernel of the latent field. Note that this new set of PDEs could be solved directly in conjuction with the main set of PDEs and the latent fields would be the uncertain terms in these additional PDEs. In this sense the method proposed in wu2019physics represents a hybrid approach with the benefit of reducing the number of coupled PDEs that need to be solved for each forward evaluation and reducing the number of latent fields that need to be inferred.
While the current work applies to any Bayesian inversion technique, an ensemblebased algorithm is used in the test cases. There are numerous ongoing efforts in improving existing ensemblebased Bayesian inversion algorithms to allow for incorporating physical constraints wu2019adding ; wu_improving_2019 ; zhang2019regularization . The current work and that of Wu et al. wu2019physics complement these efforts. The common overarching goal of these works is inferring latent physical fields with limited observation data.
The rest of the paper is organized as follows. Section 2 presents the method of embedding known boundary behavior into the prior statistical model. Section 3 presents the results from a series of test cases, where it is shown that all prior and posterior realizations of the random field satisfy the boundary conditions. Finally, Section 4 concludes the paper and offers some discussions.
2 Methodology
The goal is to ensure that the realizations of the random latent field satisfy any known values and derivatives, be that at internal locations or boundary conditions, by appropriate choice of statistical model for the random field. Internal observations come from experimental measurements and are fundamentally different from boundary conditions which come from physical constraints and modeling choices. However, they can be enforced in a mathematically similar way and thus both are considered here. We consider the cases where the latent field is treated as a Gaussian process (Equation (3)) or a function of a Gaussian process. Incorporating known values and derivatives into a Gaussian process consists of modifying the mean and covariance of the Gaussian process solak_derivative_2003 ; rasmussen_gaussian_2006 . That is, after observations Equation (3) becomes
(5) 
where the modified mean conforms to all known values and derivatives, and the modified covariance ensures all realizations of the random field do too, within observation errors. This theory is summarized in Section 2.1 and is used to enforce internal observations. In Section 2.2 the application to enforcing boundary conditions is discussed.
In the case where the latent field is modeled as a function of a Gaussian process as
(6a)  
(6b) 
the desired constraints on the latent field first need to be mapped to constraints on the Gaussian process . These modified constraints are then enforced in the same manner as for the general Gaussian process. In Section 2.3 we describe this procedure for the general case (Equation (6)) and show the results for the lognormal process (Equation (4)) in particular. The lognormal process not only serves as a specific example but is a particularly useful and prevalent formulation in fluid mechanics problems, including the problem we use in the results section to showcase this methodology.
To verify the approach, a KarhunenLoève (KL) representation of the latent field is used and it is shown that the KL modes satisfy the boundary conditions. A stochastic process can be represented as an infinite linear combination of orthogonal functions through the KarhunenLoève theorem. The latent field in Equation (5) can then be written as
(7a)  
(7b) 
where and
are the eigenvalues and unit eigenfunctions (modes) from the KL decomposition and
is the standard normal distribution. Realizations of the random process can be obtained by sampling the coefficients
, which have independent standard normal distributions. When the field is discretized, the limit of the summation becomes the number of cells. However, the modes corresponding to the largest eigenvalues account for most of the covariance and it is only necessary to retain a small subset of modes corresponding to the largest eigenvalues. The number of retained modescan be chosen based on the desired coverage (e.g. 99%) of the variance of the field. From Equation (
7) it can be seen that for the realizations to exactly satisfy known values at some locations the KL modes must be zero at those locations. To allow for observation errors, the KL modes must have small magnitudes at these locations. Similarly, the derivatives of the KL modes must be near zero at locations of known derivatives. Therefore, it is sufficient to simply verify that the KL modes satisfy these conditions to guarantee all realizations of the latent field will satisfy the desired conditions. Note that in the case of boundary conditions the values or derivatives of KL modes at these location must be exactly zero, since there are no observation errors.In addition to showing that the modified mean and covariance guarantee any realization of the latent field satisfies the boundary conditions and internal observations, we also solve the inverse problem. The inverse problem consists of inferring the latent field from sparse observations of the output fields. There are many frameworks for solving the Bayesian inversion problem, from minimization of an objective function using adjoint methods to ensemblebased methods. Here we adopt the framework from Xiao et al. xiao_quantifying_2016 , which is based on the iterative ensemble Kalman method developed by Iglesias et al. iglesias_ensemble_2013 . The details of the Bayesian inversion are presented in A
. The mean of the posterior distribution is used as an estimate for inferred latent field, and correspond to the maximum a posteriori estimate.
For the Bayesian inversion we used a reduced order model consisting of the KL representation in Equation (7) with a truncated set of modes to capture
of the variance. The state vector to be inferred now consists of the truncated coefficients
rather than the discretized values of the latent field directly. It is clear that with this representation the inferred latent field will satisfy the boundary conditions since the reduced order model is also a weighted sum of KL modes and the mean (Equation (7)). This representation is, however, not necessary since the ensemble Kalman filter ensures that the posterior samples are linear combinations of the prior samples
iglesias_ensemble_2013 . Note that if the observation error is not zero, the inferred state vector could result in a latent field that does not satisfies the known values or derivatives. This could be remedied with regularization to penalize the inferred state vector based on the likelihood that it came from the prior distribution, i.e. where each coefficient has a standard normal distribution. For simplicity, in the test cases presented, we solve the inverse problem enforcing boundary conditions alone, without any internal observations of the latent field.2.1 General Constraints on Gaussian Processes
If the values of a Gaussian process or its derivatives are known at some points, the Gaussian process can be modified to reflect this prior knowledge solak_derivative_2003 ; rasmussen_gaussian_2006 . To incorporate the known values and derivatives, first an augmented vector is defined as where is the subvector of the values of the Gaussian process (Equation (3)) at the inference locations and is the subvector of known values and derivatives of the process at some points. From this point on all fields are discretized unless otherwise noted but the same symbols are used for simplicity. The locations for consist of the cell centers of the discretized domain, and the locations for are not necessarily colocated with points in . Known derivatives in are defined with respect to one of three orthogonal coordinate directions , , or . To account for derivatives in general directions, the derivative in the three coordinate directions are included in separately for the same physical location. A Gaussian process is then assumed for the vector as
(8) 
where
are the values of the baseline solution or its derivatives at the observation points, possibly obtained by interpolation, and
, , , and are submatrices of the covariance matrix. Finally, the distribution for inference values conditioned on known values and derivative is then also a Gaussian process given by solak_derivative_2003 ; rasmussen_gaussian_2006(9a)  
(9b)  
(9c) 
where is the vector of known values and derivatives and the corresponding error (covariance) matrix. Equation (9) now corresponds to Equation (5) but enforces the known values and derivatives, within specified error, on any realization of the latent field.
The covariance matrix in Equation (8) is created from a chosen continuous covariance kernel for the covariance between the values of at different physical locations and . The covariance between two values, between a value and a derivative, and between two derivatives are given by solak_derivative_2003 ; rasmussen_gaussian_2006
(10a)  
(10b)  
(10c) 
respectively. The most common covariance kernel is the square exponential kernel , which depends only on the physical distance between two points. The square exponential covariance between two points is given by Equation (11a), where is the variance and for are the length scales in each of the three spatial coordinate directions. With this choice of kernel, Equation (10) becomes
(11a)  
(11b)  
(11c) 
where is the Kronecker delta.
2.2 Boundary Conditions on Gaussian Processes
When enforcing Dirichlet, Neumann, and mixed boundary conditions there are several simplifications to the process described in Section 2.1. In Equation (8) now represent the known boundary values and derivatives at the centers of the boundary faces in the discretization. Since boundary conditions come from the modeled physics and are also enforced exactly on the output fields in the forward model, no observation errors are used, i.e. . The baseline solution is assumed to satisfy the boundary conditions, and the difference between the observed value and the baseline value at the boundaries is then zero, i.e. . With these simplifications Equation (9) becomes
(12) 
Equation (12) shows that only the covariance needs to be updated for enforcing boundary conditions and the mean is unchanged.
Periodic boundary conditions are also very common in computational physics, where they are used in modeling problems with certain types of symmetry. Periodic boundary conditions cannot be implemented through the modifications presented in Section 2.1, instead they can be implemented by choice of covariance kernel . For instance, if the problem is periodic in the direction , a mixed kernel can be used with periodic covariance in direction and square exponential in the other two directions as
(13) 
where is the periodicity. In order for the length scale in the periodic direction to have a similar effect as a length scale in the square exponential kernel, it is chosen as
2.3 Constraints on Functions of Gaussian Processes
If the latent field is described by a function of a Gaussian process (Equation (6)), the constraints on the latent field need to be mapped to constraints on the Gaussian process. The function is assumed to be invertible and differentiable. If the value of the latent field at some point is known to be , the corresponding value for the Gaussian process at that point is
(14) 
Similarly, if a derivative in direction of the latent field at some point is known to be , the corresponding constraints on the Gaussian process is
(15) 
where is defined as
However, Equation (15) does not have a unique solution since there are two variables, and its derivative , and only one relation. Enforcing observations of the derivative of the latent field requires an additional constraint. Specifically, if the value at location is also specified, Equation (14) can be substituted into Equation (15). It is therefore possible to enforce observations of the value of the latent field at a point (e.g. Dirichlet B.C.) and to enforce combined observations of the value and derivative at a point (e.g. mixed B.C.) but not to enforce a derivative observation alone (e.g. Neumann B.C.). The last type of boundary conditions considered, periodic boundary conditions, is enforced on the latent field with periodicity by requiring the covariance of to be periodic with the same period.
The lognormal process () is useful in problems where the physical latent field must be nonnegative. If the process is lognormal then it has the form in Equation (4), and Equations (14) and (15) become
(16) 
and
(17) 
respectively. In general the derivative constraint still has infinite solutions and requires both the value and slope of the process to be defined at that point. There is however a unique solution for the special case of zerogradient constraint, which is by far the most common derivative constraint in fluid mechanics. In this case and Equation (17) becomes
(18) 
One issue with the lognormal process is that it is singular when the field has a zerovalue condition, i.e. in Equation (16), and the common case of a zerovalued Dirichlet boundary is explored here in more detail. The problem could be addressed numerically by having the baseline solution have small value at the boundary. However, this would still cause the value of the Gaussian process in Equation (16) to be very large, albeit finite. This is also a problem for incorporating internal observations near such boundaries, where the baseline solution is close to zero. In such a case, a small difference between known values and the baseline solution corresponds to a large magnitude of the Gaussian process, approaching infinity as the baseline solution approaches zero. This is most obvious when the normalization in Equation (4) is made explicit and the baseline is used for normalization. This gives
(19a)  
(19b) 
which emphasizes the multiplicative nature of the lognormal formulation. At a location where the baseline is close to zero, even a slightly different value of the latent field would require a large value of the Gaussian process . To address this singularity at the boundary we use the formulation in Equation (19) and take the baseline solution to satisfy the boundary condition. In this case the value of the observation and Equation (16) becomes
(20) 
for boundaries with zerovalue Dirichlet conditions. Internal observations in the vicinity such boundaries should still be avoided. Another consequence of this singularity and the representation in Equation (19) is that since the baseline solution is already close to zero near the boundary, the boundary condition is enforced without need to modify the covariance. Here, however we do enforce the appropriate condition, Equation (20), at such boundaries. In summary, when using a lognormal distribution as in Equation (19) a zerovalued Dirichlet condition on the latent field implies the following:

no internal observations of the latent field should be enforced near the boundary,

the equivalent boundary condition for the Gaussian process is as in Equation (20), however

the boundary condition is defacto enforced with no need to modify the covariance matrix.
3 Results
Two fluid flow problems are used as test cases: flow through a channel and flow over periodic hills. The latent field in both cases is the eddy viscosity, which has a positivity requirement and therefore a lognormal process is used to represent it. In Section 3.2 we perform an a priori study to highlight the methodology. The boundary conditions are enforced for both cases and an internal observation is enforced in the periodic hills case. In Section 3.3 Bayesian inversion is performed on both cases using sparse observations.
3.1 Test Cases
Two test cases are used to showcase the methodology. In both cases we solve the steady, incompressible, Reynoldsaveraged NavierStokes (RANS) equations describing the mean properties of fluid flows. The linear eddy viscosity approximation is used to represent the Reynolds stress tensor field with a single scalar field, the eddy viscosity field. The RANS equations are then
(21a)  
(21b) 
and correspond to the model in Equations (1) and (2), where the output fields () are the components of velocity () and a pseudo pressure term () and the latent field is the eddy viscosity (). The indices and refer to spatial direction, e.g. for a twodimensional problem. The dynamic viscosity is a known physical property of the fluid.
Both cases consist of wallbounded incompressible flow driven by a pressure gradient source to achieve the desired bulk velocity. The first case is a steady, fully developed, onedimensional channel flow with top and bottom walls, shown in Figure 0(a). The case has a Reynolds number of based on half channel height and bulk velocity . To model this problem a single cell is used in the direction of flow with periodic boundary conditions on the boundaries perpendicular to the flow. Because of symmetry only the bottom half of the channel is modeled, using a symmetry boundary condition on the top boundary. The domain is discretized using cells in the direction normal to flow, as shown in Figure 0(b).
. The discretized onedimensional simulation domain (b) consists of only one cell and periodic boundary conditions in the streamwise direction since the interest is in the developed velocity profile. Because of symmetry only the bottom half of the domain is modeled. The baseline (initial guess) and true solutions for the latent eddy viscosity field are shown in (c) and the problem consists of inferring the multiplicative discrepancy (d) between these. The discrepancy in (d) is normalized based on the standard deviation of the prior distribution at each cell.
The second case is the twodimensional infinite hills geometry with top and bottom walls created by Mellen et al. mellen_large_2000 and shown in Figure 1(a). The case has a Reynolds number of based on the hill height and bulk velocity through the section at the hill top . This is modeled as a single hill with periodic boundary conditions and the domain was discretized into cells, as shown in Figure 1(b). To highlight the periodicity, the upstream half of the computational domain is plotted starting after the downstream half in the rest of the paper. Furthermore, a sequential colormap is used for plots of quantities that are nonnegative (e.g. eddy viscosity), and a diverging colormap is used for plots of signed quantities (e.g. dicrepancy, KL modes).
The periodic boundary conditions require all fields to have the same values and derivatives on both boundaries while the symmetric boundary condition requires the derivative in the direction normal to the boundary to be zero. These requirements are for all fields including the eddy viscosity field. The wall boundary condition translates to different requirements for different fields. The velocities and are zero at the wall while the pressure has zero gradient normal to the wall. The wall imposes a zerovalue Dirichlet condition on the eddy viscosity field since the Reynolds stresses are zero at the wall and the velocity gradient is not necessarily zero. The eddy viscosity has a positivity constraint, and this is enforced by modeling it as a lognormal process as in Equation (4). The boundary conditions are therefore not implemented directly as described in Section 2.2 but with the modifications described in Section 2.3.
For both the channel and periodic hills cases two solutions are required, the baseline solution and a synthetic truth. The baseline solution is used as the median for the prior statistical model describing the latent field. The synthetic truth is used to create synthetic observations based on perturbations of this truth and to evaluate the performance of the Bayesian inference. The RANS equations with linear eddy viscosity approximation are also used in creating the synthetic truth. This is because the RANS equations are considered to contribute little to the uncertainty of the solution, justifying the embedded discrepancy formulation. Solving the RANS equations requires choosing a turbulence model for the eddy viscosity. Different turbulence models were chosen to create baseline and synthetic truth solutions that were sufficiently different to illustrate the methodology, and these are summarized in Table 1. For the channel case the Spalart–Allmaras spalart1992oneequation and – jones1972prediction turbulence models are used for creating the baseline and synthetic truth solutions, respectively. For the periodic hills case the – turbulence model wilcox2006turbulence is used for both the baseline and synthetic truth solutions, but the coefficient is modified to for the latter instead of its standard value of . The synthetic observations are obtained by taking a realization of a multivariate normal distribution with mean equal to the synthetic truth and with assumed uncorrelated errors.
The baseline solution, synthetic truth, and logarithm of the multiplicative discrepancy are shown in Figure 1 for the channel case and Figure 2 for the periodic hills. The baseline solutions are used for the formulation in Equation (19) and the logarithm of the multiplicative discrepancy is the value of the Gaussian process that needs to be inferred. Note that the baseline solutions satisfy all boundary conditions. As discussed in Section 2.3, the lognormal formulation has singularities whenever the baseline solution is zero, as is the case on walls on both test cases. This can be seen in the large values for the discrepancies for both cases. Because of this we will avoid making observations near the singularities, which would drive the entire inferred field to an unrealistic solution. The inferred solutions are presented in Section 3.3, as well as quantitative measures of the discrepancies. To quantify the discrepancy between a field and the true field we use
(22) 
where denotes the norm of the field.
Baseline  Truth  

Channel  SpalartAllmaras   
Hills  standard    with 
3.2 Enforcing Boundary Conditions and Internal Observations
This section presents the a priori results of enforcing boundary conditions and internal observations. Section 3.2.1 presents the results of enforcing boundary conditions on the channel case. Section 3.2.2 presents the results of enforcing boundary conditions on the periodic hills case, a more complex, twodimensional case. Finally, Section 3.2.3 presents the results of enforcing an internal observation, in addition to the boundary conditions, for the periodic hills case. In all cases it is shown that all realizations of the modified fields satisfy the boundary conditions and the internal observations.
3.2.1 Boundary Conditions  Channel
For the channel case the eddy viscosity has (1) a zerovalue Dirichlet condition at the wall, (2) a zerogradient Neumann condition at the symmetry plane, and (3) a periodic boundary condition in the flow direction. The boundary conditions are enforced by modifying the covariance of the Gaussian process in Equation (19). In general, the Dirichlet boundary condition at the wall requires the value of the Gaussian process to be zero at the boundary. As discussed in Section 2.3 it is not necessary to enforce this condition but we still chose to do so. The zerogradient Neumann condition at the symmetry plane requires that the gradient of normal to the symmetry plane be zero. The periodic boundary condition is implicitly enforced since there is a single cell in the direction of periodicity.
The original covariance matrix consists of a square exponential kernel with variance and length scale . This is then modified to enforce the Dirichlet and Neumann boundary conditions, by using Equation (12). The original (square exponential) and modified covariance matrices are shown in Figures 2(a) and 2(b), respectively. Close to the bottom wall (), where a Dirichlet boundary condition is enforced, the modified covariance is near zero. This results in the value of the baseline near the wall being used for all realizations. Note that the boundary locations (faces) are not in the covariance matrices in Figures 2(a)2(b), since the inferred quantities are the internal cell values. Near the symmetry plane (), where a Neumann boundary condition is enforced, there is strong covariance with the immediate neighbors and a quick decrease as you move further away. Since the baseline solution has zero gradient at this boundary, this ensures that the gradient of any realization is also zero at the boundary.
As discussed earlier, ensuring that the KL modes satisfy their corresponding conditions is enough to ensure that any realization of the latent field satisfies the boundary conditions. The original and modified covariances require and modes respectively to capture of the variance. These number of modes are used for the truncation. Figures 2(c) and 2(d) show the first five modes using the original and modified covariance matrices, respectively. It can be seen that the modified covariance does enforces zero gradients on the modes whereas the original covariance does not, as expected. While it was shown that this condition is sufficient to ensure the realizations of eddy viscosity satisfy the zero gradient boundary conditions, it is helpful to actually visualize this. The realizations are created as in Equation (7) using the truncated modes. Figures 2(e) and 2(f) show created realizations using the original and modified covariance matrices, respectively. It can be seen that the gradient at the symmetry plane is zero for all realizations when using the modified covariance matrix, but is not the case when using the original covariance matrix.
3.2.2 Boundary Conditions  Periodic Hills
For the periodic hills case the eddy viscosity has a periodic boundary condition in the direction of the flow and zerovalue Dirichlet conditions at the top and bottom walls. The periodic boundary condition is enforced by making the covariance kernel periodic in the direction, as in Equation (13), with periodicity of . Making periodic ensures that the eddy viscosity is also periodic with the same periodicity. The Dirichlet condition is enforced on the bottom and top walls using Equation (12). Again, it is not necessary to enforce this condition but we still chose to do so.
The original and modified covariance matrices are constructed as before, using and . Since this case is twodimensional the covariance matrix is difficult to interpret and depends on the cell ordering. Instead, Figures 3(a) and 3(b) show the covariance at a selected point for the original and modified covariance matrices, respectively. The modified covariance can be seen to be periodic and the covariance with the wall is zero, as expected. Note that the plotted results are shown with the first half of the simulation domain shifted behind the second half to highlight when the periodicity constraint is met.
Ensuring that the KL modes satisfy their corresponding conditions is enough to ensure that any realization of the latent field satisfies the boundary conditions. Figures 3(c) and 3(d) show a selected KL mode using the original and modified covariance matrices. As can be seen the mode of the modified covariance matrix satisfies the periodicity and Dirichlet conditions, while that for the original covariance matrix do not. All other modes using the modified covariance matrix similarly satisfy the boundary conditions. For the truncation the number of modes where chosen as and for original and modified covariances, respectively, to cover of the covariance. Figures 3(e) and 2(f) show a realization of the eddy viscosity field using the original and modified covariance matrices. It can be seen that the modification to the covariance matrix results in the boundary conditions being enforced on the realizations of the eddy viscosity field. Specifically, the periodic and zerovalue boundary conditions are satisfied for all realizations of the eddy viscosity field. Note that the zerovalue boundary conditions at the bottom and top walls are also satisfied when using the original covariance as disussed in Section 2.3.

3.2.3 Internal Observations  Periodic Hills
As a last test case the observed value of the eddy viscosity at an internal point is enforced in the periodic hills case. The observation is made at and the true value of the eddy viscosity is used. The observation is considered exact with no observation error. The mean and covariance are modified using Equation (9). The modified baseline eddy viscosity is shown in Figure 4(a) for the crosssection passing through the observation point. The modified baseline now satisfies the observation value. The first five KL modes for the modified covariance are shown in Figure 4(b) for the same crosssection, and can all be seen to have a value of zero at the observation location. Figure 4(c) shows ten realizations of the eddy viscosity field at the same cross section. All realizations satisfy the internal observation, as desired.
3.3 Bayesian Inversion
The inverse problem was solved for both the channel case and the periodic hills case using the Bayesian ensemble Kalman scheme presented in A. The boundary conditions where enforced during the inversion by using the modified covariance matrices in Sections 3.2.1 and 3.2.2 for the channel and periodic hills cases respectively. The inverse problem consists of inferring the eddy viscosity based on the prior distribution and sparse observations of the velocity field. As discussed in Section 2 a reduced order model is used and the state vector to be inferred consists of the coefficients of the first KL modes. All components of velocity are observed at sparse locations. The observation errors are assumed independent and hence the covariance matrix is diagonal. The error for each observation is constructed using a relative error of of the true value of velocity plus an absolute error of .
For the channel case two observations of velocity are used at locations . The inversion was done twice using both the original covariance matrix and the modified covariance matrix that enforces the boundary conditions. The inferred streamwise velocity and eddy viscosity fields for both cases are shown in Figure 6, and the discrepancies are presented in Table 2. In both cases the velocity is improved nearly everywhere. However, the predicted eddy viscosity does not show the same level of improvement over the baseline solution, and are actually significantly worst when considering the qualitative shape. The good match between velocities in spite of the mismatch in eddy viscosity is due to the illposedness of the problem, i.e. the same velocity field can come from different eddy viscosity fields. This situation has been common in previous research. Although at first the two inferred eddy viscosities might seem similarly wrong, the inferred eddy viscosity using the modified covariance matrix satisfies not only the positivity constraint but all boundary conditions. This is not the case for the inferred eddy viscosity using the original covariance matrix. This can be seen at the symmetry plane location () where the inferred eddy viscosity has zero gradient when using the modified covariance but not when using the original covariance. This means that the inferred eddy viscosity using the modified covariance lies within a smaller subspace of possible solutions that includes more of the physics. The inferred field could be further improved by using regularization to enforce expected qualities such as smoothness, but it should be noted that the boundary conditions are hard physical constraints that should be enforced regardless of any additional regularization.
Base  Inferred (Original)  Inferred (Modified)  

63.8%  43.9%  48.8%  
5.1%  1.3%  1.7% 
For the periodic hills case a single observation of velocity was used at location . The inferred eddy viscosity field is shown in Figure 6(c) and profiles of both the inferred eddy viscosity and the propagated velocity are shown in Figures 6(a) and 6(b), respectively. The discrepancy of the baseline and inferred fields with the truth, using Equation (22), are presented in Table 3. Again the results show improvements in the predicted velocity nearly everywhere but the predicted latent field does not show the same level of improvement. However, the inferred eddy viscosity does satisfy all boundary conditions, including periodicity, as desired.

Base  Inferred  

4 Conclusions
Problems in computational physics often consist of a forward model that solves a set of partial differential equations and associated boundary conditions for the output fields of interest. The forward model typically has latent physical fields that need to be estimated. For instance, the Reynoldsaveraged Navier–Stokes equations describing the mean velocities and pressure of fluid flows contain the unclosed Reynolds stress field. In the inverse problem these latent fields can be treated as model parameters to be inferred. The inverse problem then consists of inferring these latent fields from sparse observations of the output fields. As is the case with the Reynolds stress, these latent fields often have clear physical meaning and constraints. However, the inherent illposedness of inverse problems means that numerous possible latent fields can result in improvements in the fields of interest and agreement with observations. Since the focus is usually on the fields of interest, improvements have often been obtained with unphysical latent fields. Particularly, the boundary conditions on these physical latent fields have been often ignored.
In this paper we demonstrate how to enforce known boundary conditions on the latent fields in Bayesian inversion of physical fields. This is done by appropriate choice of the statistical model used to represent the random latent fields. While this does not necessarily result in the true latent field (i.e. the problem is still illpossed) it does ensure the inferred latent fields satisfy these important physical constraints. In particular, we start with a choice of statistical model consisting of a Gaussian process or a function of a Gaussian process with an initial covariance matrix and modify this matrix in a way that ensures all realizations of the field satisfy the boundary conditions. We demonstrate how to enforce Dirichlet, Neumann, and mixed boundary conditions in this fashion. Particular attention is given to the lognormal process which is a common choice of statistical model for scalar fields with positivity constraint. Some idiosincracies of the lognormal model are explored, particularly the implications of the singularity at boundaries with zerovalue Dirichlet conditions. One consequence of this singularity is that observations should be avoided at locations near such boundaries. Additionally we demonstrate how to enforce periodic boundary conditions by choice of the initial covariance kernel and how to enforce internal observations of the latent field by modifying not just the covariance but also the mean of the Gaussian process.
As a test case, an iterative ensemble Kalman method was used as the Bayesian framework, and the Reynoldsaveraged Navier–Stokes equations describing the mean velocity and pressure of fluid flows were used as the forward model. The inverse problem was solved for the eddy viscosity field from sparse observations of the velocity. The method was tested on two different flows: a onedimensional channel flow and twodimensional flow over periodic hills. The results were compared to Bayesian inversion using the initial unmodified covariance, which does not enforce the boundary conditions. This comparison shows that the method is successful in enforcing boundary conditions on the inferred latent field while still resulting in similar improvements on the fields of interest and similar agreement with observations.
Appendix A Framework for Bayesian Field Inversion
This appendix presents a summary of the Bayesian field inversion framework used in this paper. The framework is adopted from Xiao et al. xiao_quantifying_2016 and is based on the iterative ensemble Kalman method developed by Iglesias et al. iglesias_ensemble_2013 . Ensemble methods are used in practice to avoid calculating the derivatives (e.g. via adjoint methods) of complex computational physical models and use a finite number of samples to represent the distribution of stochastic fields. The ensemble Kalman method is presented first followed by the iterative ensemble Kalman method used in this work.
In the ensemble Kalman method the distribution for the state vector describing the field is represented by a finite number of samples. The samples of the prior distribution are denoted . The observations consist of observation values and observation error (covariance matrix). The distribution representing the observations is also represented using samples
obtained from the multivariate Gaussian distribution
. The operator is the general mapping from state space to observation space and is linearized as , i.e.(23) 
Given these observations each of the prior samples is updated as
(24) 
to obtain the posterior distribution described by the samples . The Kalman gain matrix is given by
(25a)  
(25b) 
where is the sample covariance, is the matrix of all samples, and is the mean over all samples. The formulation in Equation (25b) is obtained by using the definition of the sample covariance and avoids creating the matrix by directly using the values of the samples in observation space . The values can often be obtained in a computationally efficient way avoiding performing this large matrix multiplication.
The ensemble Kalman method assumes linearity in the mapping from state space to observation space as shown in Equation (23). In the problems considered here the mapping from state space to observation space consists of two steps: (1) a set of PDEs for the output fields denoted by the operator , and (2) an observation operator that maps from the output fields to the observation space. That is,
(26) 
We can avoid linearization of if the model is nonlinear and use directly in Equations (24) and (25b) in place of . These values are obtained as in Equation (26) or with the linearziation of the operator as
(27) 
In this case an iterative approach is used where the posterior samples are used as the new prior samples and the observation data ensemble is resampled each iteration. This is done until the mean value of the posterior samples converges. This process is summarized as
(28) 
with , where the subscript indicates the iteration step. The Kalman gain at each iteration is given as
(29) 
where a tilde indicates mean subtracted values, i.e.
and the operator operating on a matrix is the matrix of operating on each column. Note that each iteration requires evaluations of the model. One disadvantage of this framework is the collapse of the sample variance, which results in an accurate posterior mean but loses any information on the posterior variance.
Acknowledgment
This material is based on research sponsored by the U.S. Air Force under agreement number FA86501922204. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon.
References
 (1) M. C. Kennedy, A. O’Hagan, Bayesian calibration of computer models, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63 (3) (2001) 425–464. doi:10.1111/14679868.00294.
 (2) T. A. Oliver, G. Terejanu, C. S. Simmons, R. D. Moser, Validating predictions of unobserved quantities, Computer Methods in Applied Mechanics and Engineering 283 (2015) 1310–1335. doi:10.1016/j.cma.2014.08.023.
 (3) K. Duraisamy, G. Iaccarino, H. Xiao, Turbulence modeling in the age of data, Annual Review of Fluid Mechanics 51 (1) (2019) 357–377. doi:10.1146/annurevfluid010518040547.
 (4) H. Xiao, P. Cinnella, Quantification of model uncertainty in RANS simulations: A review, Progress in Aerospace Sciences 108 (2019) 1–31.
 (5) H. Xiao, J.L. Wu, J.X. Wang, R. Sun, C. Roy, Quantifying and reducing modelform uncertainties in Reynoldsaveraged Navier–Stokes simulations: A datadriven, physicsinformed Bayesian approach, Journal of Computational Physics 324 (2016) 115–136. doi:10.1016/j.jcp.2016.07.038.
 (6) E. Dow, Q. Wang, Quantification of structural uncertainties in the – turbulence model, in: 52nd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference 19th AIAA/ASME/AHS Adaptive Structures Conference 13t, 2011, p. 1762.
 (7) A. P. Singh, K. Duraisamy, Using field inversion to quantify functional errors in turbulence closures, Physics of Fluids 28 (4) (2016) 045110. doi:10.1063/1.4947045.
 (8) S. H. Cheung, T. A. Oliver, E. E. Prudencio, S. Prudhomme, R. D. Moser, Bayesian uncertainty analysis with applications to turbulence modeling, Reliability Engineering & System Safety 96 (9) (2011) 1137–1149. doi:10.1016/j.ress.2010.09.013.
 (9) M. A. Iglesias, K. J. H. Law, A. M. Stuart, Ensemble Kalman methods for inverse problems, Inverse Problems 29 (4) (2013) 045001. doi:10.1088/02665611/29/4/045001.
 (10) W. Edeling, P. Cinnella, R. P. Dwight, H. Bijl, Bayesian estimates of parameter variability in the – turbulence model, Journal of Computational Physics 258 (2014) 73–94.
 (11) W. Edeling, P. Cinnella, R. P. Dwight, Predictive RANS simulations via Bayesian modelscenario averaging, Journal of Computational Physics 275 (2014) 65–91.
 (12) J. Ray, S. Lefantzi, S. Arunajatesan, L. Dechant, Bayesian parameter estimation of a – model for accurate jetincrossflow simulations, AIAA Journal 54 (8) (2016) 2432–2448.
 (13) J. Ray, L. Dechant, S. Lefantzi, J. Ling, S. Arunajatesan, Robust Bayesian calibration of a – model for compressible jetincrossflow simulations, AIAA Journal 56 (12) (2018) 4893–4909.
 (14) J. Ray, S. Lefantzi, S. Arunajatesan, L. Dechant, Learning an eddy viscosity model using shrinkage and Bayesian calibration: A jetincrossflow case study, ASCEASME Journal of Risk and Uncertainty in Engineering Systems, Part B: Mechanical Engineering 4 (1) (2018) 011001.
 (15) M. Meldi, A. Poux, A reduced order model based on Kalman filtering for sequential data assimilation of turbulent flows, Journal of Computational Physics 347 (2017) 207–234.
 (16) M. Meldi, Augmented prediction of turbulent flows via sequential estimators, Flow, Turbulence and Combustion 101 (2) (2018) 389–412.
 (17) W. N. Edeling, G. Iaccarino, P. Cinnella, Datafree and datadriven RANS predictions with quantified uncertainty, Flow, Turbulence and Combustion 100 (3) (2018) 593–616.

(18)
A. Vollant, G. Balarac, C. Corre, Subgridscale scalar flux modelling based on optimal estimation theory and machinelearning procedures, Journal of Turbulence 18 (9) (2017) 854–878.
 (19) M. Ma, J. Lu, G. Tryggvason, Using statistical learning to close twofluid multiphase flow equations for a simple bubbly system, Physics of Fluids 27 (9) (2015) 092101.
 (20) M. Ma, J. Lu, G. Tryggvason, Using statistical learning to close twofluid multiphase flow equations for bubbly flows in vertical channels, International Journal of Multiphase Flow 85 (2016) 336–347.
 (21) M. N. Gibbs, Bayesian Gaussian processes for regression and classification, Ph.D. thesis, University of Cambridge (1998).

(22)
C. E. Rasmussen, Evaluation of Gaussian processes and other methods for nonlinear regression., Ph.D. thesis, University of Toronto (1999).
 (23) B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, N. De Freitas, Taking the human out of the loop: A review of Bayesian optimization, Proceedings of the IEEE 104 (1) (2015) 148–175.
 (24) J. Hensman, N. Fusi, N. D. Lawrence, Gaussian processes for big data, arXiv preprint arXiv:1309.6835.
 (25) E. Solak, R. Murraysmith, W. E. Leithead, D. J. Leith, C. E. Rasmussen, Derivative observations in Gaussian process models of dynamic systems, in: Advances in Neural Information Processing Systems 15, MIT Press, 2003, pp. 1057–1064.
 (26) C. E. Rasmussen, C. K. I. Williams, Gaussian processes for machine learning, Adaptive computation and machine learning, MIT Press, Cambridge, Mass, 2006, oCLC: ocm61285753.

(27)
S. Särkkä, Linear operators and stochastic partial differential equations in Gaussian process regression, in: International Conference on Artificial Neural Networks, Springer, 2011, pp. 151–158.
 (28) T. Graepel, Solving noisy linear operator equations by Gaussian processes: Application to ordinary and partial differential equations, in: ICML, 2003, pp. 234–241.
 (29) M. Raissi, P. Perdikaris, G. E. Karniadakis, Inferring solutions of differential equations using noisy multifidelity data, Journal of Computational Physics 335 (2017) 736–746.
 (30) M. Raissi, P. Perdikaris, G. E. Karniadakis, Machine learning of linear differential equations using Gaussian processes, Journal of Computational Physics 348 (2017) 683–693.

(31)
F. Dondelinger, D. Husmeier, S. Rogers, M. Filippone, ODE parameter inference using adaptive gradient matching with Gaussian processes, in: Artificial Intelligence and Statistics, 2013, pp. 216–228.
 (32) J.L. Wu, C. MichelénStröfer, H. Xiao, Physicsinformed covariance kernel for modelform uncertainty quantification with application to turbulent flows, Computers & Fluids (2019) 104292.
 (33) J. Wu, J.X. Wang, S. C. Shadden, Adding constraints to Bayesian inverse problems, in: Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 33, 2019, pp. 1666–1673.
 (34) J. Wu, J.X. Wang, S. C. Shadden, Improving the convergence of the itterative ensemble Kalman filter by resampling, arXiv:1910.04247 [math].
 (35) X.L. Zhang, C. MichelénStröfer, H. Xiao, Regularization of ensemble Kalman methods for inverse problems, arXiv preprint arXiv:1910.01292.
 (36) C. Mellen, J. Fröhlich, W. Rodi, Large eddy simulation of the flow over periodic hills, in: 16th IMACS World Congress, 2000, pp. 21–25.
 (37) P. Spalart, S. Allmaras, A oneequation turbulence model for aerodynamic flows, in: 30th Aerospace Sciences Meeting and Exhibit, 1992, p. 439.
 (38) W. Jones, B. E. Launder, The prediction of laminarization with a twoequation model of turbulence, International Journal of Heat and Mass Transfer 15 (2) (1972) 301–314.
 (39) D. C. Wilcox, Turbulence Modeling for CFD, 3rd Edition, DCW Industries, 2006.