Near‑boundary error reduction with an optimized weighted Dirichlet–Neumann boundary condition for stochastic PDE‑based Gaussian random field generators

Random field generation through the solution of stochastic partial differential equations is a computationally inexpensive method of introducing spatial variability into numerical analyses. This is particularly important in systems where material heterogeneity has influence over the response to certain stimuli. Whilst it is a convenient method, spurious values arise in the near boundary of the domain due to the non-exact nature of the specific boundary condition applied. This change in the correlation structure can amplify or dampen the system response in the near-boundary region depending on the chosen boundary condition, and can lead to inconsistencies in the overall behaviour of the system. In this study, a weighted Dirichlet–Neumann boundary condition is proposed as a way of controlling the resulting structure in the near-boundary region. The condition relies on a weighting parameter which scales the application to have a more dominant Dirichlet or Neumann component, giving a closer approximation to the true correlation structure of the Matérn autocorrelation function on which the formulation is based on. Two weighting coefficients are proposed and optimal values of the weighting parameter are provided. Through parametric investigation, the weighted Dirichlet–Neumann approach is shown to yield more consistent correlation structures than the common boundary conditions applied in the current literature. We also propose a relationship between the weighting parameter and the desired length-scale parameter of the field such that the optimal value can be chosen for a given problem.

Thi s v e r sio n is b ei n g m a d e a v ail a bl e in a c c o r d a n c e wit h p u blis h e r p olici e s. S e e h t t p://o r c a . cf. a c. u k/ p olici e s. h t ml fo r u s a g e p olici e s. Co py ri g h t a n d m o r al ri g h t s fo r p u blic a tio n s m a d e a v ail a bl e in ORCA a r e r e t ai n e d by t h e c o py ri g h t h ol d e r s .

Introduction
A degree of randomness is often included in numerical analyses to account for uncertainty within systems. Numerical simulation techniques rely on having a description of randomness through spatial variability across the computational domain [6,11,21,27,28]. In many cases, this can be introduced by including random fields. In its most basic form, initial conditions, material parameters, or the physical geometry can be represented by spatial values assigned by a scaled uncorrelated normal distribution having no inherent structure. However, this is often not physically meaningful, and having a correlation structure imposed on the field would result in a structure more consistent with the patterns of spatial continuity seen in physical systems, such as for example in soil bodies [19]. Existing methods of generating correlated random fields include Karhunen-Loève's expansion [12,20], local averaging subdivision methods [8], covariance matrix decomposition [5,13,23,26] and the solution of stochastic partial differential equations (PDEs) derived from Whittle-Matérn's autocorrelation function (ACF) [18,25]. The latter of these methods is computationally efficient due to sparse matrix linear algebra, and is well suited to existing FEM codes due to the construction of the PDE components and their solution. Another key aspect is the strong theoretical basis in the desired domain properties of the fields that can be produced, being due to the clear distinction between the construction of the theoretical model and the numerical methods used in the solution process [17]. The approach also allows for many generalisations of stationary Matérn fields such as non-stationary fields and those generated over less idealised domains [2,10,18]. For practical simulations, the domain must be reduced to a bounded domain of interest, requiring boundary conditions that are not generally known [14]. By applying a non-exact condition at the boundary, the approximation of the ACF that the stochastic PDE represents will breakdown, resulting in spurious values in the near-boundary region. Often, the homogenous Neumann or homogeneous Dirichlet condition is chosen due to the ease of implementation. To deal will this, the computational domain Ω is often extended, such that solving the stochastic PDE over the extended domain and extracting Ω will result in minimal effects from the applied the boundary condition. The reduction in error at the boundary is relative to the size of extension, and has been previously studied [14].
An alternative approach is to apply the Robin boundary condition through careful choice of the Robin coefficient , which can be thought of as a tuning parameter. The choice of can be problem dependent, but a choice of = 1.42l , where l is defined as the length-scale parameter, has been found to be an adequate approximation [25]. As such, the use of an extended computational domain is not strictly necessary when applying this condition. It is worth noting that this choice of = 1.42l was not established in a rigorous way, as it was deemed out of the scope of the study [25], and it was suggested that should vary as a function of the boundary. This was later considered by Daon and Stadler [7], who utilised a spatially varying Robin coefficient to provide domain Green's functions that are close to the free-space Green's functions of the Matérn covariance function. Through numerical experiments, the approach was seen to reduce the observed boundary effects, but the computation of the spatially dependent coefficient posed difficulties due to integral singularities and a required prior knowledge of pointwise variance over the domain which may not always be available. Other approaches have been taken, such as the use of a partial Dirichlet-to-Neumann operator on the extended boundary [3]. The mapping depends on the unknown correlation structure of the extended domain, and as such is an unknown itself that needs to be estimated. It was later shown to have suitable representation by a lower dimensional truncated Karhunen-Loève expansion based on information about the desired correlation structure [4]. The mapping was seen to result in reductions in computational complexity and matched well with simulated and real data.
In this study, an alternative approach is presented for reducing spurious values in the near-boundary region through a weighted Dirichlet-Neumann (D-N) boundary condition. Two variants of the condition are proposed based upon adopting different dependencies to define the Neumann coefficient used. The proposed conditions are both dependent on a weighting parameter α that controls the ratio of Dirichlet-to-Neumann components applied to the boundary, with the second also implementing a dependency based on l . Through a detailed parametric investigation, optimal values of α are found based on the length-scale parameter l . To allow for comparison between the two weighted D-N approaches, the homogeneous Neumann boundary condition, with and without domain extension, and the Robin boundary condition with = 1.42l are also applied. The error reduction at the boundary based on the applied condition is then evaluated through computing the covariance functions of the corresponding generated fields, and comparing against the true autocorrelation function. Finally, a relation between and l is proposed to generalise the condition, allowing for simpler application when considering wider problems.
The layout of the remainder of the paper is as follows; Sect. 2 presents the theory of Gaussian random field generation through the solution of stochastic PDEs; Sect. 3 introduces the weighted Dirichlet-Neumann boundary condition; Sect. 4 outlines the test cases and testing regime; Sect. 5 presents and discusses the results of the parametric investigation; and Sect. 6 presents the main conclusions of the study.

Theory and numerical discretisation
Let ∈ ℝ d be a Gaussian random field whose contents are a parameterised collection of Gaussian random variables { ( )} ∈ℝ d . The field is assumed to be stationary, such that the covariance function defining the correlated structure of the field is a function of spatial distance alone. In this way, the resulting field structure can be defined by the standard autocorrelation form, where here the Matérn autocorrelation function ACF X ( ) is chosen for ∈ ℝ d , where | | is the Euclidean distance, ν > 0 is the smoothness parameter, Γ is the gamma function, and K ν is the Bessel function of the second kind of order ν [24]. The parameter l > 0 is the length-scale parameter where δ=l √ 8ν is the distance for correlations near 0.1 [18], and controls the correlated structure of the generated fields. Equation (1) can be approximated by posing the function as a stochastic PDE, the full details of which are given in Roininen et al. [25], such that where d =1, 2,3 , is white noise on ℝ d , and α is a constant defined as The PDE in Eq. (2) is rendered elliptic by fixing the smoothness parameter as ν=2 − d∕2 , resulting in the equivalent matrix equation where is the standard identity matrix. Alternative choices of ν would lead to different -potentially fractional-differential operator powers, and would require a significant change in approach due to the methods employed to discretise the problem [9]. As such, the results that follow are associated with the value of ν implied by the regularised Laplacian operator, being a common choice in spatial statistics [2,24].
Consider a generalised random variable, i.e. a continuous linear mapping from the space of rapidly decreasing smooth functions S ℝ d to square-integrable random variables. In this way, if is also a generalised random variable, then is a square-integrable random variable for all , ∈ S ℝ d . Similarly, ⟨ , ⟩ is a Gaussian random variable by its own definition, such that Finally, by extending to be a linear function on L 2 ℝ d [16], Eq. (2) can be reduced to: find such that for all ∈ S ℝ d .
To solve the problem numerically, it must be reduced to a bounded domain. Let Ω ⊂ ℝ d be a bounded Lipshitz domain, and as such, the problem becomes: find on Ω such that Eq. (7) holds for all C ∞ 0 (Ω) . Here, the solution is non-unique, so boundary conditions need to be supplied. The well-known Dirichlet, Neumann and Robin conditions, respectively, are specified where is the unit normal to the boundary, and the Robin coefficient (a scalar value). The choice of condition imposed will result in slight changes in the resulting matrix equation to be solved.
Here, the Neumann condition Eq. (9) is considered first. To derive a weak bilinear approximation of the problem that lies in the defined problem space, we employ the finite element approximation where ψ j are the basis functions in H 1 (Ω) (Sobolev space), and X j are random variables. Applying Green's first identity to Eq. (7) results in and by making the usual Galerkin choice =ψ i , the problem can be approximated as where a is a bilinear functional defined as As such, the following matrix equation can be solved where = X j is our generated Gaussian field and and and the vector are When considering the Robin condition, following the application of Green's theorem, the bilinear functional becomes where = Ω being the domain boundary, such that the matrix equation becomes Considering the Dirichlet condition results in an equation that shares the form of Eq. (15) through choice of the function space to be H 1 0 (Ω).

Weighted Dirichlet-Neumann boundary condition
Here, two alternative coefficients 1 , 2 are proposed by formulating the Robin condition in Eq. (10) as a weighted boundary condition between the Dirichlet and Neumann components, denoted as the weighted Dirichlet-Neumann (D-N) condition. The conditions will depend on the weighting parameter α∈[0, 1] which controls the ratio of the Dirichlet and Neumann components. Let us first consider 1 , where the weighted D-N condition is Taking Eq. (19) in the form of Eq. (10), then it can be seen that Similarly, when considering λ 2 , the given boundary condition is where in this case Both conditions are functionally the same as the Robin condition, where can be used to tune the influence of its components. The inclusion of l in Eq. (22) follows the functional dependency of the proposed Robin coefficient in Roininen et al. [25], and whilst similar to Eq. (20), the additional dependence will result in a different optimal range when tuning α to reduce errors in the near-boundary region. Equations (19) and (21) can be applied by simply exchanging the expression for λ in Eq. (18) to that of 1 , 2 . To illustrate the effects of changingα , random fields were generated by solving Eq. (18) after applying Eq. (21) for a 1 m cube consisting of hexahedral elements with l = 0.2 m. Figure 1 shows a partition of the resulting fields generated over a cube for (a) = 0 , (b) = 0.5 , (c) → 1 , and (d) the result of implementing the standard Neumann condition Eq. (9) where the fields are scaled to have zero mean and a standard deviation of 1.
As → 1 , the Neumann component of Eq. (21) tends to zero, relating to a fixed (Dirichlet) condition Eq. (8) on the boundary as seen in Fig. 1c. On the other hand, when (20)  α=0 , the Dirichlet component is removed, leading to the resulting field matching the corelation structure of the pure Neumann enforcement Eq. (9) (see Fig. 1a, d). It can also be seen in (a) and (d) that variability within the domain is more diffuse than the variation over the boundary. This is commonly observed when considering Neumann boundary conditions. The converse is apparent in Fig. 1c, as → 1 , with no variability over the boundary and more pronounced structures within the domain. Figure 1b is the generated field for = 0.5 , highlighting that combing the components of both conditions can lead to correlation structures that are more consistent with those expected over the whole domain (this is more explicitly quantified in the following sections). Whilst = 0.5 may not the optimal value of for this l , there is a distinct reduction in error in the near-boundary region with the field sharing its correlation structure over the full domain.

Testing procedure
To determine the relationship between α and l for 1 and 2 , a detailed parametric testing regime was devised to find the optimal value of for a given range of l . Here, l is in relative terms, where l = 0.1 would relate to a length-scale parameter of 10% of the domain length. The fields were generated for values of l = 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, and 0.4. For each l , is varied from 0 to near 1 in increments of 0.1, with 10 realisations being generated at each to allow for the average covariance function for each to be calculated and compared with the true ACF in Eq. (1). The average is taken since it is common when modelling a stochastic problem to do so for many realisations, and find a range of responses. Once the optimal value of α is found for this level of precision, the search is continued between the two best choices of for each l . Here, the maximal value of l is taken as 0.4. When applying spatial variation using Gaussian random fields in numerical models, it is unlikely that higher values will be used due to low levels of variation across the domain defeating the purpose of imposing variability in the system. It is more likely to be at most 0.3 in relative terms. Similarly, due to the break down in the formulation, taking l larger than 0.4 will not give a reasonable approximation of the desired correlation structure for all boundary conditions considered.
To compare the generated field structures with Eq. (1), the semivariogram of each field is computed. The semivariogram is used to measure the level of spatial continuity for a given domain, being most commonly applied in the area of geostatistics [1,29]. It can be calculated as half the average squared difference of values separated by a location vector where h is the lag distance, N(h) is the number of pairs of points in Ω of separation h , is the location vector, and z( ) is the value in the domain at the vector . The covariance function and semivariogram are directly related by where σ 2 STD is the standard deviation of the field values. In this way, we can compute the covariance function of each field, average over the 10 realisations for each choice of , and compare this with Eq. (1) to determine an optimal value of .
The standard Neumann and Robin conditions, with = 1.42l , are also given the same treatment, with 10 realisations being generated for each l . In doing so, the optimal choice of boundary condition can be evaluated, and further insight to the choice of = 1.42l can be used to determine its suitability in matching the true covariance function Eq. (1). Furthermore, use of an extended domain, Ω EXT and a Neumann condition is also considered (denoted as Neumann-Extended) again with 10 realisations being generated for each l.
Here, the covariance functions are computed over the full domain Ω , as well as the near-boundary region. We define Ω B as the near-boundary region which includes all points that are a distance of less than or equal to l from Ω . As l grows, the approximation of Eq. (1) in the formulation breaks down in Ω B due to its increasing size relative to Ω . When implementing the SPDE approach, it is advisable to minimise Ω B in relation to Ω to avoid heavily constraining the problem. It is possible for complex geometry that Ω=Ω B , leading to the choice of boundary condition applied becoming irrelevant due to the propagation of error through the whole domain.

Results and discussion
In this section, results relating to the application of the following boundary conditions are presented: Neumann, Neumann-Extended, Robin with = 1.42l , weighted D-N with = 1 , and weighted D-N with = 2 . The computational domain of the Neumann-Extended case-Ω EXT -is a 1.5 m cube consisting of 97,336 nodes discretised into 91,125 8-noded hexahedral elements, the computational domain of the standard case ( Ω ) is a 1 m cube sub-domain of Ω EXT , aligned concentrically, consisting of 29,791 nodes discretised with 27,000 elements. For both domains, the relative element length is 0.03333. Once the field has been generated over Ω EXT , the 1 m cube internal domain Ω is extracted and used to compare with the other methods whose fields are generated directly over Ω . The computational domains can be seen in Fig. 2, being Ω (a) and Ω EXT (b). It is worth noting that the strongly enforced Dirichlet boundary condition is neglected here as upon its enforcement, the variation on the boundary will vanish (as in Fig. 1c), defeating the purpose of having a random field to represent variable spatial continuity. It could be argued that the domain extension principle should be applied in this case, but doing so would lead to an almost identical field as that of the Neumann-Extended condition, so it has been neglected in this study. Figure 3 shows the range of average covariance functions obtained when applying the weighted D-N boundary condition with = 1 for α from 0 to near 1 with increasing correlation length (a-d) over Ω . In (a-d), the maximum and minimum alpha will define the range over which the covariance functions lie for linearly increasing , and are compared against the ACF Eq. (1).
The covariance functions over Ω can also be visualised when considering = 2 , and is shown in Fig. 4.
In both Figs. 4 and 5, the shaded curves plotted between the maximal and minimal values relate to linear increases in between the given region.
Finally, Fig. 5 compares the average covariance function over Ω of all considered boundary conditions with the true ACF Eq. (1), where the optimal values of at each l have been chosen for both weighted D-N conditions. The coefficient of determination R 2 of the covariance functions and (1) were used to determine which resulted in the best fit, by maximising R 2 . Figure 3 depicts a clustering of curves as → 0 , suggesting that the optimal values should have a more non-linear relationship with increasing l . What is more enlightening are the results seen in Fig. 5 when comparing all tested boundary conditions. In (a), almost all applied boundary conditions result in a correlation structure that matches well with the ACF, with the exception being Neumann without extension due to much larger error in Ω B . As l increases, the performance of the proposed method progressively degrades due to Ω B becoming larger, meaning that the approximation of Eq. (1) will weaken in all cases. However even though there is a reduction in accuracy, the weighted D-N with optimal α appears to give better matches to the ACF than the standard applied boundary conditions. This can be quantified by computing the R 2 values of the curves presented in Fig. 5 to evaluate their performance in matching the correlation structure of Eq. (1) with increasing l . Consequently, the optimal values for α at each l used in Fig. 5 were determined through maximising R 2 when measuring the goodness of fit between the weighted D-N covariance function curves over Ω with the ACF Eq. (1). Figure 6 shows the full comparison of all considered boundary conditions in terms of their R 2 value for average covariance functions computed over Ω and Ω B , where the suffix "-F" and "-B" denote calculation over the full domain and near boundary, respectively.
It can be seen in Fig. 6 that the weighted D-N approach yields a more consistent match to the ACF. Almost all applied boundary conditions result in fields that follow a similar pattern of a decreased wellness of fit as l increases, and similarly after a certain point as l → 0 . An exception is seen when the Neumann boundary condition is applied. As l → 0.5 , the correlation structure becomes unchanging, where further increasing l will have no effect on the static structure. By this, we can assume that the correlation structure is beginning to stabilise at l = 0.4 , which can be seen in the difference in the covariance function of the Neumann boundary condition in Fig. 5c, d. Thus, as l increases, the ACF will shift closer to the static Neumann covariance curve, resulting in a larger R 2 value but not necessarily a truer correlation structure. It is also worth noting that when Fig. 2 The computational domains used for field generation, being: a Ω , a 1 m cube with 29,791 nodes, and b Ω EXT , a 1.5 m cube with 97,336 nodes l = 0.5 , our domain Ω=Ω B , suggesting that the error in the correlation structure at the boundary can be seen throughout the whole domain regardless of the applied boundary condition. Figure 6 also shows that the Robin condition with = 1.42l may not always be the best choice for approximating the ACF due to its sharper fall in R 2 . This supports the suggestion in Roininen et al. [25] that could be given as a function on the boundary.
The effect of the applied boundary condition on the level of error in Ω B can also be quantified in this way, being visualised in Fig. 6 as the separation between a given boundary conditions R 2 curves when calculated over Ω and Ω B respectively. The largest level of disparity is seen when the Neumann boundary condition is applied without extension, where the covariance functions calculated over Ω B appear to have a better fit. This can be further seen in Fig. 1d, where the centre of the domain appears more diffuse. On the other hand, Neumann-Extended results in the least error over Ω B . This is due to the computational domain being extended, and the errors seen in the near boundary of will not be carried through to the inner cube Ω . The weighted D-N approach also has marginal differences in R 2 over Ω and Ω B , where the choice of = 1 or = 2 has negligible effects on the level of difference. As l decreases, the near-boundary effects can be seen to reduce to the point where the curves appear identical, suggesting that an appropriate choice of can mitigate error seen over Ω B . This is in contrast to the Robin condition, whose near-boundary effects begin to increase as l decreases from 0.15 to 0.1, and the Neuman-Extended condition whose nearboundary effects increase over the full domain but decrease over the near-boundary region with decreasing l in this range (see Fig. 7). When generating fields with Neumann-Extended boundary conditions, the combination of domain extraction and a decreasing l suggests the error in the near boundary and full domain will converge due to the shrinking of near-boundary region as l reduced. After a certain point, the difference between the boundary regions of the weighted D-N conditions begins to grow as seen in Fig. 7, and could be mitigated with slight variation in the chosen α or using a finer mesh.
Finally, the relationship between the optimal values of α with respect to l can be determined as seen in Fig. 8.
The functions can be described as (25) (l) = al 2 + bl + c This suggests that upon implementing the weighted D-N condition, small variations in the chosen α from its optimal value with have a less detrimental effect on the resulting correlation structure when = 2 as opposed to = 1 . As l increases, the optimal value of α reduces in both cases. This relates to a more dominant Neumann component in the condition. This also agrees with the increase in R 2 for the Neumann condition, suggesting that at larger l , a larger Neumann component would result in a better matching correlation structure. It is also worth highlighting that both fitted curves converge to a common root l = 0.445 . Here, the Dirichlet component of the weighted D-N condition vanishes, suggesting that the pure Neumann condition will give just as sufficient an approximation.

Mesh convergence
A similar testing procedure was carried out as above for l = 0.1 over the same 1 m cube as in Fig. 2a for different mesh sizes to determine if the solution is mesh converged. Here, the weighted D-N boundary condition was applied with = 2 . The mesh was divided into relative element lengths L e of 0.1, 0.05, 0.03333, and 0.025 for the regular hexahedral elements. The testing range of was chosen as 0.41-0.49 to assess the capabilities of the relationship between and l given in Eq. (25) and shown in Fig. 8, where the functions exact value gives as 0.45 . Here, the goodness of fit was determined by Root Mean Squared Error (RMSE) and R 2 values of the covariance function plots to the true ACF. The RMSE was used here as an additional metric as it quantifies the error as absolute values as opposed to the R 2 which is a percentage based indicator and can be less interpretable for ill-fitting curves. The covariance functions compared with the ACF for different mesh sizes can be seen in Fig. 9.
Apart from the element length 0.1 case (Fig. 9a), the covariance plots match well with the true function, showing that for all mesh sizes, the mesh is converged and the 2 relationship Eq. (25) gives a sufficiently accurate choice of . This can be seen further in Table 1, where the R 2 and RMSE values are presented for the considered mesh sizes when = 0.45 , with the RMSE values being visualised in Fig. 10. The convergence of the mesh-as well as the confirmation of the 2 relationships in Eq. (25) applicability-suggests that the values of fit through the parametric investigation are reasonable, and are done so on an appropriate mesh.

Dog bone example
To further test the applicability of the weighting parameter, random fields were generated over a dog bone shaped specimen, whose shape is synonymous with the experimental determination of tensile properties of cement composites [15,22,30]. The domain Ω was discretised to contain 91,203 nodes, where the mesh and relative dimensions can be seen in Fig. 11a.
The length scale given is relative to the central portion of the dog bone, whose dimensions are of 1 in both the x-and z-axis. The size and shape of the domain was chosen such that the percentage of the domain which is dominated by Ω B is relatively low. This ensured that the problem formulation is not heavily constrained due to the domain shape, and the error seen in the near-boundary region does not propagate through the whole domain. Similar to the parametric regime, 10 fields were generated with a length scale of l = 0.25relative to the dimensions of the central portion-for each applied boundary condition, with = 0.32 for the weighted D-N condition with = 2 . A sample field is shown in Fig. 11b when using the weighted D-N condition with = 0.32 . The average covariance function for each applied boundary condition was determined, as seen in Fig. 12.
As expected, the Neumann case offers the worst fit to the function, with the Robin and weighted D-N conditions offering a marginal under and over estimation respectively. The choice of adopted here is based on a function obtained from a domain where the relative length scale is well defined in all directions. For the case of the dog bone, the true relative length scale differs in each portion of the domain, where the top and tail sections have a smaller associated l than that considered. This will be reflected in the chosen , where a smaller l suggests taking a larger . In this way, taking a larger alpha will shift the covariance function closer to the true ACF, converging to the best possible fit in terms of the error estimation for this domain and applied boundary condition. Both the robin and weighted D-N provide an R 2 > 0.99 from Fig. 12, suggesting that both are more than sufficient to apply in this case. The benefit of applying the weighted D-N condition here is the choice of whether or not to tune it by small variations in α to further reduce the error.
The near-boundary error associated with complex geometries is related to the size of the near-boundary region relative to the full domain, as well as the distinct portions that it may encompass. Having complex geometry could lead to regions of the domain that are entirely made up by the near-boundary region due to its size being dependent on the length scale. Having a smaller relative length-scale value would suggest a reduction of the near-boundary region, so the chance of having zones that are purely made up of the near-boundary region would decrease. However, it is possible that small relative length-scales will not be sufficient to reduce this as there could be thin sections of domain that contain points that are always close enough to the boundary to be considered inside the near-boundary region. In extreme cases, where much of the domain is considered as the near boundary, a convenient but less computationally efficient choice may be to apply the Neumann-Extended boundary condition to an extended domain which envelopes the desired numerical domain, such that it is largely out of the near-boundary region and can be extracted for further use.

Concluding remarks
An optimised boundary condition has been proposed and applied to the formulation of the stochastic PDEs solved to generate correlated Gaussian random fields. The condition aimed at reducing spurious values in the near-boundary region, and was found to perform well inside and outside of the range of the parametric analysis that followed. It was found that the near-boundary error of the weighted D-N approach reduced with decreasing l when employing optimal   values of , and is more consistent with the correlation structure of the full domain. Similarly, the weighted D-N approach provides an overall better match to the autocorrelation function when compared with all other applied boundary conditions. This was also tested for complex geometry and was found to perform well. The functions for based on the parametric study enable its optimal value to be determined for other domains and desired length-scales. As such, the weighted D-N approach is recommended as being the applied boundary condition when formulating the stochastic PDEs to be solved. Whilst both = 1 and = 2 perform equally well, = 2 is a more practical choice. When = 2 , the resulting covariance functions show that changing leads to more uniform variation in resulting covariance structure, providing a more linear relationship than if = 1 . This suggests that choosing an α with small variations from its fitted value will have less of an impact on the final correlation structure than if = 1 , making it more consistent when applying the condition outside of the tested range or on less regular domains. The other main advantage, as opposed to the Neumann boundary condition with extension, is the lack of dependence on an extended domain. In this case, the computational domain was chosen as a cube with relatively low numbers of elements, so computational expense did not need to be considered when solving over an extended domain. If the field generation method was conducted over a finer mesh, then this dependence could cause complications. Finally, the weighted D-N approach was shown to be significantly more accurate for larger values of l as opposed to the Robin condition with = 1.42l , thus providing better control over the resulting correlation structure where the formulation begins to break down.
The relationships presented are given in relative terms, being highly applicability to a wider range of domains and engineering problems. The proposed boundary condition is aimed at giving a more sufficient approximation of the correlation structure of said parameters, allowing for a more robust quantification of uncertainty through numerical analyses.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, Fig. 11 a Dog bone unstructured mesh consisting of hexahedral elements, b random field realisation when using the weighted D-N boundary conditions with = 0.32 for l = 0.25 (scaled to zero mean and a standard deviation of 1)

Fig. 12
Average correlation function comparison across the dog bone domain when applying the Robin, Neumann, and Weighted D-N boundary condition with = 0.32 provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.