A domain mapping approach for elliptic equations posed on random bulk and surface domains

In this article, we analyse the domain mapping method approach to approximate statistical moments of solutions to linear elliptic partial differential equations posed over random geometries including smooth surfaces and bulk-surface systems. In particular, we present the necessary geometric analysis required by the domain mapping method to reformulate elliptic equations on random surfaces onto a fixed deterministic surface using a prescribed stochastic parametrisation of the random domain. An abstract analysis of a finite element discretisation coupled with a Monte-Carlo sampling is presented for the resulting elliptic equations with random coefficients posed over the fixed curved reference domain and optimal error estimates are derived. The results from the abstract framework are applied to a model elliptic problem on a random surface and a coupled elliptic bulk-surface system and the theoretical convergence rates are confirmed by numerical experiments.


Introduction
In the mathematical characterization of numerous scientific and engineering systems, the topology of the domain may not be precisely described.The main sources of uncertainty are usually insufficient data, measurement errors or manufacturing variability.This uncertainty in the geometry often naturally appears in many applications including surface imaging, manufacturing of nano-devices, material science and biological systems.As a result, the analysis of uncertainty in the computational domain has become an interesting and rich mathematical field.A comprehensive summary concerning the first directions in the treatment of elliptic partial differential equations (PDEs) in random domains can be found in [4,18,26,28,32] and recently [10] for a parabolic equation on a randomly evolving domain.Aside from the fictitious domain method [4,27,28], the main approaches utilize a probabilistic framework by describing the random boundary of the domain with a random field.This probabiltistic approach is usually proceeded with one of two main techniques: the perturbation approach and the domain mapping method.The perturbation approach [17,20] exploits a shape Taylor expansion with respect to the boundary random field to represent the solution, however as a result is limited to consideration of only small random deformations.The domain mapping approach [5,19,32] on the other hand allows does not suffer the same limitations.The key idea behind this method is to define an extension of the random boundary process into the interior domain to form a complete random mapping for the whole domain and then to use this domain mapping to transform the original partial differential equation on the random domain onto the fixed deterministic reference domain resulting in an partial differential equations with random coefficients.For the latter formulation, there is a wealth of literature available on numerical techniques to compute any quantities of interest, see for example [16,23,24].The aim of this paper, is to incorporate the domain mapping method with the well-developed field of surface PDEs [8,13,14] which has so far only considered uncertainty in the coefficients of the considered PDEs, see [11].This will lead to more realistic geometric description of many of the situations previously dicussed.Note that while the domain mapping method will be applicable to domains with random rough surfaces, we will only choose to focus on sufficiently smooth random surfaces and leave the rough case for future considerations.The layout of the article is as follows.In Section 2, we provide an overview of the domain mapping method for partial differential equations in flat random domains and furthermore discuss suitable notions for the expectation of a family of random domains.In Section 3, we introduce the necessary geometric analysis and computations required to apply the domain mapping method to elliptic partial differential equations posed on random surfaces.In Section 4, we present a model elliptic problem on a random surface and a coupled elliptic system on a random bulk-surface, and analyse weak formulations in both the stochastic and spatial variables for the reformulated equations with stochastic coefficients on the fixed reference surface and bulksurface respectively.Section 5 provides an abstract analysis of a finite element discretisation incorporating a pertubation to the variational set-up due to a first order approximation of the curved reference domain, and couples with a Monte-Carlo sampling to approximate the first moment of the solution.An optimal error estimate is derived and subsequently applied in Section 6, to two discretisations of the proposed reformulated problems.We conclude in Section 7, by presenting numerical results confirming the theoretical rate of convergence.

The domain mapping method
We begin with a brief introduction on spaces of random fields.For further details on these spaces, we refer the reader to [24].Note throughout this paper, we will let (Ω, F, P) denote a complete, separable probability space consisting of a sample space Ω, a σ−algebra of events F and a probability measure P.
2.1.Random field notation.For a given Banach space V and p ∈ [1, ∞], the Lebesgue-Bochner space L p (Ω; V ) consists of all strongly F−measurable functions f : Ω → V for which the norm is finite.For convenience, we will express the parameters of a given random field (f (ω))(x) by f (ω, x).In the case that V is a separable Hilbert space, it follows that L 2 (Ω; V ) is also a separable Hilbert space and furthermore is isomorphic to the tensor product For details, see [29].
2.2.The domain mapping method.To illustrate the key concepts of the domain mapping method, consider the following boundary value problem posed on an open, connected, bounded domain D(ω) ⊂ R 2 with a random boundary Γ(ω) = ∂D(ω).Here the prescribed random field f (ω) : D(ω) → R and additionally the boundary, will be assumed to be sufficiently regular to ensure well-posedness for a.e.ω.The first essential feature of the domain mapping method is the representation of the stochastic boundary via a random field.More precisely, in the above context we will assume that there exists a random field φ ∈ L ∞ (Ω; C 0 (Γ 0 ; R 2 )), that maps a fixed closed curve Γ 0 ⊂ R 2 onto realisations of the random boundary φ(ω, •) : Γ 0 → Γ(ω), see figure 1.The next step in the method is to define an extension of the boundary process into the interior to form a stochastic mapping φ(ω, •) : D 0 → D(ω) for the whole domain.For instance, [32] proposed an extension based on the solution of the Laplace equation over the unit square with boundary conditions prescribed by segments of the random boundary.However, alternative approaches may wished to be considered depending on the application in question and the geometry of the computational reference domain.With a complete domain mapping at hand, the random domain problem (2.2) can now be reformulated as a partial differential equation with random coefficients over the fixed deterministic domain D 0 , where the specific random coefficients for this particular problem are given by G(ω) = ∇φ (ω)∇φ(ω) g(ω) = det G(ω).
We now have access to a wide breadth of numerical techniques, including Monte-Carlo [6,21] and the stochastic Galerkin method [2,25], to compute any statistical quantities of interest.
Remark 2.1.Note that the choice of the reference domain D 0 , for the stochastic domain mapping φ describing the complete random geometry in question, is arbitrary and should be chosen in such a way that simplifies the computation at hand.Furthermore in practice, only statistical properties such as the expectation and two-point covariance function of the stochastic mapping φ will be known.As a result, an approximation of the true process will instead be used, commonly taking the form of a truncated series with centered, uncorrelated random coefficients Y k with unit variance, such as a truncated Karhunen-Loève expansion.Considerations of the induced error is beyond the scope of this paper and we instead refer the reader to [18].

2.3.
Expected domain and quantity of interest.In order to give a precise definition of our quantity of interest, which for our purpose shall be some notion of a mean solution, we will first need to fix a suitable domain of definition.A natural choice would be the parametrisation based expected domain, introduced in [7] for random star-shaped domains, which we shall generalise as follows.
Remark 2.2.Note that there are other alternative methods in which to define the expected value of a family of random sets.For example, we could characterise the random set D(ω) as an indicator function 1 D(ω) and then use its average, the so-called coverage function p(x) = P(x ∈ D(ω)) to define the expected value to be set ≥ λ}, where the parameter λ > 0 is selected in a such a way that the volume of E V [D] is close as possible to the expected volume of the random sets D(ω).This is known as the Vorob'ev expectation and was shown in [7] not to coincide with the parameterisation based expectation.Although there is no canonical definition of the expected value of a random domain, the parametrisation based expected domain fits naturally in the setting of the domain mapping method and thus will be adopted.Assumption 2.1.We will assume that the expected value of the stochastic mapping is bi-Lipschitz continuous.This will ensure that the parametrisation based expected domain E[D] is also Lipschitz continuous and furthermore of the same dimension as D 0 and D(ω).
We will denote the induced zero-mean stochastic mapping between the parametrisation based expected domain E[D] and realisations of the random domain D(ω) by (2.5) See figure 2 for an illustration of the different mappings and domains.Our quantity of interest can now be defined on the expected domain as follows.
Definition 2.2 (QoI).Given a random field u(ω, •) : D(ω) → R defined over the family of random Lipschitz domains given in (2.3), the expected value of the random field is given by As previously discussed our aim is to apply the domain mapping method for random domains which involve random surfaces.We will therefore now proceed with some preliminary computations of geometric quantities as well as tangential derivatives of functions given over parametrised hypersurfaces in terms of quantites of the reference surface and derivatives of the domain mapping and corresponding pull-back function.This will provide a basis for the domain mapping method to be employed to several model PDEs over random surfaces.

Computations for the pull-back of tangential differential operators and geometric quanitities of parametrised hypersurfaces
Let us first introduce some notation for hypersurfaces that will be adopted throughout this paper.For a more detailed introduction, see [14].
3.1.Hypersurface notation.A set Γ ⊂ R n+1 is said to be a C k -hypersurface for k ∈ N ∪ {∞}, provided that for every x ∈ Γ there exists an open set U ⊂ R n+1 containing x and a smooth function ϕ ∈ C k (U ) such that ∇ϕ(x) = 0 on U ∩ Γ and The unit normal vector field ν Γ to the hypersurface Γ can be computed via with a choice of orientation.For a differentiable function f : Γ → R, we define the tangential gradient by is the projection operator mapping onto the tangent space T Γ to the hypersurface Γ and f is a smooth extension of f to an open neighbourhood in R n+1 .It can be shown that the tangential gradient is independent of the extension chosen [14, Lemma 2.4] and we shall denote its components by ) .We shall further denote the tangential derivative of the unit normal by H Γ = ∇ Γ ν Γ and will refer to this matrix as the extended Weingarten map.It can be shown that H Γ is symmetric with a zero eigenvalue corresponding to the unit normal vector ν Γ and furthermore agrees with the Weingarten map when restricted to the tangent space T Γ, see [8] for details.The Laplace-Beltrami operator is then defined for twice differentiable functions as follows We next introduce the Fermi coordinates with the following well-known lemma [14,Lemma 2.8].These are a global coordinate system defined in an open neighbourhood around Γ in which every point can be uniquely expressed in terms of its signed distance d Γ (x) and its closest point a Γ (x) on the surface Γ .
Lemma 3.1.Let d Γ denote the signed distance function to Γ oriented in the chosen direction of the unit normal vector field ν Γ .Then there exists δ > 0 such that for every x ∈ U δ := {y ∈ R n+1 | |d Γ (y)| < δ} there exists a unique point a Γ (x) ∈ Γ that satisfies 3.2.Geometric settings.As a point of reference, we will now describe the deterministic geometric settings that will be considered for the parametrised surfaces in the subsequent calculations.In each case, the reference surface Γ 0 ⊂ R n+1 will be assumed to be of class at least C 2 and oriented by the unit normal vector field ν Γ0 .The general geometric setting for the parametrised surface will be as follows.
We will further consider the special case, where the parametrised surface has the following graph-like representation over the reference surface.
Geometric setting 2 (Graph-like surface).The surface Γ ⊂ R n+1 will be prescribed by for a given height function h : Γ 0 → R defined over the reference surface.
Additionally, we will consider the case where the surface is compact (and thus without a boundary) and encloses an open bulk domain.

Geometric setting 3 (Parametrised bulk-surface).
The open bulk domain D ⊂ R n+1 and its boundary Γ = ∂D which is a surface, will be given by Note that in each case, the given parametrisation φ will be assumed to be a sufficiently smooth diffeomorphism for the calculation in question.Furthermore, we shall denote the associated pull-back of a given function f defined over the parametrised domain onto the reference domain by

3.3.
The tangential gradient and Laplace-Beltrami operator.Considering a general parametrised hypersurface Γ as described in (3.8), we will now compute expressions for the pull-back of the tangential gradient ∇ Γ and Laplace-Beltrami operator ∆ Γ onto the reference surface Γ 0 under the domain mapping φ.
As a motivation for these calculations, let us first recall that for a given local parametrisation of the hypersurface Γ, where U ⊂ R n and W ⊂ R n+1 denote open sets, we can express the tangential gradient and Laplace-Beltrami operator in local coordinates as follows where F = f • X and the first fundamental form G : U → R n×n is defined as G = ∇X ∇X with g = detG.
In deriving expressions for the pull-back onto the reference surface Γ 0 instead of the local coordinates, we will see similar expresssions to (3.12) and (3.13) but with the first fundamental form replaced by the following tensor G Γ0 : Γ 0 → R (n+1)×(n+1) defined by where we will similarly denote its determinant by g Γ0 = detG Γ0 .This tensor can be seen to arise by considering a local parametrisation σ : U → V ∩ Γ 0 of the reference surface Γ 0 , with V ⊂ R n+1 denoting an open set, and the induced local parametrisation of the hypersurface Γ.By computing the first fundamental form with the chain rule, we observe that Since ∇ Γ0 φ ∇ Γ0 φν Γ0 = 0 and its restriction to the tangent space maps ∇ Γ0 φ ∇ Γ0 φ : T Γ 0 → T Γ 0 , we are able to extend in the normal direction as in (3.14) to form an invertible matrix.Furthermore as ∇σ ∈ T Γ 0 , it follows that we have Note that the given extension (3.14) in the normal direction is a natural choice since the surface measures dA Γ and dA Γ0 of the respective surfaces can be shown to satisfy the relation dA Γ = √ g Γ0 dA Γ0 under the domain transformation mapping φ.We now continue by proving that a similar expression to (3.12) holds for the pull-back of the tangential gradient.
Lemma 3.2 (Tangential gradient).Given any differentiable function f : Γ → R, the pull-back of the tangential gradient onto the reference surface Γ 0 is given by Proof.Differentiating the associated pull-back function f = f • φ and applying the chain rule for tangential derivatives gives (3.16) Since the tangential gradient of the surface parametrisation bijectively maps ∇ Γ0 φ : Γ and additionally has kernel equal to span{ν Γ0 }, we see that in order to invert the matrix ∇ Γ0 φ, we must first modify the corresponding linear map to bijectively map the space span{ν Γ0 } into span{ν Γ • φ}.One possible solution is to add the linear map L : R n+1 → R n+1 characterised by which translates to adding the following tensor product (3.17) and thus leads to (3.15).For the second equality, we again use the property that the restriction ∇ Γ0 φ : Hence we deduce α = G −1 Γ0 ∇ Γ0 f and obtain the second equality.Remark 3.1.Note that the chain rule for tangential gradients (3.15) holds for any choice of orientation of the unit normals ν Γ0 and ν Γ as a result of (3.17).
Let us denote the given extension of the tangential gradient of the surface parametrisation appearing in (3.15) and furthermore denote its determinant by b = det B and the entries of its inverse by B −1 = b ij i,j .We observe with the orthogonality result ∇ Γ0 φ (ν Γ • φ) = 0 which follows from the property that the restriction maps ∇ Γ0 φ : Consequently, we have We can now compute the pull-back of the Laplace-Beltrami operator onto the reference surface as follows.
Lemma 3.3 (Laplace-Beltrami operator).Given any f : Γ → R twice differentiable, the pull-back of the Laplace-Beltrami operator is given by Proof.By the chain rule for tangential gradients (3.15), we can express the Laplace-Beltrami operator as Writing in divergence form gives The last step follows from the observations (3.19) and (3.20).We continue by proving that the remaining terms vanish.Recalling Jacobi's formula for the derivative of a determinant D Γ0 j detB = detB trace B −1 D Γ0 j B and computing the derivative of the inverse matrix It therefore follows after relabelling indices that By the symmetry of the Weingarten map D Γ0 l ν Γ0 j = D Γ0 j ν Γ0 l , we see that the last two terms cancel.We next interchange tangential derivatives [14, Lemma 2.6] (3.22), we arrive at the following expression for the remaining terms Examining the first term, we have , the first term vanishes.For the second and third term, we observe that which can be seen by We next compute the specific form of the coefficients appearing in the pull-back of the tangential gradient (3.15) and the Laplace-Beltrami operator (3.21), for the particular case of a graph-like parametrisation over the reference surface.Lemma 3.4 (Graphical case).Assuming that the parametrisation of the hypersurface Γ has the particular graph-like representation described in (3.9) for a given height function h : Γ 0 → R, then the inverse and determinant of the tensor G Γ0 defined in (3.14) simplify to give

24)
Here A := I + hH Γ0 −1 and {κ Γ0 j } j denotes the non-zero eigenvalues of the extended Weingarten map H Γ0 .Proof.Differentiating the given surface parametrisation φ(x) = x + h(x)ν Γ0 (x), we obtain Expanding the tensor G Γ0 = ∇ Γ0 φ ∇ Γ0 φ + ν Γ0 ⊗ ν Γ0 and cancelling the orthogonal terms with the tensor product identity (a Taking the inverse with the identity (I + a ⊗ b) −1 = I − a⊗b 1+a•b we obtain (3.23).For (3.24), we take the determinant and apply det( and thus obtain the stated result for √ g Γ0 = detG Γ0 . 3.4.The unit normal and extended Weingarten map.We continue by computing expressions for the pull-back onto the reference surface Γ 0 , of the unit normal ν Γ and extended Weingarten map H Γ for a general parametrised hypersurface Γ as given in (3.8).To obtain an expression for the unit normal, we smoothly extend the given surface parametrisation φ : Γ 0 → Γ to a C 1 −diffeomorphic mapping φ : U → V between some open sets U and V containing Γ 0 and Γ respectively.The existence of such a mapping is gauranteed by the Whitney extension theorem [31].We now have a level-set description of Γ consequently leading to the following expression for the unit normal vector field due to (3.2).
Lemma 3.5 (Unit normal).The pull-back of the unit normal vector field ν Γ of the parametrised surface Γ described in (3.8) onto the reference surface Γ 0 , is given by Note that (3.25) can be shown to be independent of the extension chosen.As an example of a possible extension of the given surface parametrisation, we now consider the case of a graph-like surface.
Corollary 3.1 (Graphical case).Assuming that the hypersurface Γ has the particular graph-like form described in (3.9), then the pull-back of the unit normal vector field ν Γ is given by Here the orientation has been chosen to coincide with the reference surface Γ 0 when the height function is identically zero.Recall that A := (I + hH Γ0 ) −1 .
Proof.We extend the given surface parametrisation φ : Γ 0 → Γ defined by For δ > 0 sufficiently small, its image V = φ(U ) is contained within the neighbourhood in which the Fermi coordinates (a Γ0 (x), d Γ0 (x)) are well defined.Consequently, the extension φ : U → V which equivalently acts upon the Fermi coordinates as follows can be seen to be a bijective mapping.Computing its derivative and evaluating on the reference surface Γ 0 , recalling that ∇d Γ0 = ν Γ0 and ∇a Γ0 = P Γ0 on Γ 0 by (3.6) and (3.7), we obtain Hence taking the inverse with the identity and thus obtain the stated result.Note that ν Γ0 − A∇ Γ0 h = 0 since the matrix We next compute the pull-back of the extended Weingarten map H Γ for a general parametrised surface Γ.Since the restriction of the derivative of the surface parametrisation maps ∇ Γ0 φ(•) : and rewrite (3.27) as It therefore follows from an application of the chain rule and the symmetry of the extended Weingarten map that We can then extend the tangential derivative ∇ Γ0 φ to an invertible matrix as previously discussed in (6.17), to obtain the following result.Lemma 3.6 (Extended Weingarten map).Let the orientation of the parametrised hypersurface Γ described in (3.8), be fixed by a choice of a unit normal vector field ν Γ .Then the pull-back of the extended Weingarten map is given by Note that the matrix L(x) given in (3.28) is symmetric even though the tangential derivatives do not necessarily commute, as by interchanging the derivatives we obtain and since D Γ0 m φ • ν Γ • φ = 0, for all m = 1, ..., n + 1, we see that the last two terms vanish.3.5.The normal derivative at the boundary.We conclude this section by computing the pull-back of the normal derivative at the boundary for functions defined over the parametrised bulk-surface described in (3.10).Lemma 3.7 (Normal derivative).Given any u : D → R sufficiently smooth, the pull-back of its normal derivative is given by where G = ∇φ ∇φ and g = det(G) denoting its determinant.
Proof.Differentiating u = û • φ −1 and substituting in the expression (3.25) for the pull-back of the unit normal ν Γ , where the orientation has been chosen to be in the outer direction to the domain D gives We next observe with the decomposition ∇φ = ∇ Γ0 φ + ∂φ ∂νΓ 0 ⊗ ν Γ0 and the orthogonality result Since φ maps the boundary Γ 0 onto Γ, it follows that ∂φ ∂νΓ 0 We now continue by showing that the normal component of ∂φ ∂νΓ 0 can be expressed as the ratio between the bulk √ g and the surface area element √ g Γ0 .This will be achieved in the context of exterior algebras.
Let τ 1 , ..., τ n represent an orthonormal basis of the tangent space T Γ 0 and thus {τ 1 , ..., τ n , ν Γ0 } forms a basis of R n+1 .The determinant of linear map corresponding to ∇φ evaluated on the boundary Γ 0 can be expressed in the notation of exterior algebras as follows Since ∇ Γ0 φτ 1 , ..., ∇ Γ0 φτ n form a basis of the tangent space T Γ and the exterior product of any set of linearly dependent vectors is zero, we are therefore able to remove the tangent component of the normal derivative yielding Observing that each term in the above exterior product is the image of the basis {τ 1 , ..., τ n , ν Γ0 } under the linear mapping We thus obtain the stated result with the following observations 4. First applications of the domain mapping method to complex random geometries involving random surfaces We will now consider two model elliptic problems posed on complex random domains involving random surfaces.In particular, the first problem will be posed on a sufficiently smooth random surface and the second on a random bulk-surface.In both cases, the complete random domain mapping will be assumed to be known.Furthermore, we will assume that the computational domain was chosen to coincide with the expected domain, and thus will assume in both cases that E[φ] = 0. We will now employ the domain mapping method, and reformulate both equations onto their corresponding expected domain and prove well-posedness as well as a regularity result.

4.
1.An elliptic equation on a random surface.Let Γ(ω) represent a random, compact We will assume that the random domain mapping φ(ω, •) : Γ 0 → Γ(ω) is a C 2 − diffeomorphism for almost every ω and furthermore satisfies the uniform bounds for some constant C > 0 independent of ω.We consider the following model elliptic equation on the random surface for a given random field f (ω, •) : Γ(ω) → R. Our goal is to analyse the mean solution defined by Reformulating (4.3) onto the expected domain with the calculation of the Laplace-Beltrami operator provided in Lemma 3.2 yields where the random coefficient is given by with g Γ0 (ω) = det G Γ0 (ω).Multiplying through by surface area element g Γ0 (ω) and integrating by parts, we arrive at the following mean-weak formulation on the fixed deterministic domain Γ 0 .
Proposition 4.1.Under the uniformity assumptions (4.2) on the random domain mapping, there exists constants C DΓ 0 , C gΓ 0 > 0 such that the singular values σ i of D Γ0 and the surface area element √ g Γ0 are bounded above and below by for all x ∈ Γ 0 and a.e.ω.
Proof.We can rewrite G Γ0 using the orthogonality ∇ Γ0 φ (ν Γ • φ) = 0, as follows Examining each term separately, we see that the inverse is given by Hence it follows Therefore with (4.2), we have uniform bounds above and below on the singular values of G Γ0 (ω) and hence obtain the estimates (4.10) and (4.11).
A direct consequence of the above uniform bounds on the random coefficients is the existence and uniqueness of a solution to (4.6) gauranteed by the Lax-Milgram theorem.
By considering the original surface equation (4.3) on Γ(ω) ∈ C 2 , we would expect from standard elliptic surface regularity results that for given f (ω) ∈ L 2 (Γ(ω)), the pathwise solution belongs to u(ω) ∈ H 2 (Γ(ω)) and therefore û(ω) ∈ H 2 (Γ 0 ) for a.e.ω.However since the H 2 a-priori estimate on u(ω) will naturally depend on the geometry of the realisation Γ(ω), it is not immediately clear whether the solution to the mean-weak formulation belongs to û ∈ L 2 (Ω; H 2 (Γ 0 )).We will therefore continue by explicitly treating all arising constants and their dependency on the geometry of the random domain.
Theorem 4.2 (Regularity).Given any f ∈ L 2 (Ω; L 2 (Γ 0 )), the solution to (4.6) belongs to û ∈ L 2 (Ω; H 2 (Γ 0 )) and furthermore satisfies the following estimate Proof.Let us consider the push-forward u = û • φ −1 of realisations of the weak solution onto Γ(ω) for almost every ω, which as a result of the tensor structure L 2 (Ω; ) and therefore û(ω) ∈ H 2 (Γ 0 ).For the a-priori estimate (4.13), it was shown in [14] through a series of integration by parts and interchanging of tangential derivatives that the H 2 semi-norm satisfies . Here H Γ(ω) = trace H Γ(ω) is the mean-curvature.Hence with the uniform bounds (4.2) on the random domain mapping and the previously calculated expression (3.29) for the Weingarten map, we obtain an upper bound on the constant c(ω) independent of ω.Thus, with the PDE (4.14) pointwise we have the bound We can now pull-back onto the expected domain, applying the norm equivalence of the pull-back transformation for k = 0, 1, 2 and a.e.ω where the constants are independent of ω due to bounds (4.2), and the stability estimate (4.12) to obtain and thus the stated result.4.2.A coupled elliptic system on a random bulk-surface.For the second problem, we consider a coupled elliptic system on a random bulk-surface motivated by the deterministic case analysed in [12].More precisely, the geometric setting is as follows.We let {Γ(ω)} denote a family of random, compact C 2 −hypersurfaces in R n+1 enclosing open domains D(ω) and will denote the outer unit normal by ν Γ(ω) .The family of random domains will be prescribed by the mapping (4.15) φ : where the reference surface Γ 0 ⊂ R n+1 will also be a compact C 2 −hypersurface with open interior D 0 .We will assume that the domain mapping is a C 2 −diffeomorphism for a.e.ω ∈ Ω and additionally satisfies for a constant C > 0 independent of ω.The proposed coupled elliptic system on the random bulk-surface reads as follows Here α, β > 0 are given positive constants and f (ω, •) : D(ω) → R and f Γ (ω, •) : Γ(ω) → R are prescribed random fields.As with our previous problem, our quantity of interest is the mean solution, that is the pair Let us continue by reformulating the system (4.17)onto the expected domain D 0 with our previously calculated expressions for the Laplace-Beltrami operator (3.21) and the normal derivative (3.30) giving Here the random coefficients are with g(ω) = detG(ω), g Γ0 (ω) = detG Γ0 (ω).For convenience, we have set fΓ0 = f Γ • φ.To derive a meanweak formulation, we follow the variational approach presented in [12].We begin by multiplying through the bulk equation (4.18a) by the area element √ g and integrating by parts which gives Similarly, for the surface equation (4.18c) we integrate by parts recalling that the hypersurface Γ 0 is without boundary, to obtain Taking the weighted sum and substituting in the reformulated Robin boundary condition (4.18b), we arrive at the following mean-weak formulation: for every φ ∈ L 2 (Ω; H 1 (D 0 )) and ξ ∈ L 2 (Ω; H 1 (Γ 0 )).Here we set D Γ0 (ω) = g Γ0 (ω)G −1 Γ0 (ω).We denote the associated bilinear form and linear functional stated above by where we have set H = L 2 (D 0 ) × L 2 (Γ 0 ) and V = H 1 (D 0 ) × H 1 (Γ 0 ) to be Hilbert spaces equipped with respective inner products The mean-weak formulation thus reads as follows (4.22) a((û, v), ( φ, ξ)) = l(( φ, ξ)).
The following uniform bounds on the random bulk coefficients follow immediately from the assumption (4.16) on the random domain mapping.Furthermore, the derived bounds on the surface coefficients presented in Proposition 4.1 also hold since the tangential derivatives of the surface parametrisation and its inverse are also uniformly bounded as a consequence of (4.16).

Proposition 4.2 (Uniform bounds).
There exists constants C g , C D > 0 such that the bulk area element g(ω) and the singular values σ i of D(ω) are uniformly bounded for all x ∈ D 0 and a.e.ω by Theorem 4.3.Given any ( f , fΓ0 ) ∈ H, there exist a unique solution (û, v) ∈ L 2 (Ω; V ) to (4.22) which satisfies the energy estimate Proof.With our uniform bounds (4.23), (4.10) on the random bulk and surface coefficients, we can now proceed in verifying all the conditions of the Lax-Milgram theorem are satisified.For a coercivity estimate, we argue . For the continuity of the bilinear form a(•, •), we apply the Cauchy-Schwarz inequality with the boundedness of the trace operator Thus we have the existence and uniqueness of a solution to (4.22).The estimate (4.25) then follows from coercivity of a(•, •).Theorem 4.4 (Regularity).Given any f ∈ L 2 (Ω; L 2 (D 0 )) and fΓ0 ∈ L 2 (Ω; L 2 (Γ 0 )), the mean-weak solution (û, v) to (4.22) satisfies Furthermore, we have where the constant C > 0 depends only the geometry of the reference domain D 0 and the uniform bound (4.16) on the random domain mapping.
Proof.Observe that for a.e.ω ∈ Ω, the solution (û, v) satisfies for every φ ∈ H 1 (D 0 ) and ξ ∈ H 1 (Γ 0 ), Setting φ = 0 gives Hence we see that v(ω) is the pathwise weak solution to the elliptic surface equation It therefore follows form the surface regularity result given in Theorem 4.2 since û(ω) ∈ L 2 (Γ 0 ), that v(ω) ∈ H 2 (Γ 0 ) for a.e.ω and furthermore where the constant C > 0 is independent of ω.To obtain higher regularity of the bulk quantity, we set ξ = 0 yielding This is precisely the weak formulation of the following elliptic boundary value problem subject to the reformulated Robin boundary condition Since the coefficients are sufficiently regular, more precisely and the boundary is sufficiently smooth Γ 0 ∈ C 2 , we can apply standard regularity results [22] to deduce û(ω) ∈ H 2 (D 0 ) for a.e.ω with the estimate Here the constant C > 0 is independent of ω since all the coefficients are uniformly bounded and furthermore, D(ω) is uniformly elliptic in ω.Combining (4.28) and (4.29) with the stability estimate (4.25) and boundedness of the trace operator leads to and hence the stated result.

An abstract numerical analysis of elliptic equations on random curved domains
We continue by considering in an abstract setting, the mean-weak formulation of general elliptic equations on random curved domains after being transformed onto the expected domain via the given stochastic domain mapping.Working in this abstract framework, we will present and analyse a finite element discretisation coupled with the Monte-Carlo method to approximate our quantity of interest, the mean solution.As the expected domain is assumed to be curved, the proposed finite element method will involve perturbations of the variational set up corresponding to the approximation of the domain.An optimal error bound in the energy norm for our non-conforming approach is derived with the help of the first lemma of Strang with suitable assumptions on the finite element space approximation and arising consistency error.Furthermore, an L 2 (Ω; L 2 )-type estimate is proved by a standard duality argument.
5.1.Abstract mean-weak formulation.Let V and H denote separable Hilbert spaces for which the embedding V → H is dense and continuous.We assume that we are in the setting where we have a sample dependent bilinear form ã(ω; •, •) : V × V → R and linear functional l(ω; •) : H → R corresponding to the path-wise weak formulation ã(ω; u(ω), ϕ) = l(ω; ϕ) of the elliptic equation after being reformulated onto the expected domain.For convenience, we will omit the pull-back notation for functions û since all the subsequent analysis will be considered on the expected domain.The mean-weak formulation will thus in general read as follows: Problem 5.1 (Mean-weak formulation).Find u ∈ L 2 (Ω; V ) such that for every ϕ ∈ L 2 (Ω; V ) we have We denote the associated bilinear form a(•, and shall assume all the requirements of the Lax-Milgram theorem are satisfied thus ensuring the existence and uniqueness of the solution. 5.2.Abstract formulation of the finite element discretisation.For a given h ∈ (0, h 0 ), let V h be a finite dimensional space that will represent a finite element space and let V h and H h denote the space V h endowed with respective norms • V h and • H h .We assume that V h and H h are Hilbert spaces and furthermore that V h → H h is uniformly embedded, that is for a constant c > 0 independent of h.In practice, the spaces V h and H h will represent equivalent Hilbert spaces to the continuous solution spaces V and H but posed over a discrete approximation of the curved domain, with h denoting the discretisation parameter.We introduce the sample-dependent bilinear form and linear functional ãh (ω; that are perturbations approximating their continuous counterparts and will assume ãh (ω : •, •) is uniformly V h -elliptic and bounded and additionally lh (ω; •) is uniformly bounded.More precisely, there exists constants c 1 , c 2 , c 3 > 0 independent of ω and h such that The finite element approximation of the mean-weak formulation (5.1) for a given a finite dimensional subspace V h ⊂ V h will then take the following form: Furthermore, for any w, ϕ ∈ L 2 (Ω; Z 0 ) with inverse lifts w −l , ϕ −l we have Our final assumption will be on the regularity of an associated dual problem that will enable us to derive an L 2 (Ω; H) error estimate using the standard Aubin-Nitsche trick.The associated dual problem reads as follows: Here (•, •) L 2 (Ω;H) denotes the inner product on the Hilbert space L 2 (Ω; H).
Assumption 5.4 (Regularity of dual problem).The solution w(g) to the dual problem belongs to space L 2 (Ω; Z 0 ) and furthermore satisfies for a constant c > 0 independent of both g and h ∈ (0, h 0 ).

5.4.
Error estimates for the semi-discrete solution.Recall that the abstract finite element space V h is not necessarily contained in the Hilbert space V .However, with the assumed existence of a lifting map (Ω; V ), we can lift the discrete bilinear form a h (•, •) and the linear functional l h (•) onto the space L 2 (Ω; V l h ) by the following relations for (5.10) thus inducing a third variational problem equivalent to (5.5).
Problem 5.4 (Lifted semi-discrete problem).Find u h ∈ L 2 (Ω; V l h ) such that for every ϕ h ∈ L 2 (Ω; V l h ) we have is contained in the solution space L 2 (Ω; V ), the lifted semi-discrete problem fits into the abstract non-conforming finite element setting considered in the first lemma of Strang [30].We will now present these results in the context of our random Hilbert space setting.Lemma 5.1 (First lemma of Strang).Let u h denote the solution to the lifted semi-discrete problem (5.11) and assume that the bilinear form for all ϕ ∈ L 2 (Ω; V l h ) and h ∈ (0, h 0 ).Then there exists a constant C > 0 independent of h such that Theorem 5.2 (Error estimates).Let u denote the solution of the continuous problem (5.1) and assume that it is sufficiently regular u ∈ L 2 (Ω; Z 0 ) and let U h be the discrete solution of (5.5) with lift u h = Λ h U h .Then with the assumptions listed in section 5.3 satisfied, there exists a constant c > 0 such that for all h ∈ (0, h 0 ) we have the error estimate Proof.It follows from the uniform ellipticity assumption (5.2) on the bilinear form a h (•, •) and the norm equivalence of the lifting map, that for any . Therefore the bilinear form a l h (•, •) is uniformly coercive and thus we can apply the first lemma of Strang.Substituting ϕ h = I h u into the estimate (5.12) and inserting the consistency bounds (P1), (P2) gives Hence with the interpolation estimate (I1) applied to u ∈ L 2 (Ω; Z 0 ) we obtain For the L 2 (Ω; H)−estimate, we use a standard duality argument.Given g ∈ L 2 (Ω; H) and an arbitrary Choosing w h = I h w(g) and applying the interpolation estimate (I1) to the solution of the dual problem which is assumed (R1) to be sufficiently regular w(g) ∈ L 2 (Ω; Z 0 ) gives We bound the consistency error in the second term with (P2) giving To obtain a bound of order h 2 for the third term, we begin by rewriting it as follows h (u, w(g)) .Now we are able to apply the estimate (P3) to the last term since both u, w(g) ∈ L 2 (Ω; Z 0 ) and can then follow a similar argument as to the previous cases for the first two terms which leads to Combining the results gives the stated result We conclude our abstract error analysis by combining our finite element discretisation with the Monte-Carlo method to estimate our quantity of interest, the mean solution E[u].Recall, that for an arbitrary Hilbert space H, the Monte-Carlo estimator of the expectation of a random variable Ŷi where M ∈ N is the chosen number of samples taken and Ŷi are independent identically distributed copies of the random variable Y .Furthermore, we have the following well-known convergence result, see [24].
Lemma 5.2 (Monte-Carlo convergence rate).For a given M ∈ N and a H-valued random variable Y ∈ L 2 (Ω; H), the Monte-Carlo estimator satisfies the convergence rate Therefore, if we consider the error between the mean solution E[u] and our discrete approximation E[u h ] in the L 2 (Ω M ; H) norm, and decompose it into the error arising from the finite element discretisation and the statistical error for the Monte-Carlo approximation, we obtain the following bound A similar argument in the L 2 (Ω; V ) leads to the following convergence rates.
Theorem 5.3.Let all the conditions from Theorem 5.2 be satisfied.Then we have the following error estimates (5.17)

Discretisation of the reformulated elliptic PDEs on their expected domains
In this section, we apply the results from the abstract theory to two finite element discretisation schemes for the reformulations of the two model elliptic equations.In each case, we will verify that all the listed assumptions in abstract setting are satisfied hence giving the stated convergence rate.
6.1.The elliptic equation on a random surface.To discretise the reformulation of the elliptic equation on the expected domain, we propose a semi-discrete scheme using linear Lagrangian surface finite elements [14].Our computational domain Γ h approximating the smooth expected hypersurface Γ 0 will be a polyhedral surface consisting of finitely many non-degenerate triangles whose vertices are taken to lie on the surface Γ 0 and have the maximum diameter bounded above by h > 0. The triangulation will be assumed to be shape regular and quasi-uniform, in the sense that the in-ball radius of each element is uniformly bounded below by ch, for some constant c > 0. In order to lift functions between the continuous and discrete surface, we shall assume that the projective mapping a : Γ h → Γ 0 decribed in (3.5) is bijective and define the lift and inverse lift of functions f and g given over Γ h and Γ 0 respectively by where x(a) denotes the inverse of the projection mapping a.We introduce the linear finite element space on Γ h (6.2) and define the lifted finite element space by (6.3) The finite element discretisation of the mean-weak formulation reads as follows.Problem 6.1 (Semi-discrete scheme).Find U h ∈ L 2 (Ω; S h ) such that (6.4) for every φ h ∈ L 2 (Ω; S h ).
In the context of the abstract framework, the finite dimensional space V h is taken to be the finite element space S h and the Hilbert spaces V h , H h are given by H 1 (Γ h ) and L 2 (Γ h ).Furthermore, the abstract sampledependent discrete bilinear form ãh (ω; With the uniform bounds on the random coefficients (4.10), (4.11), we deduce that ãh (ω : )-elliptic and bounded, and additionally l(ω; •) is uniformly bounded as presumed in (5.2 -5.4), and hence obtain existence and uniqueness of a semi-discrete solution to (6.4).We continue by checking the stated assumptions in the abstract error analysis.In particular, we begin with the norm equivalence (L1),(L2) of the lifting map Λ h : V h → V given by Λ h χ h = χ l h .A proof of these estimates can be found in [14, Lemma 4.2].Lemma 6.1 (Equivalence in norms of lifts).There exists constants c 1 , c 2 > 0 independent of h such that for any For the interpolation assumption (I1), we set the Hilbert space Z 0 consisting of functions of higher regularity to be H 2 (Γ 0 ).It follows from the Sobolev embedding that H 2 (Γ 0 ) ⊂ C 0 (Γ 0 ) for n ≤ 3 and therefore we can introduce the interpolation operator I h : H 2 (Γ 0 ) → S l h defined by (6.5) where Îh : C 0 (Γ h ) → S h denotes the standard Lagrangian interpolatant defined element-wise on Γ h .The following estimate was proved in [14,Lemma 4.3].
Lemma 6.2 (Interpolation estimate).Given any η ∈ H 2 (Γ 0 ), there exists a constant c > 0 independent of h such that To derive the assumed bounds (P1),(P2) and (P3) on the approximation of the discrete bilinear forms, we first need a preliminary result on the order of approximation of the geometry, see [14,Lemma 4.1].
Lemma 6.3 (Geometric error bounds).Let δ Γ0 h denote the surface element corresponding to the transformation from Γ 0 to Γ h under the lifting map dσ(a(x)) = δ h (x)dσ h (x) and define where P h := I − ν h ⊗ ν h is the projection operator mapping onto the tangent space of the discrete surface Γ h defined element-wise.Then we have the estimates We can now bound the consistency error as follows.

Lemma 6.4 (Consistency error). Given any (W
Proof.Lifting the discrete integral in the linear functional l h (•) onto the smooth surface Γ 0 with the projective mapping a(•) leads to Hence with the uniform bound (4.11) on the random coefficient g Γ0 (ω) and the order h 2 approximation of the geometric pertubation (6.9), we obtain the estimate (6.11).For (6.12), we begin by applying the chain rule to lift W h (ω, x) = w h (ω, a(x)) Suppressing the parameter x, we deduce . Therefore, we can express the pertubation error in the approximation of the bilinear form a(•, •) by and hence with the uniform bounds (4.10), (4.11) on the random coefficients and the geometric estimates (6.9), (6.10) we obtain (6.12).
6.2.The coupled elliptic system.We next apply the results from the abstract framework to the second model problem of the coupled elliptic system on a random bulk-surface.Our proposed finite element discretisation of the system reformulated on the expected domain and the subsequent analysis will be based on the approach presented in [12].For the computational domain, we approximate the open bulk D 0 ⊂ R n+1 by a polyhedral domain consisting of closed (n+1)−simplices with maximum diameter uniformly bounded above by positive constant h > 0 and will assume that the triangulation T h is quasi-uniform.We denote the induced discrete surface Γ h = ∂D h and the associated triangulation by and impose the same assumptions on T h as were listed in the previous example.A piece-wise diffeomorphic mapping G h : D h → D 0 from the discrete bulk to the continuous can be constructed by fixing the interior simplices (simplices with at most one vertex on the boundary Γ 0 ) and using the projective mapping a Γ0 (•) to define a diffeomorphism Λ h,k : K → K e between the boundary simplices K (simplices with at least two vertices on Γ 0 ) and the exact curved simplices K e , (6.13) Details on the precise form of Λ h,K can be found in [12].We are therefore able to define lifts and inverse lifts of functions on the bulk domain by Note that, the diffeomorphism Λ h,K is chosen such that the mapping G h coincides with the projective mapping (6.16) on the boundary of the discrete bulk and hence the bulk lift agrees with the surface lifting map described in (6.1) on ∂D h .For convenience, we will denote the sub-triangulation consisting of all boundary simplices by and define the corresponding sets (6.17) where the lifting maps G h , G −1 h differ from the identity mapping.We introduce the linear finite element spaces on the discrete bulk and discrete surface by and denote the corresponding lifted finite element spaces by (6.20) An important feature of our finite element spaces is that the trace of a function φ h ∈ V h belongs to S h and similarly the trace of ϕ h ∈ V l h belongs to S l h as a result of (6.16).The finite element discretisation of the mean-weak formulation then reads as follows.
Here the abstract finite dimensional space is V h = V h × S h and the Hilbert spaces V h , H h are given by H 1 (D 0 ) × H 1 (Γ 0 ) and L 2 (D 0 ) × L 2 (Γ 0 ) respectively.We denote the associated bilinear form and linear functional a to be the respective left hand side and right hand side of the semi-discrete variational problem 6.2.By the uniform bounds on the random coefficients (4.23), (4.10), we deduce the existence and uniqueness of a semi-discrete solution using a similar argument to the continuous problem.We proceed in a similar manner and check that the assumptions of the abstract analysis are satisfied.The norm equivalence (L1), (L2) of the lifting mapping which in this setting Λ h : , follows from the estimates on the surface lifting map given Lemma 6.1 in combination with the following bulk lifting norm equivalence derived in [12, Proposition 4.9].Lemma 6.5 (Bulk lift estimates).There exists constants c 1 , c 2 > 0 independent of h, such that for any For the interpolation assumption (I1), we set the abstract function space Z 0 = H 2 (D 0 ) × H 2 (Γ 0 ) and define the interpolation operator component-wise with Ĩh denoting the standard Lagrangian intepolation operator and have the following estimate .

Lemma 6.6 (Interpolation estimate).
There exists a well-defined interpolation operator The next step will entail bounding the consistency error arising from the geometric approximation of the domain.Estimates for the surface pertubation have previously been given in Lemma 6.3.For the bulk approximation, we recall that the lifting mapping G h : D h → D 0 is defined to be the identity on interior simplices and a C 1 −diffeomorphism for simplices near the boundary.Therefore the corresponding bulk error will be comprised of two parts; the first part will be related to the smallness of the neighbourhood around Γ 0 in which the lifted boundary simplices lie in and the second part is the associated geometric error of the boundary simplices approximating the corresponding exact curved simplex.We begin with the latter and state geometric bulk estimates on the diffeomorphic mapping G h , for which a proof of the bounds (6.24) and (6.
Then we have the following estimates for a constant c > 0 independent of ω, Proof.The estimate (6.26) follows from the observation Then for any η ∈ H 1 (D 0 ) we have The consistency error can now be bounded as follows.
Examining the bulk term, we see that we are unable to apply the narrow band inequality Lemma 6.8, to the derivative of ϕ h (ω) and w h (ω) since the functions only belong to the space V h ⊂ H 1 (D 0 ), resulting in the bound of order h given in (6.29).However, considering sufficiently regular functions ϕ, w ∈ L 2 (Ω; H 2 (D 0 )), we are able to employ Lemma 6.8 attaining the estimate of order h 2 given in (6.30).
The regularity assumption (R1) on the associated dual problem follows again from the symmetry of the bilinear for a(•, •) and the previously derived regularity result given in Theorem 4.4.Hence all the assumptions of the abstract theory are satisfied and we have the stated convergence rate given in Theorem 5.3.

Numerical results
In this section, we numerically verify the stated convergence rates of the two proposed finite element discretisations of the reformulated model elliptic problems.In both cases, the numerical scheme has been implemented in DUNE [3,9].denotes the spherical harmonic function of degree l and order m, which correspond to the eigenvalues of the Laplace-Beltrami operator.For further details on exact form of the spherical harmonics, we refer the reader to [1,15].Realisations of the random surface for different samples are given below in Figure 3.To numerical verify the convergence rate, we set the exact pull-back solution to be given by û(ω, x) = sin(π(x 2 − 1)y(z − 1)) + σ tol ν 1 (ω)cos(πz(y + 1)) + σ tol ν 2 (ω)sin(π(x + y)z 2 ) with ν 1 , ν 2 ∼ U (−1, 1) and σ tol > 0 a constant controlling the largest deviation of pathwise solution.This in turn determines the random data f given in the reformulated elliptic equation (4.4).We observe the following errors for the approximation E[û] − E M [û h ] in L 2 (Ω M ; L 2 (Γ 0 )) and L 2 (Ω M ; H 1 (Γ 0 )) and thus the stated convergence results.7.2.Random bulk-surface.For the coupled-elliptic system on a random bulk-surface, we adopt a similar approach to the random surface numerical example and prescribe the curved boundary to the random bulk D(ω) which for simplicity is taken to lie in R 2 , as a graph λ n (ω)cos(nθ) + λn (ω)sin(nθ) x = (cos(θ), sin(θ)) ∈ S 1 , with independent, uniformly distributed random coefficients λ n , λn ∼ U (−1, 1).We extend the given boundary process in the normal direction into the interior with a sufficiently smooth blending function to form the stochastic domain mapping (7.4) φ(x, ω) = x + L δ (|x − a Γ0 (x)|)h(a Γ0 (x), ω)ν Γ0 (a Γ0 (x)) x ∈ B 1 (0).
Here the precise form of the chosen blending function L δ (•) : R ≥0 → R ≥0 is given by Realisations of the image of the reference domain mappped under the random domain mapping (7.4) are provided in Figure 4. g Γ0 (ω) G −1 (ω)ν Γ0 • ∇û(ω) = 0 on Γ 0 , from which the data f and fΓ0 can then be computed.Note that in practice, the expectation E[v] and its surface derivative are approximated with Monte-Carlo sampling to sufficiently high accuracy.We observe the following errors and experimental order of convergence for the approximations of the bulk

Figure 1 .
Figure 1.A realisation of the random domain mapping.

Figure 2 .
Figure 2. The computational domain, parametrisation based expected domain and a realisation of the random domain.
and the uniform bounds (4.23) on the random coefficient D(ω).To obtain a bound on the open neighbourhood containing the boundary simplices, we have the subsequent narrow band inequality[12, Lemma 4.10].Lemma 6.8 (Narrow band trace inequality).Given any δ < δ Γ0 , let N δ be a narrow band in the interior domain D 0 around the boundary Γ 0 defined by(6.27)

Figure 3 .
Figure 3. Realisations of the path-wise solution on the associated realisation of the random surface.

Figure 4 .
Figure 4. Realisations on pathwise solution on the random bulk-surface.