{'fraction_non_alphanumeric': 0.03569234475697163, 'fraction_numerical': 0.017538998976088658, 'mean_word_length': 4.342772557600721, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 5, 'lorem ipsum': 0, 'www.': 5, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bi-fidelity Modeling of Uncertain and Partially Unknown Systems using DeepONets

Recent advances in modeling large-scale complex physical systems have shifted research focuses towards data-driven techniques. However, generating datasets by simulating complex systems can require significant computational resources. Similarly, acquiring experimental datasets can prove difficult as well. For these systems, often computationally inexpensive, but in general inaccurate, models, known as the low-fidelity models, are available. In this paper, we propose a bi-fidelity modeling approach for complex physical systems, where we model the discrepancy between the true system's response and low-fidelity response in the presence of a small training dataset from the true system's response using a deep operator network (DeepONet), a neural network architecture suitable for approximating nonlinear operators. We apply the approach to model systems that have parametric uncertainty and are partially unknown. Three numerical examples are used to show the efficacy of the proposed approach to model uncertain and partially unknown complex physical systems.

Recently, several neural network architectures capable of approximating mesh-independent solution operators for the governing equations of physical systems have been proposed [29,30,31,32,33,34,35], drawing inspiration from the universal approximation theorem for operators in Chen and Chen [36]. These architectures overcome two shortcomings of the physics-informed use of neural networks to learn a system's behavior: firstly, they do not depend on the mesh used to generate the training dataset. Secondly, they are solely data-driven and do not require the knowledge of the governing equations. Among these network architectures, the deep operator network (DeepONet) proposed in Lu et al. [29,34] uses two separate networks -one with the external source term as input and the other with the coordinates where the system response is sought as input. The former network is known as the branch network while the latter is known as the trunk network. Outputs from these two networks are then used in a linear regression model to produce the system response. Lu et al. [37] showed how the computational domain can be decomposed into many small subdomains to learn the governing equation with only one solution of the equation. While Cai et al. [38] used DeepONets for a multiphysics problem involving modeling of field variables across multiple scales., Ranade et al. [39] and Priyadarshini et al. [40] employed DeepONets for approximating joint subgrid probability density function in turbulent combustion problems and to model non-equilibrium flow physics respectively. Additionally, Wang et al. [41,42] used the error in satisfying the governing equations to train the DeepONets, with a variational formulation of the governing equations being applied by Goswami et al. [43] to train DeepONets to model brittle fracture. Further comparisons of DeepONet with the iterative neural operator architecture proposed in Li et al. [32] were performed in Kovachki et al. [44] and Lu et al. [45]. The approximation errors in the latter neural operator architectures were subsequently discussed in Kovachki et al. [46], Lanthaler et al. [47], Deng et al. [48], and Marcati and Schwab [49].
The training datasets required to construct the DeepONets for a complex physical system can often be limited due to the computational burden associated with simulating such systems. Similarly, repeating experiments to obtain measurement data can become infeasible. In the presence of a small dataset from the true high-fidelity system, we propose herein a modeling strategy with DeepONets that uses a similar but computationally inexpensive model, known as low-fidelity, that captures many aspects of the true system's behavior. The discrepancy between the response from the true physical system and the low-fidelity model is estimated using DeepONets. The assumption used here is that the modeling of the discrepancy requires only a small training dataset of the true system's response, and as a result, bi-fidelity approaches are often advantageous with scarce training data from the high-fidelity simulations [50,51,52]. In principle, we train a DeepONet that uses the low-fidelity response as the input of the branch network to predict the discrepancy between the low-and high-fidelity data. However, we do not incorporate physics loss or error in satisfying the governing equations (as previously employed in physics informed neural networks) to circumvent potential bias incorporation in partially known engineering systems subject to the presence of structural uncertainty.
In this paper, we employ three numerical examples to show the efficacy of the proposed bi-fidelity approach to modeling the response of an uncertain and partially unknown physical system. The first two examples consider a nonlinear oscillator and a heat transfer problem, respectively, where uncertainty is present in some of the parameters of the governing equations. The low-fidelity model used in these examples ignores the nonlinearity or some of the uncertain components of the governing equations. In the third example, we model power generated in a wind farm, , with the two and three dimensional simulations respectively serving as high and low-fidelity model systems.
The rest of the paper is organized as follows: first, we provide a brief background on DeepONets in Section 2. Then, we present our proposed approach for modeling the discrepancy between a low-fidelity system response and the true response for uncertain and partially unknown systems in Section 3. Three numerical examples involving a dynamical system, a heat transfer problem, and a wind farm application are subsequently used in Section 4 to showcase the efficacy 2 A PREPRINT - AUGUST 22, 2022 of the proposed approach. Finally, we conclude the paper with a brief discussion on future applications of the proposed approach.

Background on DeepONets
Let us assume the governing equation of a physical system in a domain Ω is given by where L is a nonlinear differential operator; x is the location; u is the response of the system; and f is the external source. For a boundary value problem, Dirichlet boundary condition u(x) = u D (x) on x ∈ ∂ Ω D and/or Neumann boundary condition L N (u)(x) = u N (x) on x ∈ ∂ Ω N are also present, where L N is the differential operator for the Neumann boundary condition; ∂ Ω D ∪ ∂ Ω N = ∂ Ω with the boundary of the domain being ∂ Ω; and ∂ Ω D ∩ ∂ Ω N = / 0. For an initial value problem, an initial condition u(x 0 ) = u 0 is used instead, where x 0 is the initial instant. The solution of (1) is given by G( f )(x), where G : X → Y denotes the solution operator. For X and Y, there exist embeddings X → L 2 (D) and Y → L 2 (Ω).
A PREPRINT - AUGUST 11, 2022 2 Background on DeepONets Let us assume the governing equation of a physical system in a domain W is given by where L is a nonlinear differential operator; x is the location; u is the response of the system; and f is the external source. For a boundary value problem, Dirichlet boundary condition u(x) = u D (x) on x 2 ∂ W D and/or Neumann boundary condition L N (u)(x) = u N (x) on x 2 ∂ W N are also present, where L N is the differential operator for the Neumann boundary condition; ∂ W D [ ∂ W N = ∂ W with the boundary of the domain being ∂ W; and ∂ W D \ ∂ W N = / 0. For an initial value problem, an initial condition u(x 0 ) = u 0 is used instead, where x 0 is the initial instant. The solution of (1) is given by G( f )(x), where G : X ! Y denotes the solution operator. For X and Y, there exist embeddings X ,! L 2 (D) and Y ,! L 2 (W).
The DeepONet architecture proposed in Lu et al. [29] approximates the solution operator G( f )(x) as .., f (x m );q q q branch )y i (x;q q q trunk ), where c 0 is a constant bias parameter; c i (·), i = 1,. .., p, are obtained from a neural network, known as the branch network, that has the input function (e.g., external source term) f (·) measured at {x i } m i=1 locations as input; y i (·), i = 1,. .., p, are obtained from a neural network, known as the trunk network, that has the location x as the input; and q q q = [q q q T branch q q q T trunk ] T is the parameter vector containing weights and biases q q q branch and q q q trunk from the branch and trunk networks, respectively, and c 0 . This approximation is based on the universal approximation theorem in Chen and Chen [36] for operators. A schematic of the DeepONet architecture used in this study is shown in Figure 1. The main advantage of using a DeepONet architecture over other standard ones, such as feedforward, recurrent, etc., is that DeepONet provides a mesh-independent approximation of the solution operator since DeepONet uses a set of parameters q q q trained only once with different mesh discretizations [29,32].
Lanthaler et al. [47] showed that the DeepONet can be thought of as a composition of three components. The first component is an encoder E that maps the input function space to m sensor inputs  [29] used in the present study.
The DeepONet architecture proposed in Lu et al. [29] approximates the solution operator G( f )(x) as where c 0 is a constant bias parameter; c i (·), i = 1, . . . , p, are obtained from a neural network, known as the branch network, that has the input function (e.g., external source term) f (·) measured at {x i } m i=1 locations as input; ψ i (·), i = 1, . . . , p, are obtained from a neural network, known as the trunk network, that has the location x as the input; and θ θ θ = [θ θ θ T branch θ θ θ T trunk ] T is the parameter vector containing weights and biases θ θ θ branch and θ θ θ trunk from the branch and trunk networks, respectively, and c 0 . This approximation is based on the universal approximation theorem in Chen and Chen [36] for operators. A schematic of the DeepONet architecture used in this study is shown in Figure 1. The main advantage of using a DeepONet architecture over other standard ones, such as feedforward, recurrent, etc., is that DeepONet provides a mesh-independent approximation of the solution operator since DeepONet uses a set of parameters θ θ θ trained only once with different mesh discretizations [32,29].
Lanthaler et al. [47] showed that the DeepONet can be thought of as a composition of three components. The first component is an encoder E that maps the input function space to m sensor inputs The   3   A PREPRINT -AUGUST 22, 2022  A PREPRINT -AUGUST 11, 2022 Encoder, E Fig. 2. DeepONet components and corresponding errors adapted from Lanthaler et al. [47].
DeepONet is the reconstructor, which maps the output of the branch network to the response u(x), thereby showing dependence on the trunk network outputs {y k } p k=1 . With approximate inverses of the encoder and reconstructor defined as decoder D and projector P, respectively, three errors associated with DeepONet approximation can be defined as encoding error E e = kD E Ik L 2 (µ) , approximation error E a = kA P G Dk L 2 (E # µ) , and reconstruction error E r = kR P Ik L 2 (G # µ) , where I is an identity matrix; µ is the measure associated with f ; E # µ and G # µ are push-forward measures associated with the encoder and operator G, respectively. Figure 2 shows the three components of a DeepONet and associated errors as defined here. The total approximation error associated with a DeepONet depends on these three errors as outlined in Theorem 3.2 in Lanthaler et al. [47].
The parameter vector q q q is estimated by solving the optimization problem where in the training dataset for each of N f realizations of the external source term , which might not be the same as the locations where the external source term is measured, are used. These sensor locations can be random or equispaced inside the domain and kept fixed during training. The optimization problem in (3) can be solved using stochastic gradient descent [53][54][55] to obtain optimal parameter values of the DeepONets. In this paper, we use a popular variant of stochastic gradient descent, known as the Adam algorithm [56], to train DeepONets. Note that in Lu et al. [29], the input function f (x) is sampled from function spaces with the locations x as domain and the external source f as range, such as Gaussian random field and Chebyshev polynomials. However, for practical reasons, in this paper, we use probability distribution of the uncertain variables associated with the problems to sample the input function. However, the number of training data for DeepONets can be limited if the computational cost of the simulation of the physical system or the cost of repeating the experiments is significant. In this study, to address the challenge of training DeepONets in the presence of limited data, we propose a bi-fidelity training approach utilizing a computationally inexpensive low-fidelity model of the physical system. In the next section, we detail this proposed approach and discuss two applications to address parametric and structural uncertainties present in the system.

Methodology
In the present work, we assume a low-fidelity model of the physical system in (1) is given by where u l (x) is the low-fidelity response; L l is the differential operator related to the low-fidelity model; and f l (x) is the external source in the low-fidelity model. Next, we denote the discrepancy between the solution of the low-fidelity 4 Fig. 2. DeepONet components and corresponding errors adapted from Lanthaler et al. [47].
second component, the approximator A, models the relation between N f sensor inputs { f (x i )} m i=1 and the coefficients {c k } p k=1 . The encoder and the approximator together comprise the branch network. The third or the final component of DeepONet is the reconstructor, which maps the output of the branch network to the response u(x), thereby showing dependence on the trunk network outputs {ψ k } p k=1 . With approximate inverses of the encoder and reconstructor defined as decoder D and projector P, respectively, three errors associated with DeepONet approximation can be defined as encoding error E e = D • E − I L 2 (µ) , approximation error E a = A − P • G • D L 2 (E # µ) , and reconstruction error E r = R • P − I L 2 (G # µ) , where I is an identity matrix; µ is the measure associated with f ; E # µ and G # µ are push-forward measures associated with the encoder and operator G, respectively. Figure 2 shows the three components of a DeepONet and associated errors as defined here. The total approximation error associated with a DeepONet depends on these three errors as outlined in Theorem 3.2 in Lanthaler et al. [47].
The parameter vector θ θ θ is estimated by solving the optimization problem where in the training dataset for each of N f realizations of the external source term , which might not be the same as the locations where the external source term is measured, are used. These sensor locations can be random or equispaced inside the domain and kept fixed during training. The optimization problem in (3) can be solved using stochastic gradient descent [53,54,55] to obtain optimal parameter values of the DeepONets. In this paper, we use a popular variant of stochastic gradient descent, known as the Adam algorithm [56], to train DeepONets. Note that in Lu et al. [29], the input function f (x) is sampled from function spaces with the locations x as domain and the external source f as range, such as Gaussian random field and Chebyshev polynomials. However, for practical reasons, in this paper, we use probability distribution of the uncertain variables associated with the problems to sample the input function. However, the number of training data for DeepONets can be limited if the computational cost of the simulation of the physical system or the cost of repeating the experiments is significant. In this study, to address the challenge of training DeepONets in the presence of limited data, we propose a bi-fidelity training approach utilizing a computationally inexpensive low-fidelity model of the physical system. In the next section, we detail this proposed approach and discuss two applications to address parametric and structural uncertainties present in the system.

Methodology
In the present work, we assume a low-fidelity model of the physical system in (1) is given by where u l (x) is the low-fidelity response; L l is the differential operator related to the low-fidelity model; and f l (x) is the external source in the low-fidelity model. Next, we denote the discrepancy between the solution of the low-fidelity model u l (x) and that of the physical system u(x) as As it is often the case that the low-fidelity model will capture many important aspects of the physical system's response, we assume that the training of a DeepONet for modeling u d (x) is easier than modeling the true response u and that generating the low-fidelity response u l (x) is computationally inexpensive. We present two applications of the proposed approach next. Note that u d (x) is estimated from the DeepONet and u l (x) from the low-fidelity model for any new x ∈ Ω to predict the response u(x). Further, this new location x need not be the same locations used during training or the places where the external source is measured due to the use of DeepONet to model u d (x). We also do not use another DeepONet or neural network to model the low-fidelity response as the simulation of the low-fidelity model is computationally inexpensive, and this avoids any approximation error introduced through a low-fidelity trained network.
A PREPRINT -AUGUST 11, 2022 model u l (x) and that of the physical system u(x) as As it is often the case that the low-fidelity model will capture many important aspects of the physical system's response, we assume that the training of a DeepONet for modeling u d (x) is easier than modeling the true response u and that generating the low-fidelity response u l (x) is computationally inexpensive. We present two applications of the proposed approach next. Note that u d (x) is estimated from the DeepONet and u l (x) from the low-fidelity model for any new x 2 W to predict the response u(x). Further, this new location x need not be the same locations used during training or the places where the external source is measured due to the use of DeepONet to model u d (x). We also do not use another DeepONet or neural network to model the low-fidelity response as the simulation of the low-fidelity model is computationally inexpensive, and this avoids any approximation error introduced through a low-fidelity trained network.

Branch Network
Trunk Network

Application I: Systems with Parametric Uncertainty
The proposed approach can be used for modeling systems with parametric uncertainty, where a computationally inexpensive low-fidelity model of a similar system is available. For example, consider a system described by where x x x is the vector of uncertain variables; L d and L s are deterministic and uncertain (or stochastic) parts of the differential operator in (1), respectively; the response u depends on the uncertain variables x x x ; and the external source may depend on x x x as well. In this case, a low-fidelity model for (6) can be defined as follows Hence, the response of the physical system is approximated as where u l (x j ;x x x ) m j=1 are used as the input of the branch network and (x,x x x ) as the input of the trunk network to model the discrepancy u d (x;x x x ). A schematic of the approximation in (8) is shown in Figure 3

Application I: Systems with Parametric Uncertainty
The proposed approach can be used for modeling systems with parametric uncertainty, where a computationally inexpensive low-fidelity model of a similar system is available. For example, consider a system described by where ξ ξ ξ is the vector of uncertain variables; L d and L s are deterministic and uncertain (or stochastic) parts of the differential operator in (1), respectively; the response u depends on the uncertain variables ξ ξ ξ ; and the external source may depend on ξ ξ ξ as well. In this case, a low-fidelity model for (6) can be defined as follows Hence, the response of the physical system is approximated as c i (u l (x 1 ;ξ ξ ξ ), . . . , u l (x m ;ξ ξ ξ );θ θ θ branch )ψ i (x,ξ ξ ξ ;θ θ θ trunk ),

5
A PREPRINT - AUGUST 22, 2022 where u l (x j ;ξ ξ ξ ) m j=1 are used as the input of the branch network and (x,ξ ξ ξ ) as the input of the trunk network to model the discrepancy u d (x;ξ ξ ξ ). A schematic of the approximation in (8) is shown in Figure 3.
As a special case, when the operators L d and L s (ξ ξ ξ ) are linear, using (5) in (6) we arrive at which is modeled using a DeepONet, where u l (x) is the solution to (7) and f u (x;ξ ξ ξ ) = f (x;ξ ξ ξ ) − f l (x). Note that we do not necessarily require the linearity of the differential operators for the proposed approach to apply. However, we make that assumption to illustrate the simplified construction of the low-fidelity model and the correction.
A PREPRINT - AUGUST 11, 2022 As a special case, when the operators L d and L s (x x x ) are linear, using (5) in (6) we arrive at which is modeled using a DeepONet, where u l (x) is the solution to (7) and . Note that we do not necessarily require the linearity of the differential operators for the proposed approach to apply. However, we make that assumption to illustrate the simplified construction of the low-fidelity model and the correction.

Branch Network
Trunk Network

Application II: Systems with Structural Uncertainty
In this application, we propose to use (5) to model the response of a partially unknown system. For example, let us assume we can rewrite (1) as follows where L k and L u are the known and unknown parts of the differential operator, respectively; and f k and f u are the known and unknown parts of the external source, respectively. A low-fidelity model for this application can be defined as The response of the physical system is approximated next as .., u l (x m );q q q branch )y i (x;q q q trunk ), where we use u l (x j ) m j=1 as the branch network input and x as the trunk network input. Figure 4 shows a schematic of the approximation in (12).
As a special case, assuming the differential operators L k and L u are linear leads to the following relation for u d (x) in (5) 6 Fig. 4. DeepONet used for bi-fidelity modeling of partially unknown systems in Section 3.2.

Application II: Systems with Structural Uncertainty
In this application, we propose to use (5) to model the response of a partially unknown system. For example, let us assume we can rewrite (1) as follows where L k and L u are the known and unknown parts of the differential operator, respectively; and f k and f u are the known and unknown parts of the external source, respectively. A low-fidelity model for this application can be defined as The response of the physical system is approximated next as where we use u l (x j ) m j=1 as the branch network input and x as the trunk network input. Figure 4 shows a schematic of the approximation in (12). 6 A PREPRINT - AUGUST 22, 2022 As a special case, assuming the differential operators L k and L u are linear leads to the following relation for u d (x) in (5) As before, we do not necessarily require linear differential operators to apply the proposed method. Note that it is straightforward to apply the proposed approach for an uncertain and partially unknown system using (x,ξ ξ ξ ) as trunk network input, as in the previous subsection and addresses both parametric and structural uncertainties. Also, the low-fidelity model can be uncertain in these two applications, but for brevity we do not explicitly write the dependence of the low-fidelity solution on ξ ξ ξ in this section.

Choice of Low-fidelity Model
The choice of the low-fidelity model can have a significant effect on the total DeepONet approximation error. As discussed in Section 2, with a three-component decomposition of the DeepONet, Lanthaler et al. [47] showed that the total error in a DeepONet approximation depends on errors in each of these components [47,Eq. (3.4)]. Among them, the error in the approximator depends on the neural network approximation error, such as in Yarotsky [57,Thm. 1].
On the other hand, the errors in encoder and reconstructor depend on the eigenvalues of covariance operator of the branch network inputs and G( f )(x,ξ ξ ξ ) [47], respectively. In particular, these errors can be bounded from below using . . are eigenvalues of the covariance operator of the branch network input for encoding error E e and eigenvalues of the covariance operator of G( f )(x,ξ ξ ξ ) for reconstruction error E r , respectively.
In the bi-fidelity DeepONet, we are replacing the external source term with the low-fidelity response in the branch network input, and the discrepancy u d (x;ξ ξ ξ ) is modeled using the DeepONet. Hence, an ideal choice for the low-fidelity model that will reduce the encoding error E e will be a low-pass filter system that will remove most of the eigenvalues for i > p of the covariance operator of the branch network input. Also, in many cases, encoding the low-fidelity model's response will be easier than encoding the external source reducing the encoding error. Similarly, to reduce the reconstruction error E r , we can select a low-fidelity model such that the eigenvalues of covariance of u d (x,ξ ξ ξ ) has smaller magnitudes for i > p. However, in practice, we may be limited in terms of choosing the low-fidelity model. A simplified physics or a coarser discretization is often used as the low-fidelity model in most of the multi-fidelity methods. In our numerical examples described in Section 4, we will show that these practical low-fidelity models can still provide significant improvements in the total DeepONet approximation error.
As an elementary example to illustrate the above discussion, we consider a duffing oscillator with governing differential equation where δ is the damping coefficient; α is the linear stiffness coefficient; β is the nonlinear stiffness coefficient; γ is the amplitude of the external force; ω is the frequency of the external force; and u 0 is the initial response of the system. We consider the two applications of bi-fidelity DeepONet discussed in Section 3. In the first application, the parameters β and ω are assumed as uncertain, and a low-fidelity model that does not include the term β u 3 (t;ξ ξ ξ ) is considered. For this case, Figure 5a shows the lower bound of the reconstruction error E r computed using are eigenvalues of the sample covariance matrix of u d (x,ξ ξ ξ ) estimated using 2000 realizations of the random parameters [47,Th. 3.4]. In the second application, we use a low-fidelity model that does not consider the term β u 3 (t;ξ ξ ξ ) and δ is replaced with a different δ nom . The lower bound of the reconstruction error E r computed similarly for this case is shown in Figure 5b. These two plots show that for smaller p in (2), the bi-fidelity DeepONet provides significantly lower bound for the reconstruction error. For large p, this lower bound for both bi-fidelity DeepONet and the DeepONet described in Section 2, denoted as standard in the rest of the paper, reaches zero.
Experimental results for these two applications of bi-fidelity DeepONet are provided in the next section. Note that, we can also choose different low-fidelity models that do not necessarily provide a smaller lower bound for the reconstruction error, but the approximation and encoding errors are smaller with bi-fidelity DeepONet compared to a standard DeepONet as is done in Example 4.2.

Numerical Examples
In this section, we illustrate the proposed approach in approximating the discrepancy between the low-fidelity response and the physical system's response with DeepONets using three numerical examples. The first example uses nonlinear ordinary differential equation for an initial value problem. In the second example, we model a nonlinear heat  and structural uncertainties, where two-dimensional and three-dimensional simulations are used as the low-and highfidelity model of the system, respectively. In these examples, we use two separate datasets, namely a training dataset . We use the training dataset D tr to train the DeepONets, and the accuracy of the trained DeepONets is measured using the validation dataset D val , with the relative validation root mean squared error (RMSE), e val , defined as follows where k·k 2 is the`2-norm of its argument, u pred is the vector of predicted response using DeepONets, and u val is from the validation dataset D val . Note that we select the number of hidden layers in the trunk and branch networks and the number of neurons per hidden layer using an iterative procedure, where we increase the number of neurons per hidden layer gradually up to an upper limit (100 in this study) and then add hidden layers to the network until we observe no reduction of e val . For the number of terms P in (8) and (12), we increase P in multiples of 10 to a maximum, which is the number of neurons in the hidden layers, to observe any reduction of e val .

Example 1: Duffing Oscillator
In this example, we use a duffing oscillator described by the governing differential equation (14). We use a neural network that has three hidden layers with 40 neurons each as the branch network and a neural network that has two hidden layers with 40 neurons each as the trunk network. We use ELU as the activation function in these networks [58]. The number of terms in the summation in (2) is fixed at p = 20 following the procedure discussed in the previous section.  and structural uncertainties, where two-dimensional and three-dimensional simulations are used as the low-and highfidelity model of the system, respectively. In these examples, we use two separate datasets, namely a training dataset . We use the training dataset D tr to train the DeepONets, and the accuracy of the trained DeepONets is measured using the validation dataset D val , with the relative validation root mean squared error (RMSE), e val , defined as follows where k·k 2 is the`2-norm of its argument, u pred is the vector of predicted response using DeepONets, and u val is from the validation dataset D val . Note that we select the number of hidden layers in the trunk and branch networks and the number of neurons per hidden layer using an iterative procedure, where we increase the number of neurons per hidden layer gradually up to an upper limit (100 in this study) and then add hidden layers to the network until we observe no reduction of e val . For the number of terms P in (8) and (12), we increase P in multiples of 10 to a maximum, which is the number of neurons in the hidden layers, to observe any reduction of e val .

Example 1: Duffing Oscillator
In this example, we use a duffing oscillator described by the governing differential equation (14). We use a neural network that has three hidden layers with 40 neurons each as the branch network and a neural network that has two hidden layers with 40 neurons each as the trunk network. We use ELU as the activation function in these networks [58]. The number of terms in the summation in (2) is fixed at p = 20 following the procedure discussed in the previous section.  transfer in a thin plate, where a partial differential equation for the boundary problem is used. We consider both parametric and structural uncertainties to illustrate the applications I and II of the proposed approach from Section 3 in these two examples. The third example considers flow through a wind farm with a combination of parametric and structural uncertainties, where two-dimensional and three-dimensional simulations are used as the low-and highfidelity model of the system, respectively. In these examples, we use two separate datasets, namely a training dataset . We use the training dataset D tr to train the DeepONets, and the accuracy of the trained DeepONets is measured using the validation dataset D val , with the relative validation root mean squared error (RMSE), ε val , defined as follows where · 2 is the 2 -norm of its argument, u pred is the vector of predicted response using DeepONets, and u val is from the validation dataset D val . Note that we select the number of hidden layers in the trunk and branch networks and the number of neurons per hidden layer using an iterative procedure, where we increase the number of neurons per hidden layer gradually up to an upper limit (100 in this study) and then add hidden layers to the network until we observe no reduction of ε val . For the number of terms P in (8) and (12), we increase P in multiples of 10 to a maximum, which is the number of neurons in the hidden layers, to observe any reduction of ε val .

Example 1: Duffing Oscillator
In this example, we use a duffing oscillator described by the governing differential equation (14). We use a neural network that has three hidden layers with 40 neurons each as the branch network and a neural network that has two hidden layers with 40 neurons each as the trunk network. We use ELU as the activation function in these networks [58]. The number of terms in the summation in (2) is fixed at p = 20 following the procedure discussed in the previous section.

Application I: Parametric Uncertainty
For the oscillator defined in (14), we use α = 1.0, δ = 1.0, and γ = 2.5. The uncertain variables are assumed as ] denotes a uniform distribution between limits a and b. We use a training dataset D tr with N tr = 200 and a validation dataset D val with N val = 1000 realizations of the uncertain variables uniformly sampled from their corresponding distributions, respectively, where the response is measured at 100 equidistant points along t. A low-fidelity model of the system can be specified without the nonlinear stiffness, as follows where u l (t;ξ ξ ξ ) is the low-fidelity response. Note that uncertainty in the low-fidelity model only comes through the uncertain variable ω. Figure 6a depicts a representative realization of u(t;ξ ξ ξ ), the low-fidelity response u l (t;ξ ξ ξ ), and the discrepancy u d (t;ξ ξ ξ ). Note that in this example the cost of generating the low-and high-fidelity datasets is insignificant.

Application I: Parametric Uncertainty
For the oscillator defined in (14), we use a = 1.0, d = 1.0, and g = 2.5. The uncertain variables are assumed as ] denotes a uniform distribution between limits a and b. We use a training dataset D tr with N tr = 200 and a validation dataset D val with N val = 1000 realizations of the uncertain variables uniformly sampled from their corresponding distributions, respectively, where the response is measured at 100 equidistant points along t. A low-fidelity model of the system can be specified without the nonlinear stiffness, as follows where u l (t;x x x ) is the low-fidelity response. Note that uncertainty in the low-fidelity model only comes through the uncertain variable w. Figure 6a depicts a representative realization of u(t;x x x ), the low-fidelity response u l (t;x x x ), and the discrepancy u d (t;x x x ). Further, Figure 6b shows the correlation between the true and low-fidelity responses estimated using 500 realizations of the uncertain variables. Note that in this example the cost of generating the low-and high-fidelity datasets is insignificant. Here, a standard use of DeepONet as described in Section 2 is denoted as 'Standard', and the proposed use of DeepONet to model the discrepancy between the low-fidelity model and true system is denoted as 'Bi-fidelity.' The proposed use of DeepONet to model the discrepancy u d (t;x x x ) between the low-fidelity and the true responses produces 4.72 ⇥ 10 3 as the validation error e val with the validation dataset D val . In contrast, the standard DeepONet results in a validation error of 1.11 ⇥ 10 2 . Figure 6c compares the histograms of the relative validation RMSE of these two approaches. Hence, in this case, the proposed approach improves the validation error by more than a factor of two, as modeling of the discrepancy u d (t;x x x ) proved to be easier using the same DeepONet architecture than modeling the response u(t;x x x ).

9
(a) Comparison of the low-fidelity response u l (t;ξ ξ ξ ), correction u d (t;ξ ξ ξ ), and true response u(t;ξ ξ ξ ) for one realization of the uncertain variables l uncertain variable w. Figure 6a depicts a representative realization of u(t;x x x ), the low-fidelity response u l ( discrepancy u d (t;x x x ). Further, Figure 6b shows the correlation between the true and low-fidelity respon using 500 realizations of the uncertain variables. Note that in this example the cost of generating high-fidelity datasets is insignificant.   Here, a sta DeepONet as described in Section 2 is denoted as 'Standard', and the proposed use of DeepONet discrepancy between the low-fidelity model and true system is denoted as 'Bi-fidelity.' The proposed use of DeepONet to model the discrepancy u d (t;x x x ) between the low-fidelity and the tr produces 4.72 ⇥ 10 3 as the validation error e val with the validation dataset D val . In contrast, the standa results in a validation error of 1.11 ⇥ 10 2 . Figure 6c compares the histograms of the relative validation R two approaches. Hence, in this case, the proposed approach improves the validation error by more than a as modeling of the discrepancy u d (t;x x x ) proved to be easier using the same DeepONet architecture than response u(t;x x x ).

9
(b) Relative validation error using DeepONets Fig. 6. The use of DeepONets to model the response of an uncertain Duffing oscillator. Here, a standard use of DeepONet as described in Section 2 is denoted as 'Standard', and the proposed use of DeepONet to model the discrepancy between the low-fidelity model and true system is denoted as 'Bi-fidelity.' The proposed use of DeepONet to model the discrepancy u d (t;ξ ξ ξ ) between the low-fidelity and the true responses produces 4.72 × 10 −3 as the validation error ε val with the validation dataset D val . In contrast, the standard DeepONet results in a validation error of 1.11 × 10 −2 . Figure 6b compares the histograms of the relative validation RMSE of these two approaches. Hence, in this case, the proposed approach improves the validation error by more than a factor of two, as modeling of the discrepancy u d (t;ξ ξ ξ ) proved to be easier using the same DeepONet architecture than modeling the response u(t;ξ ξ ξ ).

Application II: Partially Unknown System (Structural Uncertainty)
In this case, we assume α = 1.0, δ = 1.0, and γ = 2.5 as before, but β = 1.0. The frequency of the external load is assumed to be uniformly distributed as before. We use a training dataset D tr with N tr = 200, where we also add 5% additive Gaussian white noise to the dataset, and a validation dataset D val with N val = 1000. The low-fidelity model is assumed as where the nominal value of the damping coefficient is δ nom = 1.5. Hence, the unknown dynamics not modeled in the low-fidelity model can be given by

Application II: Partially Unknown System (Structural Uncertainty)
In this case, we assume a = 1.0, d = 1.0, and g = 2.5 as before, but b = 1.0. The frequency of the external load is assumed to be uniformly distributed as before. We use a training dataset D tr with N tr = 200, where we also add 5% additive Gaussian white noise to the dataset, and a validation dataset D val with N val = 1000. The low-fidelity model is assumed as where the nominal value of the damping coefficient is d nom = 1.5. Hence, the unknown dynamics not modeled in the low-fidelity model can be given by Figure 7a shows a comparison of the response from the low-fidelity model u l (t;x x x ), true response u(t; x ), and their difference u d (t; x ). When DeepONet is used to model the discrepancy u d (t, x ) between the true response and the low-fidelity model response, it produces a validation error e val of 9.98 ⇥ 10 3 . Note that the standard use of DeepONet produces an e val of 2.22 ⇥ 10 2 with the same training and validation datasets. Figure 7b shows a comparison of histograms from the proposed use of DeepONet to model u d (t; x ) and the standard use of DeepONet as described in Section 2. The figure shows that, in this case as well, modeling u d (t, x ) as a function of the low-fidelity response u l (t;x x x ) using DeepONet can achieve better accuracy compared to the standard use. Hence, this example shows the effectiveness of the proposed approach to model the discrepancy using DeepONet for uncertain and partially unknown dynamical systems, which reduces the validation error by a factor of more than two.

Example 2: Nonlinear Steady-state Heat Transfer in a Thin Plate
In our second example, we use a two-dimensional steady-state heat transfer problem in a thin L-shaped plate with dimensions shown in Figure 8a. The left boundary is assumed to have a constant temperature of 100 K. The steady-state temperature at any location inside the plate, T (x), is given by where -: ; h t is the thickness of the plate; k(x) is the thermal conductivity of the plate; Q c (x) is the amount of heat transferred through one face of the plate per unit area using convection; Q r (x) is the amount of heat transferred through one face of the plate per unit area using radiation; and Q i (x) is the amount of heat added to the plate per unit area from an internal heat source. Heat transferred by convection can be given by Newton's law of cooling as where h c is the specified convection coefficient and T a is the ambient temperature. For heat transferred due to radiation, the Stefan-Boltzmann law can be used,

Application II: Partially Unknown System (Structural Uncertainty)
In this case, we assume a = 1.0, d = 1.0, and g = 2.5 as before, but b = 1.0. The frequency of the external load is assumed to be uniformly distributed as before. We use a training dataset D tr with N tr = 200, where we also add 5% additive Gaussian white noise to the dataset, and a validation dataset D val with N val = 1000. The low-fidelity model is assumed as where the nominal value of the damping coefficient is d nom = 1.5. Hence, the unknown dynamics not modeled in the low-fidelity model can be given by Figure 7a shows a comparison of the response from the low-fidelity model u l (t;x x x ), true response u(t; x ), and their difference u d (t; x ). When DeepONet is used to model the discrepancy u d (t, x ) between the true response and the low-fidelity model response, it produces a validation error e val of 9.98 ⇥ 10 3 . Note that the standard use of DeepONet produces an e val of 2.22 ⇥ 10 2 with the same training and validation datasets. Figure 7b shows a comparison of histograms from the proposed use of DeepONet to model u d (t; x ) and the standard use of DeepONet as described in Section 2. The figure shows that, in this case as well, modeling u d (t, x ) as a function of the low-fidelity response u l (t;x x x ) using DeepONet can achieve better accuracy compared to the standard use. Hence, this example shows the effectiveness of the proposed approach to model the discrepancy using DeepONet for uncertain and partially unknown dynamical systems, which reduces the validation error by a factor of more than two.

Example 2: Nonlinear Steady-state Heat Transfer in a Thin Plate
In our second example, we use a two-dimensional steady-state heat transfer problem in a thin L-shaped plate with dimensions shown in Figure 8a. The left boundary is assumed to have a constant temperature of 100 K. The steady-state temperature at any location inside the plate, T (x), is given by where -: ; h t is the thickness of the plate; k(x) is the thermal conductivity of the plate; Q c (x) is the amount of heat transferred through one face of the plate per unit area using convection; Q r (x) is the amount of heat transferred through one face of the plate per unit area using radiation; and Q i (x) is the amount of heat added to the plate per unit area from an internal heat source. Heat transferred by convection can be given by Newton's law of cooling as where h c is the specified convection coefficient and T a is the ambient temperature. For heat transferred due to radiation, the Stefan-Boltzmann law can be used, 10 (b) Relative validation error using DeepONets produces an ε val of 2.22 × 10 −2 with the same training and validation datasets. Figure 7b shows a comparison of histograms from the proposed use of DeepONet to model u d (t; ξ ) and the standard use of DeepONet as described in Section 2. The figure shows that, in this case as well, modeling u d (t, ξ ) as a function of the low-fidelity response u l (t;ξ ξ ξ ) using DeepONet can achieve better accuracy compared to the standard use. Hence, this example shows the effectiveness of the proposed approach to model the discrepancy using DeepONet for uncertain and partially unknown dynamical systems, which reduces the validation error by a factor of more than two.

Example 2: Nonlinear Steady-state Heat Transfer in a Thin Plate
In our second example, we use a two-dimensional steady-state heat transfer problem in a thin L-shaped plate with dimensions shown in Figure 8a. The left boundary is assumed to have a constant temperature of 100 K. The steady-state temperature at any location inside the plate, T (x), is given by where ∇ : ; h t is the thickness of the plate; k(x) is the thermal conductivity of the plate; Q c (x) is the amount of heat transferred through one face of the plate per unit area using convection; Q r (x) is the amount of heat transferred through one face of the plate per unit area using radiation; and Q i (x) is the amount of heat added to the plate per unit area from an internal heat source. Heat transferred by convection can be given by Newton's law of cooling as where h c is the specified convection coefficient and T a is the ambient temperature. For heat transferred due to radiation, the Stefan-Boltzmann law can be used, where ε is the emissivity of the plate surface and σ is the Stefan-Boltzmann constant of the plate. Given (20) and (21), we rewrite (19) as In this example, we set h t = 0.01 m; h c = 1 W(m 2 K) −1 ; ε = 0.025; and σ = 5.670373 × 10 −8 W(m 2 K 4 ) −1 (similar to copper) to generate the datasets. We also assume the internal heat source to be an uncertain squared exponential source centered at the (0, 0) corner given by where Q 0 = 5 × 10 3 Wm −2 and ξ Q is a standard Gaussian random variable. We use the finite element method with the mesh shown in Figure 8a to solve (22). Figure 8b shows a realization of the steady-state temperature in the plate. In this example, we use a neural network that has three hidden layers with 40 neurons each and ELU activation as the branch network, and a neural network that has two hidden layers with 40 neurons each and ELU activation as the trunk network. The number of terms in the summation in (2) is fixed at p = 20 to produced the smallest validation error.
where Q 0 = 5 ⇥ 10 3 Wm 2 and x Q is a standard Gaussian random variable. We use the finite element method with the mesh shown in Figure 8a to solve (22). Figure 8b shows a realization of the steady-state temperature in the plate. In this example, we use a neural network that has three hidden layers with 40 neurons each and ELU activation as the branch network, and a neural network that has two hidden layers with 40 neurons each and ELU activation as the trunk network. The number of terms in the summation in (2) is fixed at p = 20 to produced the smallest validation error.

Application I: Uncertain Thermal Conductivity
In our first application, we model the thermal conductivity as where k 0 = 400 W(m · K) 1 and z(x) is a zero mean Gaussian random field with a square exponential covariance function K(x, x 0 ) where l i , i = 1, 2, are length scale parameters, which we choose as 0.5 for both dimensions. A Karhunen-Loève (KL) expansion truncated at the 25th term that captures 97% of the energy (defined as the sum of all eigenvalues) is used to generate realizations of the random field z(x). For the low-fidelity model, we consider where T l is the temperature given by the low-fidelity model; and Q nom i (x) = Q 0 exp Â 2 i=1 x 2 i . Hence, the low-fidelity model is deterministic and we only need to solve it once. Here, the input to the trunk network of the bi-fidelity DeepONet that models the discrepancy (see Figure 3) contains 25 standard normal random variables used in the KL expansion of z(x), x Q in (23), and the location x for predicting the temperature. We use N tr = 30 realizations of the 11 (a) Finite element mesh used where Q 0 = 5 ⇥ 10 3 Wm 2 and x Q is a standard Gaussian random variable. We use the finite element method with the mesh shown in Figure 8a to solve (22). Figure 8b shows a realization of the steady-state temperature in the plate. In this example, we use a neural network that has three hidden layers with 40 neurons each and ELU activation as the branch network, and a neural network that has two hidden layers with 40 neurons each and ELU activation as the trunk network. The number of terms in the summation in (2) is fixed at p = 20 to produced the smallest validation error.

Application I: Uncertain Thermal Conductivity
In our first application, we model the thermal conductivity as where k 0 = 400 W(m · K) 1 and z(x) is a zero mean Gaussian random field with a square exponential covariance function K(x, x 0 ) where l i , i = 1, 2, are length scale parameters, which we choose as 0.5 for both dimensions. A Karhunen-Loève (KL) expansion truncated at the 25th term that captures 97% of the energy (defined as the sum of all eigenvalues) is used to generate realizations of the random field z(x). For the low-fidelity model, we consider where T l is the temperature given by the low-fidelity model; and Q nom i (x) = Q 0 exp Â 2 i=1 x 2 i . Hence, the low-fidelity model is deterministic and we only need to solve it once. Here, the input to the trunk network of the bi-fidelity DeepONet that models the discrepancy (see Figure 3) contains 25 standard normal random variables used in the KL expansion of z(x), x Q in (23), and the location x for predicting the temperature. We use N tr = 30 realizations of the 11 (b) Temperature for one realization of the uncertain variables Fig. 8. Nonlinear steady-state heat transfer in a thin plate is studied in Example 2.

Application I: Uncertain Thermal Conductivity
In our first application, we model the thermal conductivity as where k 0 = 400 W(m · K) −1 and z(x) is a zero mean Gaussian random field with a square exponential covariance function K(x, x ) where l i , i = 1, 2, are length scale parameters, which we choose as 0.5 for both dimensions. A Karhunen-Loève (KL) expansion truncated at the 25th term that captures 97% of the energy (defined as the sum of all eigenvalues) is used to generate realizations of the random field z(x). For the low-fidelity model, we consider where T l is the temperature given by the low-fidelity model; and Q nom Hence, the low-fidelity model is deterministic and we only need to solve it once. Here, the input to the trunk network of the bi-fidelity DeepONet that models the discrepancy (see Figure 3) contains 25 standard normal random variables used in the KL expansion of z(x), ξ Q in (23), and the location x for predicting the temperature. We use N tr = 30 realizations of the uncertain variables in the training dataset D tr and N val = 1000 realizations for the validation dataset D val . For each of these realizations, we measure the temperature of the plate at 200 locations marked using circles and crosses in Figure  9a. During training, we use the temperature at 100 of these locations marked with circles and validate the DeepONets by predicting the temperature at the remaining 100 locations marked with crosses. We compare the result to the standard use of DeepONet as discussed in Section 2, where a sample of Q i (x) and the right hand side of (22), measured at the training locations, marked with circles in Figure 9a, are used as input to the branch network.
Using the bi-fidelity DeepONet, the validation error ε val , as defined in (15), is 6.72 × 10 −3 . However, the validation error ε val in the standard use of DeepONet is 8.70 × 10 −2 . Figure 9b shows the histograms comparing the errors from these two approaches, denoted as 'standard' and 'bi-fidelity', respectively. The validation error ε val and the histograms show that the proposed approach improves the accuracy of the predictions by almost an order of magnitude. This is primarily due to the fact that modeling of the temperature discrepancy can be performed relatively accurately using DeepONets in the presence of a small training dataset.

Application II: Partially Unknown System (Structural Uncertainty)
Next, we assume that the low-fidelity model is unaware of the radiation component in (22), i.e., Using the bi-fidelity DeepONet, the validation error e val , as defined in (15), is 6.72 ⇥ 10 3 . However, the validation error e val in the standard use of DeepONet is 8.70 ⇥ 10 2 . Figure 9b shows the histograms comparing the errors from these two approaches, denoted as 'standard' and 'bi-fidelity', respectively. The validation error e val and the histograms show that the proposed approach improves the accuracy of the predictions by almost an order of magnitude. This is primarily due to the fact that modeling of the temperature discrepancy can be performed relatively accurately using DeepONets in the presence of a small training dataset.       Using the bi-fidelity DeepONet, the validation error e val , as defined in (15), is 6.72 ⇥ 10 3 . However, the validation error e val in the standard use of DeepONet is 8.70 ⇥ 10 2 . Figure 9b shows the histograms comparing the errors from these two approaches, denoted as 'standard' and 'bi-fidelity', respectively. The validation error e val and the histograms show that the proposed approach improves the accuracy of the predictions by almost an order of magnitude. This is primarily due to the fact that modeling of the temperature discrepancy can be performed relatively accurately using DeepONets in the presence of a small training dataset.        Using the bi-fidelity DeepONet, the validation error e val , as defined in (15), is 6.72 ⇥ 10 3 . However, the validation error e val in the standard use of DeepONet is 8.70 ⇥ 10 2 . Figure 9b shows the histograms comparing the errors from these two approaches, denoted as 'standard' and 'bi-fidelity', respectively. The validation error e val and the histograms show that the proposed approach improves the accuracy of the predictions by almost an order of magnitude. This is primarily due to the fact that modeling of the temperature discrepancy can be performed relatively accurately using DeepONets in the presence of a small training dataset.       Using the bi-fidelity DeepONet, the validation error e val , as defined in (15), is 6.72 ⇥ 10 3 . However, the validation error e val in the standard use of DeepONet is 8.70 ⇥ 10 2 . Figure 9b shows the histograms comparing the errors from these two approaches, denoted as 'standard' and 'bi-fidelity', respectively. The validation error e val and the histograms show that the proposed approach improves the accuracy of the predictions by almost an order of magnitude. This is primarily due to the fact that modeling of the temperature discrepancy can be performed relatively accurately using DeepONets in the presence of a small training dataset.        where Q nom Hence, the low-fidelity model is deterministic and we only need to solve it once. Here, we use x and ξ Q as the input to the trunk network. The training dataset D tr and the validation dataset D val consist of N tr = 3 and N val = 1000 data points, respectively, similar to the previous subsection. For each of these realizations, we measure the temperature of the plate at 200 locations marked using circles and crosses in Figure 10a. During training, we use the temperature at 20 of these locations marked with circles and validate the DeepONets by predicting the temperature at the remaining 180 locations marked with crosses. When using DeepONet for modeling the discrepancy between the low-fidelity model and the true system, the validation error ε val is 5.32 × 10 −3 . The standard use of DeepONet, however, gives ε val as 1.46 × 10 −2 , one order of magnitude larger than the proposed bi-fidelity DeepONet. Figure 10b compares the histograms of the errors from these two approaches, illustrating the advantage of the bi-fidelity approach.

Example 3: Wind Farm
In our third example, we model the power generated in a wind farm and consider uncertain wind speed, inflow direction, and yaw angle. Estimating the expected power with respect to the joint distribution of wind speed and direction is a common calculation in determining the annual energy production of a wind plant; however, the number of function evaluations required often precludes the use of high-fidelity models in industry [59]. Deliberately, yawing wind turbine rotors out of perpendicular to incoming wind is an emerging flow control technique called wake steering, which can improve the overall plant output by directing wind turbine wakes away from downwind turbines [60]. The wake steering optimization problem is fraught with uncertainty due to misalignment in the nacelle-mounted wind vanes located just downwind of the rotor [61]. We consider two layouts of a wind farm with the number of turbines six and 64, respectively. For the probability distribution of uncertain wind speed and direction, we use 10 minute-average measurement data from the Princess Amalia wind farm in the Netherlands. The joint probability density function of the wind speed and direction estimated using these measurements is shown in Figure 11.
During training, we use the temperature at 20 of these locations marked with circles and validate the DeepONets by predicting the temperature at the remaining 180 locations marked with crosses. When using DeepONet for modeling the discrepancy between the low-fidelity model and the true system, the validation error e val is 5.32 ⇥ 10 3 . The standard use of DeepONet, however, gives e val as 1.46 ⇥ 10 2 , one order of magnitude larger than the proposed bi-fidelity DeepONet. Figure 10b compares the histograms of the errors from these two approaches, illustrating the advantage of the bi-fidelity approach.

Example 3: Wind Farm
In our third example, we model the power generated in a wind farm and consider uncertain wind speed, inflow direction, and yaw angle. Estimating the expected power with respect to the joint distribution of wind speed and direction is a common calculation in determining the annual energy production of a wind plant; however, the number of function evaluations required often precludes the use of high-fidelity models in industry [59]. Deliberately, yawing wind turbine rotors out of perpendicular to incoming wind is an emerging flow control technique called wake steering, which can improve the overall plant output by directing wind turbine wakes away from downwind turbines [60]. The wake steering optimization problem is fraught with uncertainty due to misalignment in the nacelle-mounted wind vanes located just downwind of the rotor [61]. We consider two layouts of a wind farm with the number of turbines six and 64, respectively. For the probability distribution of uncertain wind speed and direction, we use 10 minute-average measurement data from the Princess Amalia wind farm in the Netherlands. The joint probability density function of the wind speed and direction estimated using these measurements is shown in Figure 11. For each realization of the wind speed, direction, and yaw angle, we use the WindSE toolkit [62], a Python package with FEniCS [63] as the backend solver, to generate the training and validation datasets. Two-dimensional (2D) Reynolds Averaged Navier-Stokes (RANS) simulations using the finite element mesh shown in Figure 13a are used to generate the low-fidelity data. Three-dimensional (3D) simulations using the mesh shown in Figure 13b are used to estimate the power generated from the true system. For the two-and three-dimensional simulations, a mixing length model [64] [65, 7.1.7] is used for closing the turbulent stress term, with a uniform mixing length l m = 15 m. In the three-dimensional simulations, a velocity log-profile is used for the atmospheric boundary layer with Von-Karman constant k = 0.4. The log profile is such that the hub-height velocity is the value of the sample drawn from the pdf shown in Fig. 11. In the two-dimensional simulations, the governing equations follow the 2D approximation of Boersma et al. [66] to approximate the flow field at hub-height, which is sufficient for estimating the power output. For each realization of the wind speed, direction, and yaw angle, we use the WindSE toolkit [62], a Python package with FEniCS [63] as the backend solver, to generate the training and validation datasets. Two-dimensional (2D) Reynolds Averaged Navier-Stokes (RANS) simulations using the finite element mesh shown in Figure 13a are used to generate the low-fidelity data. Three-dimensional (3D) simulations using the mesh shown in Figure 13b are used to estimate the power generated from the true system. For the two-and three-dimensional simulations, a mixing length model [64] [65, 7.1.7] is used for closing the turbulent stress term, with a uniform mixing length l m = 15 m. In the three-dimensional simulations, a velocity log-profile is used for the atmospheric boundary layer with Von-Karman constant κ = 0.4. The log profile is such that the hub-height velocity is the value of the sample drawn from the pdf shown in Fig. 11. In the two-dimensional simulations, the governing equations follow the 2D approximation of Boersma et al. [66] to approximate the flow field at hub-height, which is sufficient for estimating the power output.
We use a training D tr and validation D val datasets consisting of N tr = 100 and N val = 1000 realizations of the uncertain variables and power generated. The realizations of wind speed and wind direction are generated using the histograms from the dataset for the Princess Amalia wind farm in the Netherlands. We use a neural network with three hidden layers and 100 neurons each with ELU activation as the branch network and another neural network with two hidden layers and 100 neurons each with ELU activation as the trunk network in this example. In (2), we truncate the summation at p = 60 that produces the smallest validation error. Results for two layouts with six and 64 turbines are discussed next.

Layout 1: Six Turbines
First, we model the power generated in a wind farm with six turbines arranged in a 3 × 2 grid as shown in Figure 12a. The wind turbines are placed 474 m apart in the x 1 direction and 387 m apart in the x 2 direction. We use the NREL 5 MW reference turbine [67], which has a 126 m rotor diameter and 90 m hub height. Figure 12b shows a single turbine with yaw angle between the wind direction and the turbine axis, which is used as an input parameter in this example.
We consider three cases here. In Case I, we only consider wind speed and direction as uncertain variables and assume yaw angle to be zero for all the turbines. In Case II, we assume the yaw angle to be uncertain and uniformly distributed 13 A PREPRINT -AUGUST 22, 2022

Layout 1: Six Turbines
First, we model the power generated in a wind farm with six turbines arranged in a 3 ⇥ 2 grid as shown in Figure 12a. The wind turbines are placed 474 m apart in the x 1 direction and 387 m apart in the x 2 direction. We use the NREL 5 MW reference turbine [67], which has a 126 m rotor diameter and 90 m hub height. Figure 12b shows a single turbine with yaw angle between the wind direction and the turbine axis, which is used as an input parameter in this example.

Layout 2: 64 Turbines
In this layout, we consider 64 turbines in total arranged in a 8 ⇥ 8 grid that is more characteristic of the size of a utility-scale wind plant. Other specifications of the turbine layout and the input parameters are the same as before. Figure 12 shows the two-and three-dimensional meshes used to generate the training and validation datasets. The results from the proposed bi-fidelity training of DeepONets are compared with a DeepONet (denoted as 'standard') trained only using a high-fidelity dataset with wind speed, direction, and yaw angle as branch network inputs and the turbine location as the trunk network input. Note that one 2D simulation of this layout takes 3.68 min. (average of 300 runs) compared to 1.15 hrs. (average of 300 runs) for one 3D simulation on the National Renewable Energy Laboratory's Eagle high-performance computing system. As before, for a fair comparison, we train the standard DeepONet with a training dataset Dtr consisting of 106 realizations of the uncertain variables and generated power to keep the cost of generating the training dataset similar to that of the bi-fidelity DeepONet. In this layout, the validation error eval using the proposed approach is 1.68 ⇥ 10 2 , which is about half of 2.97 ⇥ 10 2 , eval obtained using the standard DeepONet. Figure 13 shows the histograms of the validation errors. We consider three cases here. In Case I, we only consider wind speed and direction as uncertain variables and assume yaw angle to be zero for all the turbines. In Case II, we assume the yaw angle to be uncertain and uniformly distributed between [ p/2, p/2], but the same for all the turbines. In the final Case III, we assume each of the six turbines has uncertain yaw angles that are independent and identically distributed uniformly between [ p/2, p/2].
We compare the the results with a standard DeepONet that uses wind speed, direction, and yaw angle (for the last two cases) as branch network input and the turbine location as the trunk network input. Note that one 2D simulation of this layout takes 4.69 s (average of 1000 runs) compared to 3.54 min (average of 1000 runs) for one 3D simulation on the National Renewable Energy Laboratory's Eagle high-performance computing system (https://www.nrel.gov/hpc/eagle- First, we model the power generated in a wind farm with six turbines arranged in a 3 ⇥ 2 grid as shown in Figure 12a. The wind turbines are placed 474 m apart in the x 1 direction and 387 m apart in the x 2 direction. We use the NREL 5 MW reference turbine [67], which has a 126 m rotor diameter and 90 m hub height. Figure 12b shows a single turbine with yaw angle between the wind direction and the turbine axis, which is used as an input parameter in this example.

Layout 2: 64 Turbines
In this layout, we consider 64 turbines in total arranged in a 8 ⇥ 8 grid that is more characteristic of the size of a utility-scale wind plant. Other specifications of the turbine layout and the input parameters are the same as before. Figure 12 shows the two-and three-dimensional meshes used to generate the training and validation datasets. The results from the proposed bi-fidelity training of DeepONets are compared with a DeepONet (denoted as 'standard') trained only using a high-fidelity dataset with wind speed, direction, and yaw angle as branch network inputs and the turbine location as the trunk network input. Note that one 2D simulation of this layout takes 3.68 min. (average of 300 runs) compared to 1.15 hrs. (average of 300 runs) for one 3D simulation on the National Renewable Energy Laboratory's Eagle high-performance computing system. As before, for a fair comparison, we train the standard DeepONet with a training dataset Dtr consisting of 106 realizations of the uncertain variables and generated power to keep the cost of generating the training dataset similar to that of the bi-fidelity DeepONet. In this layout, the validation error eval using the proposed approach is 1.68 ⇥ 10 2 , which is about half of 2.97 ⇥ 10 2 , eval obtained using the standard DeepONet. Figure 13 shows the histograms of the validation errors. We consider three cases here. In Case I, we only consider wind speed and direction as uncertain variables and assume yaw angle to be zero for all the turbines. In Case II, we assume the yaw angle to be uncertain and uniformly distributed between [ p/2, p/2], but the same for all the turbines. In the final Case III, we assume each of the six turbines has uncertain yaw angles that are independent and identically distributed uniformly between [ p/2, p/2].
We compare the the results with a standard DeepONet that uses wind speed, direction, and yaw angle (for the last two cases) as branch network input and the turbine location as the trunk network input. Note that one 2D simulation of this layout takes 4.69 s (average of 1000 runs) compared to 3.54 min (average of 1000 runs) for one 3D simulation on the National Renewable Energy Laboratory's Eagle high-performance computing system (https://www.nrel.gov/hpc/eagle-

Layout 1: Six Turbines
First, we model the power generated in a wind farm with six turbines arranged in a 3 ⇥ 2 grid as shown in Figure 12a. The wind turbines are placed 474 m apart in the x 1 direction and 387 m apart in the x 2 direction. We use the NREL 5 MW reference turbine [67], which has a 126 m rotor diameter and 90 m hub height. Figure 12b shows a single turbine with yaw angle between the wind direction and the turbine axis, which is used as an input parameter in this example.

Layout 2: 64 Turbines
In this layout, we consider 64 turbines in total arranged in a 8 ⇥ 8 grid that is more characteristic of the size of a utility-scale wind plant. Other specifications of the turbine layout and the input parameters are the same as before. Figure 12 shows the two-and three-dimensional meshes used to generate the training and validation datasets. The results from the proposed bi-fidelity training of DeepONets are compared with a DeepONet (denoted as 'standard') trained only using a high-fidelity dataset with wind speed, direction, and yaw angle as branch network inputs and the turbine location as the trunk network input. Note that one 2D simulation of this layout takes 3.68 min. (average of 300 runs) compared to 1.15 hrs. (average of 300 runs) for one 3D simulation on the National Renewable Energy Laboratory's Eagle high-performance computing system. As before, for a fair comparison, we train the standard DeepONet with a training dataset Dtr consisting of 106 realizations of the uncertain variables and generated power to keep the cost of generating the training dataset similar to that of the bi-fidelity DeepONet. In this layout, the validation error eval using the proposed approach is 1.68 ⇥ 10 2 , which is about half of 2.97 ⇥ 10 2 , eval obtained using the standard DeepONet. Figure 13 shows the histograms of the validation errors. We consider three cases here. In Case I, we only consider wind speed and direction as uncertain variables and assume yaw angle to be zero for all the turbines. In Case II, we assume the yaw angle to be uncertain and uniformly distributed between [ p/2, p/2], but the same for all the turbines. In the final Case III, we assume each of the six turbines has uncertain yaw angles that are independent and identically distributed uniformly between [ p/2, p/2].
We compare the the results with a standard DeepONet that uses wind speed, direction, and yaw angle (for the last two cases) as branch network input and the turbine location as the trunk network input. Note that one 2D simulation of this layout takes 4.69 s (average of 1000 runs) compared to 3.54 min (average of 1000 runs) for one 3D simulation on the National Renewable Energy Laboratory's Eagle high-performance computing system (https://www.nrel.gov/hpc/eagle-

Layout 1: Six Turbines
First, we model the power generated in a wind farm with six turbines arranged in a 3 ⇥ 2 grid as shown in Figure 12a. The wind turbines are placed 474 m apart in the x 1 direction and 387 m apart in the x 2 direction. We use the NREL 5 MW reference turbine [67], which has a 126 m rotor diameter and 90 m hub height. Figure 12b shows a single turbine with yaw angle between the wind direction and the turbine axis, which is used as an input parameter in this example.

Layout 2: 64 Turbines
In this layout, we consider 64 turbines in total arranged in a 8 ⇥ 8 grid that is more characteristic of the size of a utility-scale wind plant. Other specifications of the turbine layout and the input parameters are the same as before. Figure 12 shows the two-and three-dimensional meshes used to generate the training and validation datasets. The results from the proposed bi-fidelity training of DeepONets are compared with a DeepONet (denoted as 'standard') trained only using a high-fidelity dataset with wind speed, direction, and yaw angle as branch network inputs and the turbine location as the trunk network input. Note that one 2D simulation of this layout takes 3.68 min. (average of 300 runs) compared to 1.15 hrs. (average of 300 runs) for one 3D simulation on the National Renewable Energy Laboratory's Eagle high-performance computing system. As before, for a fair comparison, we train the standard DeepONet with a training dataset Dtr consisting of 106 realizations of the uncertain variables and generated power to keep the cost of generating the training dataset similar to that of the bi-fidelity DeepONet. In this layout, the validation error eval using the proposed approach is 1.68 ⇥ 10 2 , which is about half of 2.97 ⇥ 10 2 , eval obtained using the standard DeepONet. Figure 13 shows the histograms of the validation errors.   We consider three cases here. In Case I, we only consider wind speed and direction as uncertain variables and assume yaw angle to be zero for all the turbines. In Case II, we assume the yaw angle to be uncertain and uniformly distributed between [ p/2, p/2], but the same for all the turbines. In the final Case III, we assume each of the six turbines has uncertain yaw angles that are independent and identically distributed uniformly between [ p/2, p/2].
We compare the the results with a standard DeepONet that uses wind speed, direction, and yaw angle (for the last two cases) as branch network input and the turbine location as the trunk network input. Note that one 2D simulation of this layout takes 4.69 s (average of 1000 runs) compared to 3.54 min (average of 1000 runs) for one 3D simulation on the National Renewable Energy Laboratory's Eagle high-performance computing system (https://www.nrel.gov/hpc/eagle- between [−π/2, π/2], but the same for all the turbines. In the final Case III, we assume each of the six turbines has uncertain yaw angles that are independent and identically distributed uniformly between [−π/2, π/2]. We compare the the results with a standard DeepONet that uses wind speed, direction, and yaw angle (for the last two cases) as branch network input and the turbine location as the trunk network input. Note that one 2D simulation of this layout takes 4.69 s (average of 1000 runs) compared to 3.54 min (average of 1000 runs) for one 3D simulation on the National Renewable Energy Laboratory's Eagle high-performance computing system (https://www.nrel.gov/hpc/eaglesystem.html). Therefore, for a fair comparison, we train this DeepONet with a training dataset D tr consisting of 103 realizations of the uncertain variables and corresponding power generated (i.e., with similar computational cost of generating the training dataset compared to the bi-fidelity training dataset used in the proposed approach).
We list the results in Table 1, which shows that the proposed use of DeepONet to model the discrepancy between the power outputs from 2D and 3D simulations provides significant improvement in the validation error. Also, the standard use of DeepONet shows improvement when we add yaw angle as another random input. This is primarily due to the fact that yaw angle has a significant effect on the power generated, and the standard DeepONet is able to model that relation with increasing connection due to the addition of a new input node in the network. A similar but slight improvement in the proposed approach can also be noticed in Cases II and III in Table 1.

Layout 2: 64 Turbines
In this layout, we consider 64 turbines in total arranged in a 8 × 8 grid that is more characteristic of the size of a utility-scale wind plant. Other specifications of the turbine layout and the input parameters are the same as before. Figure 14 shows the two-and three-dimensional meshes used to generate the training and validation datasets. The results from the proposed bi-fidelity training of DeepONets are compared with a DeepONet (denoted as 'standard') trained only using a high-fidelity dataset with wind speed, direction, and yaw angle as branch network inputs and the turbine location as the trunk network input. Note that one 2D simulation of this layout takes 3.68 min. (average of 300 runs) compared to 1.15 hrs. (average of 300 runs) for one 3D simulation on the National Renewable Energy Laboratory's Eagle high-performance computing system. As before, for a fair comparison, we train the standard DeepONet with a training dataset D tr consisting of 106 realizations of the uncertain variables and generated power to keep the cost of generating the training dataset similar to that of the bi-fidelity DeepONet. In this layout, the validation error ε val using the proposed approach is 1.68 × 10 −2 , which is about half of 2.97 × 10 −2 , ε val obtained using the standard DeepONet. Figure 15 shows the histograms of the validation errors.
with increasing connection due to the addition of a new input node in the network. A similar but slight improvement in the proposed approach can also be noticed in Cases II and III in Table 1.

Layout 2: 64 Turbines
In this layout, we consider 64 turbines in total arranged in a 8 ⇥ 8 grid that is more characteristic of the size of a utility-scale wind plant. Other specifications of the turbine layout and the input parameters are the same as before. Figure 14 shows the two-and three-dimensional meshes used to generate the training and validation datasets. The results from the proposed bi-fidelity training of DeepONets are compared with a DeepONet (denoted as 'standard') trained only using a high-fidelity dataset with wind speed, direction, and yaw angle as branch network inputs and the turbine location as the trunk network input. Note that one 2D simulation of this layout takes 3.68 min. (average of 300 runs) compared to 1.15 hrs. (average of 300 runs) for one 3D simulation on the National Renewable Energy Laboratory's Eagle high-performance computing system. As before, for a fair comparison, we train the standard DeepONet with a training dataset D tr consisting of 106 realizations of the uncertain variables and generated power to keep the cost of generating the training dataset similar to that of the bi-fidelity DeepONet. In this layout, the validation error e val using the proposed approach is 1.68 ⇥ 10 2 , which is about half of 2.97 ⇥ 10 2 , e val obtained using the standard DeepONet. Figure 15 shows the histograms of the validation errors.

15
(a) Two-dimensional mesh with increasing connection due to the addition of a new input node in the network. A similar but slight improvement in the proposed approach can also be noticed in Cases II and III in Table 1.

Layout 2: 64 Turbines
In this layout, we consider 64 turbines in total arranged in a 8 ⇥ 8 grid that is more characteristic of the size of a utility-scale wind plant. Other specifications of the turbine layout and the input parameters are the same as before. Figure 14 shows the two-and three-dimensional meshes used to generate the training and validation datasets. The results from the proposed bi-fidelity training of DeepONets are compared with a DeepONet (denoted as 'standard') trained only using a high-fidelity dataset with wind speed, direction, and yaw angle as branch network inputs and the turbine location as the trunk network input. Note that one 2D simulation of this layout takes 3.68 min. (average of 300 runs) compared to 1.15 hrs. (average of 300 runs) for one 3D simulation on the National Renewable Energy Laboratory's Eagle high-performance computing system. As before, for a fair comparison, we train the standard DeepONet with a training dataset D tr consisting of 106 realizations of the uncertain variables and generated power to keep the cost of generating the training dataset similar to that of the bi-fidelity DeepONet. In this layout, the validation error e val using the proposed approach is 1.68 ⇥ 10 2 , which is about half of 2.97 ⇥ 10 2 , e val obtained using the standard DeepONet. Figure 15 shows the histograms of the validation errors.

Conclusions
In this paper, we propose a bi-fidelity approach for modeling uncertain and partially unknown system responses using DeepONets. In this approach, we utilize a similar and computationally less expensive model of the physical system as a low-fidelity model. The difference between the response from the low-fidelity model and the true system is subsequently modeled using DeepONets. As is often the case, with low-fidelity model we can capture the important characteristics of the system's behavior. Therefore, we make a reasonable assumption in this study that the modeling of the discrepancy between the low-fidelity model and the system's response can be done with a smaller high-fidelity dataset than a standard DeepONet to model the system's response in full to achieve similar validation error. In particular, we use DeepONets with low-fidelity model response, the locations, and the realizations of the uncertain variables as inputs. We apply the proposed approach to problems with parametric as well as with structural uncertainty, where some of the physics is unknown or not modeled. First, a nonlinear oscillator's response is modeled using this proposed approach, followed by the modeling of nonlinear heat transfer in a thin plate. Both examples show that the proposed approach provides an improvement in the validation error by almost an order of magnitude. Next, we model a wind farm with six wind turbines, where power generated from two-dimensional simulations is used as the low-fidelity response and the true system's response is obtained using three-dimensional simulations. Again, the proposed approach provides significant improvement in validation error, showcasing the efficacy of the proposed approach. These three case studies demonstrate that the bi-fidelity DeepONet approach is a powerful tool in situations with limited high-fidelity data. In the future, we plan to develop a multi-fidelity extension of this approach to determine the effect of regularization in the loss function and modeling of complex multiphysics phenomena.