Design and analysis of finite volume methods for elliptic equations with oblique derivatives; application to Earth gravity field modelling

TL;DR

This paper develops and analyzes finite volume methods for elliptic equations with oblique boundary conditions, providing a generic framework, convergence proofs, and practical applications including Earth gravity field modeling.

Contribution

It introduces a novel generic finite volume framework for oblique boundary conditions, with convergence analysis and specific flux schemes validated by numerical tests and real-world data.

Findings

Optimal convergence rate of O(h) in energy norm.

Validated methods through extensive 3D numerical tests.

Successful application to Earth gravity potential data.

Abstract

We develop and analyse finite volume methods for the Poisson problem with boundary conditions involving oblique derivatives. We design a generic framework, for finite volume discretisations of such models, in which internal fluxes are not assumed to have a specific form, but only to satisfy some (usual) coercivity and consistency properties. The oblique boundary conditions are split into a normal component, which directly appears in the flux balance on control volumes touching the domain boundary, and a tangential component which is managed as an advection term on the boundary. This advection term is discretised using a finite volume method based on a centred discretisation (to ensure optimal rates of convergence) and stabilised using a vanishing boundary viscosity. A convergence analysis, based on the 3rd Strang Lemma \cite{DPD18}, is conducted in this generic finite volume framework, and yields the expected $\mathcal O(h)$ optimal convergence rate in discrete energy norm. We then describe a specific choice of numerical fluxes, based on a generalised hexahedral meshing of the computational domain. These fluxes are a corrected version of fluxes originally introduced in \cite{Medla.et.al2018}. We identify mesh regularity parameters that ensure that these fluxes satisfy the required coercivity and consistency properties. The theoretical rates of convergence are illustrated by an extensive set of 3D numerical tests, including some conducted with two variants of the proposed scheme. A test involving real-world data measuring the disturbing potential in Earth gravity modelling over Slovakia is also presented.

Tables

Table 1. Table 1: The regularity parameters ( 24 ), ( 42 ), ( 44 )and ( 47 ) for a non-uniform mesh of the cube.

$h$	${\mathrm{reg}}_{𝔗}$	${\mathrm{reg}}_{𝔗,\Gamma }$	${\mathrm{reg}}_{𝔗,\Omega }$	${\varrho }_{𝔗,\Omega }$
8.511e-01	7.311	5.536	3.322	6.768e-01
3.660e-01	7.427	6.394	3.279	4.295e-01
1.685e-01	7.920	5.949	3.402	3.508e-01
8.309e-02	7.879	5.967	3.383	2.798e-01
4.084e-02	8.267	6.870	3.510	2.089e-01
2.041e-02	8.655	6.584	3.585	1.669e-01
1.014e-02	8.356	6.854	3.589	1.426e-01

Table 2. Table 2: EOC for the non-uniform mesh of the cube with 𝑽 ​ ( 𝒙 ) = ( − 1 , − 1 , − 1 ) 𝑽 𝒙 1 1 1 \boldsymbol{V}(\boldsymbol{x})=(-1,-1,-1)

Scheme (19)
$h$	$L_{\Omega }^2$ error	EOC	$L_{\Gamma }^2$ error	EOC	$V_h$ error	EOC	$V_{h,\Gamma }$ error	EOC
8.511e-01	2.092e-02		4.600e-02		1.886e-01		1.844e-01
3.660e-01	9.412e-03	0.946	2.634e-02	0.660	1.022e-01	0.726	1.462e-01	0.275
1.685e-01	3.922e-03	1.128	1.270e-02	0.940	4.767e-02	0.982	9.298e-02	0.583
8.309e-02	1.958e-03	0.982	6.756e-03	0.893	2.475e-02	0.927	6.358e-02	0.537
4.084e-02	1.003e-03	0.942	3.570e-03	0.898	1.294e-02	0.913	4.248e-02	0.567
2.041e-02	5.155e-04	0.959	1.867e-03	0.934	6.661e-03	0.957	2.678e-02	0.665
1.014e-02	2.560e-04	1.000	9.360e-04	0.986	3.287e-03	1.009	1.587e-02	0.747

Table 3. Table 3: EOC for the non-uniform mesh of the cube with 𝑽 ​ ( x , y , z ) = ( x , y , − 1 ) 𝑽 𝑥 𝑦 𝑧 𝑥 𝑦 1 \boldsymbol{V}(x,y,z)=(x,y,-1)

Scheme (19)
$h$	$L_{\Omega }^2$ error	EOC	$L_{\Gamma }^2$ error	EOC	$V_h$ error	EOC	$V_{h,\Gamma }$ error	EOC
8.478e-01	1.508e-02		2.908e-02		1.066e-01		1.021e-01
3.606e-01	9.891e-03	0.493	2.639e-02	0.113	8.779e-02	0.226	1.208e-01	-0.197422
1.694e-01	5.583e-03	0.757	1.671e-02	0.605	5.182e-02	0.697	9.645e-02	0.298902
8.197e-02	3.237e-03	0.751	1.001e-02	0.705	2.929e-02	0.786	7.008e-02	0.439947
4.068e-02	1.803e-03	0.835	5.652e-03	0.815	1.571e-02	0.889	4.541e-02	0.619154
2.032e-02	9.617e-04	0.905	3.036e-03	0.895	8.098e-03	0.954	2.699e-02	0.749373
1.031e-02	5.041e-04	0.951	1.598e-03	0.945	4.148e-03	0.985	1.534e-02	0.832171

Table 4. Table 4: EOC for the non-uniform mesh of the cube with 𝑽 ​ ( x , y , z ) = ( − x , z , − 1 ) 𝑽 𝑥 𝑦 𝑧 𝑥 𝑧 1 \boldsymbol{V}(x,y,z)=(-x,z,-1)

Scheme (19)
$h$	$L_{\Omega }^2$ error	EOC	$L_{\Gamma }^2$ error	EOC	$V_h$ error	EOC	$V_{h,\Gamma }$ error	EOC
8.594e-01	1.402e-02		3.084e-02		1.103e-01		1.070e-01
3.651e-01	7.731e-03	0.695	2.187e-02	0.401	7.795e-02	0.405	1.084e-01	-0.015
1.703e-01	4.769e-03	0.633	1.455e-02	0.534	4.839e-02	0.625	9.136e-02	0.224
8.271e-02	2.620e-03	0.828	8.360e-03	0.767	2.636e-02	0.840	6.458e-02	0.480
4.070e-02	1.389e-03	0.894	4.520e-03	0.867	1.347e-02	0.947	4.045e-02	0.659
2.037e-02	7.382e-04	0.913	2.428e-03	0.897	6.985e-03	0.948	2.454e-02	0.721
1.018e-02	3.768e-04	0.969	1.248e-03	0.960	3.490e-03	1.000	1.378e-02	0.832

Table 5. Table 5: The regularity parameters ( 24 ), ( 42 ), ( 44 ) and ( 47 ) for a non-uniform mesh of a tesseroid mesh.

$h$	${\mathrm{reg}}_{𝔗}$	${\mathrm{reg}}_{𝔗,\Omega }$	${\mathrm{reg}}_{𝔗,\Gamma }$	${\varrho }_{𝔗,\Omega }$
1.092	2.538e+01	5.468	3.296	6.074e-01
4.981e-01	1.336e+01	5.970	3.468	4.051e-01
2.375e-01	1.045e+01	6.228	3.297	3.691e-01
1.158e-01	1.017e+01	6.352	3.442	2.463e-01
5.733e-02	0.998e+01	6.324	3.426	1.460e-01
2.847e-02	1.029e+01	6.582	3.598	1.392e-01
1.454e-02	1.032e+01	6.761	3.580	8.204e-02

Table 6. Table 6: EOC for a non-uniform mesh of a tesseroid with 𝑽 ​ ( 𝒙 ) = ( 0.3 , 0.2 , 0.1 ) − 𝒙 𝑽 𝒙 0.3 0.2 0.1 𝒙 \boldsymbol{V}(\boldsymbol{x})=(0.3,0.2,0.1)-\boldsymbol{x}

Scheme (19)
$h$	$L_{\Omega }^2$ error	EOC	$L_{\Gamma }^2$ error	EOC	$V_h$ error	EOC	$V_{h,\Gamma }$ error	EOC
1.092	7.054e-03		9.644e-03		3.776e-02		3.620e-02
4.981e-01	1.307e-03	2.148	2.189e-03	1.890	9.516e-03	1.756	1.205e-02	1.402
2.375e-01	3.656e-04	1.720	7.681e-04	1.414	3.249e-03	1.451	5.267e-03	1.117
1.158e-01	1.294e-04	1.446	3.335e-04	1.161	1.322e-03	1.252	2.698e-03	0.931
5.733e-02	5.493e-05	1.219	1.578e-04	1.065	5.788e-04	1.175	1.400e-03	0.932
2.847e-02	2.518e-05	1.114	7.580e-05	1.047	2.643e-04	1.119	7.200e-04	0.950
1.454e-02	1.248e-05	1.044	3.843e-05	1.011	1.285e-04	1.073	3.805e-04	0.949

Table 7. Table 7: The regularity parameters ( 24 ), ( 42 ), ( 44 ) and ( 47 ) for a non-uniform mesh of a section of a perturbed ball.

$h$	${\mathrm{reg}}_{𝔗}$	${\mathrm{reg}}_{𝔗,\Omega }$	${\mathrm{reg}}_{𝔗,\Gamma }$	${\varrho }_{𝔗,\Omega }$
1.089	3.228e+01	5.017	3.571	3.068e-01
4.928e-01	1.896e+01	6.270	3.551	-6.303e-01
2.299e-01	1.207e+01	6.259	3.777	-9.516e-01
1.159e-01	1.120e+01	6.819	3.838	-1.138
5.813e-02	1.083e+01	7.246	3.904	-1.524
2.867e-02	1.114e+01	6.890	4.253	-1.403
1.462e-02	1.118e+01	6.914	4.232	-1.642

Table 8. Table 8: EOC for Ω Ω \Omega the section of a perturbed ball. 𝑽 ​ ( 𝒙 ) = ( 0.3 , 0.2 , 0.1 ) − 𝒙 𝑽 𝒙 0.3 0.2 0.1 𝒙 \boldsymbol{V}(\boldsymbol{x})=(0.3,0.2,0.1)-\boldsymbol{x}

Scheme (19)
$h$	$L_{\Omega }^2$ error	EOC	$L_{\Gamma }^2$ error	EOC	$V_h$ error	EOC	$V_{h,\Gamma }$ error	EOC
1.089	7.577e-03		1.212e-02		5.851e-02		6.169e-02
4.928e-01	2.595e-03	1.351	6.735e-03	0.740	4.232e-02	0.408	7.120e-02	-0.180
2.299e-01	1.519e-03	0.7021	5.123e-03	0.358	2.766e-02	0.557	6.351e-02	0.150
1.159e-01	8.678e-04	0.817	3.163e-03	0.704	1.465e-02	0.928	4.295e-02	0.571
5.813e-02	4.763e-04	0.868	1.776e-03	0.835	7.250e-03	1.019	2.550e-02	0.754
2.867e-02	2.548e-04	0.885	9.569e-04	0.875	3.594e-03	0.993	1.433e-02	0.815
1.462e-02	1.335e-04	0.959	5.028e-04	0.955	1.793e-03	1.032	7.777e-03	0.907

Table 9. Table 9: EOC for Ω Ω \Omega a cube and 𝑽 ​ ( 𝒙 ) = ( 11.4301 , 0 , − 1 ) 𝑽 𝒙 11.4301 0 1 \boldsymbol{V}(\boldsymbol{x})=(11.4301,0,-1)

Scheme (19)
$h$	$L_{\Omega }^2$ error	EOC	$L_{\Gamma }^2$ error	EOC	$V_h$ error	EOC	$V_{h,\Gamma }$ error	EOC
8.605e-01	2.247e-02		4.945e-02		1.686e-01		1.610e-01
3.606e-01	1.537e-02	0.436	4.232e-02	0.179	1.437e-01	0.184	2.014e-01	-0.257
1.692e-01	1.088e-02	0.455	3.301e-02	0.328	1.060e-01	0.401	2.012e-01	0.001
8.347e-02	6.582e-03	0.712	2.109e-02	0.634	6.669e-02	0.655	1.647e-01	0.283
4.054e-02	3.756e-03	0.776	1.255e-02	0.718	4.010e-02	0.704	1.248e-01	0.384
2.043e-02	2.046e-03	0.886	7.078e-03	0.836	2.333e-02	0.790	9.148e-02	0.453
1.017e-02	1.100e-03	0.889	3.913e-03	0.848	1.342e-02	0.792	6.596e-02	0.468

Table 10. Table 10: The GNSS-leveling test [ m ] delimited-[] 𝑚 [m] at 58 points in area of Slovakia.

	method 19	method 48	method 49
Min	0.229	0.237	0.239
Mean	0.326	0.330	0.331
Max	0.449	0.458	0.459
Range	0.22	0.221	0.220
STD	0.052	0.050	0.050

Equations

$$-\Delta \overline{T}(\vec{x})=0,$$

$$\nabla \overline{T}(\vec{x})\cdot \vec{V}(\vec{x})=g(\vec{x}),$$

$$\overline{T}(\vec{x})=0,$$

$$\vec{V}(\vec{x})=\vec{n}(\vec{x})+\vec{W}(\vec{x}),\forall \vec{x}\in \Gamma .$$

$$\Delta G(\vec{x})=2{\omega }^2,$$

$$|\nabla G(\vec{x})|=\overline{g}(\vec{x}),$$

$$G(\vec{x})=U(\vec{x})+\overline{T}(\vec{x}).$$

$$\Delta \overline{T}(\vec{x})=0,$$

$$|\nabla (\overline{T}+U)(\vec{x})|=\overline{g}(\vec{x}),$$

$$\overline{T}(\vec{x})={\overline{T}}_{\text{mbox}sat}(\vec{x}),$$

$$\frac{\nabla (\overline{T}+U)(\vec{x})}{|\nabla (\overline{T}+U)(\vec{x})|}\cdot \nabla (\overline{T}+U)(\vec{x})=\overline{g}(\vec{x}).$$

$$\vec{v}(\vec{x})=\frac{\nabla (\overline{T}+U)(\vec{x})}{|\nabla (\overline{T}+U)(\vec{x})|}$$

$$\vec{v}(\vec{x})\cdot \nabla \overline{T}(\vec{x})=\overline{g}(\vec{x})-\vec{v}(\vec{x})\cdot \nabla U(\vec{x}),\vec{x}\in \Gamma .$$

$$\vec{V}(\vec{x})\cdot \nabla \overline{T}(\vec{x})=\overline{g}(\vec{x})-\overline{\gamma }(\vec{x}),$$

$${\overline{T}}_n(\vec{x})+(\alpha \overline{T}{)}_t(\vec{x})=g(\vec{x}),$$

$$\displaystyle\mbox{Each control volume $p\in\mathfrak{T}$ has at most one face $\sigma$ in $\mathfrak{S}_{\Gamma}$ and, in that case, $\vec{x}_{p}\in\sigma$}.$$

$$V_{h}:=\{\varphi=((\varphi_{p})_{p\in\mathfrak{T}},(\varphi_{\sigma})_{\sigma\in\mathfrak{S}_{\rm Dir}},(\varphi_{e})_{e\in\mathfrak{E}_{\Gamma}})\,:\,\varphi_{p}\in\mathbb{R}\,,\;\varphi_{\sigma}=0\,,\;\varphi_{e}\in\mathbb{R}\,,\;\varphi_{e}=0\mbox{ if $e\not\in\mathfrak{E}_{\Gamma\!,{\rm int}}$}\}.$$

$$\Vert \phi {\Vert }_{V_h}:={\left(|\phi {|}_{V_h,\Omega }^2+h_{\Gamma }|\phi {|}_{V_h,\Gamma }^2\right)}^{1/2},$$

$$|\phi {|}_{V_h,\Omega }:={\left(\sum _{\sigma \in \mathfrak{S}\backslash {\mathfrak{S}}_{\Gamma }}\frac{|\sigma |}{d_{pq}}({\phi }_p-{\phi }_q{)}^2\right)}^{1/2},$$

$$|\phi {|}_{V_h,\Gamma }:={\left(\sum _{\sigma \in {\mathfrak{S}}_{\Gamma }}\sum _{e\in \mathfrak{E}(\sigma )}\frac{|e|}{d_{pe}^{\perp }}({\phi }_p-{\phi }_e{)}^2\right)}^{1/2}$$

$$I_h\phi =(({\phi }_p{)}_{p\in \mathfrak{T}},({\phi }_{\sigma }{)}_{\sigma \in {\mathfrak{S}}_{\mathrm{Dir}}},({\phi }_e{)}_{e\in {\mathfrak{E}}_{\Gamma }})\text{mbox}with$$

$${\phi }_p=\phi ({\vec{x}}_p)\forall p\in \mathfrak{T},{\phi }_{\sigma }=\phi ({\vec{x}}_q)\forall q=\sigma \in {\mathfrak{S}}_{\mathrm{Dir}},{\phi }_e=\frac{1}{|e|}{\int }_e\phi \forall e\in {\mathfrak{E}}_{\Gamma }.$$

$$0={\iiint }_p-\Delta \overline{T}=$$

$$=$$

$$0=\sum _{\sigma \in \mathfrak{S}(p)\backslash {\mathfrak{S}}_{\Gamma }}{\overline{F}}_{p,\sigma }(\overline{T})-\sum _{\sigma \in \mathfrak{S}(p)\cap {\mathfrak{S}}_{\Gamma }}{\iint }_{\sigma }g-\nabla \overline{T}\cdot \vec{W}.$$

$$\sum _{\sigma \in \mathfrak{S}(p)\cap {\mathfrak{S}}_{\Gamma }}{\iint }_{\sigma }g=$$

$$=$$

$$\sum _{\sigma \in \mathfrak{S}(p)\backslash {\mathfrak{S}}_{\Gamma }}{\overline{F}}_{p,\sigma }(\overline{T})+\sum _{\sigma \in \mathfrak{S}(p)\cap {\mathfrak{S}}_{\Gamma }}\sum _{e\in \mathfrak{E}(\sigma )}{\left[\overline{T}\vec{W}\cdot \vec{n}\right]}_{\sigma ,e}-\sum _{\sigma \in \mathfrak{S}(p)\cap {\mathfrak{S}}_{\Gamma }}{\left[\overline{T}{\nabla }_{\Gamma }\cdot \vec{W}\right]}_{\sigma }=\sum _{\sigma \in \mathfrak{S}(p)\cap {\mathfrak{S}}_{\Gamma }}{\iint }_{\sigma }g.$$

$${\left[\overline{T}\vec{W}\cdot \vec{n}\right]}_{\sigma ,e}\approx T_e{\left[\vec{W}\cdot \vec{n}\right]}_{\sigma ,e}\text{mbox}and{\left[\overline{T}{\nabla }_{\Gamma }\cdot \vec{W}\right]}_{\sigma }\approx T_p{\left[{\nabla }_{\Gamma }\cdot \vec{W}\right]}_{\sigma },$$

$${\mathcal{F}}_{p,\sigma }^{\Omega }(\phi )+{\mathcal{F}}_{q,\sigma }^{\Omega }(\phi )=0.$$

$$\sum _{\sigma \in \mathfrak{S}(p)\backslash {\mathfrak{S}}_{\Gamma }}{\mathcal{F}}_{p,\sigma }^{\Omega }(T)+\sum _{\sigma \in \mathfrak{S}(p)\cap {\mathfrak{S}}_{\Gamma }}$$

$$+R{h}_{\Gamma }\sum _{\sigma \in \mathfrak{S}(p)\cap {\mathfrak{S}}_{\Gamma }}\sum _{e\in \mathfrak{E}(\sigma )}{\mathcal{F}}_{\sigma ,e}^{\Gamma }(T)=\sum _{\sigma \in \mathfrak{S}(p)\cap {\mathfrak{S}}_{\Gamma }}{\iint }_{\sigma }g,\forall p\in \mathfrak{T},$$

$$\mathcal{F}^{\Gamma}_{\sigma,e}(T)+\mathcal{F}^{\Gamma}_{\tau,e}(T)=0,\quad\forall e\in\mathfrak{E}_{\Gamma\!,{\rm int}}\mbox{ with $\sigma,\tau\in\mathfrak{S}_{\Gamma}$ the two faces on each side of $e$}.$$

$$\sum _{\sigma \in \mathfrak{S}\backslash {\mathfrak{S}}_{\Gamma }}{\mathcal{F}}_{p,\sigma }^{\Omega }(\phi )\left({\phi }_p-{\phi }_q\right)\ge {\rho }_{\Omega }|\phi {|}_{V_h,\Omega }^2,$$

$$\sum _{\sigma \in {\mathfrak{S}}_{\Gamma }}\sum _{e\in \mathfrak{E}(\sigma )}{\mathcal{F}}_{\sigma ,e}^{\Gamma }(\phi )\left({\phi }_p-{\phi }_e\right)\ge {\rho }_{\Gamma }|\phi {|}_{V_h,\Gamma }^2.$$

$$\left|{\mathcal{F}}_{p,\sigma }^{\Omega }(I_{h}u)+{\iint }_{\sigma }\nabla u\cdot {\vec{n}}_{p,\sigma }\right|\le C_{\mathrm{cons}}^{\Omega }h|\sigma |\Vert u{\Vert }_{C^2(\overline{\Omega })},\forall p\in \mathfrak{T},\forall \sigma \in \mathfrak{S}(p)\backslash {\mathfrak{S}}_{\Gamma },$$

Full text

We then describe a specific choice of numerical fluxes, based on a generalised hexahedral meshing of the computational domain. These fluxes are a corrected version of fluxes originally introduced in [29]. We identify mesh regularity parameters that ensure that these fluxes satisfy the required coercivity and consistency properties. The theoretical rates of convergence are illustrated by an extensive set of 3D numerical tests, including some conducted with two variants of the proposed scheme. A test involving real-world data measuring the disturbing potential in Earth gravity modelling over Slovakia is also presented.

1. 1 Introduction

2. 2 Generic finite volume scheme

3. 2.1 Mesh, space of unknowns and interpolant

4. 2.2 Prolegomena to the scheme

5. 2.3 Scheme

6. 2.4 Error estimate

7. 3 Numerical tests

8. 3.1 Description of the scheme

9. 3.1.1 Inner fluxes

10. 3.1.2 Surface fluxes on $\Gamma$

11. 3.1.3 Properties of the fluxes

12. 3.2 Alternative schemes

13. 3.3 Results

14. 3.3.1 Cubic domain and non-uniform mesh

15. 3.3.2 Tesseroid domain with non-planar $\Gamma$

16. 3.3.3 Spherical section domain with perturbed bottom $\Gamma$

17. 3.3.4 Cubic domain, almost tangential vector field $\vec{V}$

18. 3.4 Local gravity field modelling

19. 4 Conclusion

20. A Proof of Theorem 6

21. A.1 Coercivity

22. A.2 Consistency

23. B Proof of Proposition 12

24. B.1 Boundary fluxes

25. B.2 Inner fluxes

26. B.2.1 Coercivity

27. B.2.2 Consistency

1 Introduction

We consider in this work a Laplace equation with oblique boundary conditions:

[TABLE]

where $\Omega$ is a bounded domain in $\mathbb{R}^{3}$ with piecewise $C^{2}$ boundary, $\Gamma$ is a relatively open subset of $\partial\Omega$ that is fully contained in a smooth component of this boundary, $g\in L^{2}(\partial\Omega)$ and $\vec{V}$ is a $C^{1}$ vector field such that $\vec{V}(\vec{x})\cdot\vec{n}(\vec{x})\not=0$ for all $\vec{x}\in\Gamma$. Here, $\vec{n}$ denotes the outer normal to $\partial\Omega$. We also assume that the $(d-1)$-dimensional measure of $\partial\Omega\backslash\Gamma$ is non-zero. On $\Gamma$, $\vec{V}$ can be decomposed into a normal and a tangential component to $\Gamma$. After renormalising $g$ we can assume that the normal component is $\vec{n}$, and thus that

[TABLE]

The properties of $\Gamma$ ensure that $\vec{W}$ is a $C^{1}$ tangential vector field on $\Gamma$.

A motivation to study boundary value problem (BVP) (1) comes from Earth gravity field modelling. The Earth gravity potential $G$ fulfils outside the Earth a non-homogeneous elliptic equation

[TABLE]

where $\omega$ is the spin velocity of the Earth [21]. The magnitude of the total gravity vector $\nabla G$ is called gravity. If the measured gravity is prescribed on the Earth surface, i.e.

[TABLE]

then Eq. (3) with BC (4) represents the so-called nonlinear geodetic BVP for the actual gravity potential $G$. The existence, uniqueness and other properties to the solution of this problem, and its variants, were studied extensively in physical geodesy community, see e.g. [1, 18, 24, 3, 17, 30, 20, 10, 11].

In Earth gravity field modelling, the actual gravity field $G$ is usually expressed as a sum of the selected model field $U$ and the remainder $\overline{T}$, i.e.

[TABLE]

If the model field $U$ is generated by an ellipsoid with the same mass as the Earth, rotating with the same spin velocity $\omega$ and with the constant surface potential equal to the geopotential $W_{0}$ (see [31] for a definition of $W_{0}$, the potential $G$ on the mean sea surface level), then $U$ is called the normal gravity potential and $\overline{T}$ is called the disturbing potential. This potential $\overline{T}$ has no centrifugal component and it is generally accepted that the disturbing potential satisfies the Laplace equation $\Delta\overline{T}=0$ outside the Earth, see e.g. [21, 22]. In the satellite era, people have been able to consider a bounded domain $\Omega$ outside the Earth where an upper part of the boundary is given as a sphere at altitude of a chosen satellite mission, and the bottom part $\Gamma\subset\partial\Omega$ is given by a subset of the Earth surface [16, 29]. On this bottom part $\Gamma$ the nonlinear BC (4) is given and, on the upper part, as well as on the side boundaries if one focuses on a tesseroid above the Earth, the Dirichlet-type BC obtained from satellite gravity missions can be prescribed. This allows us to fix a solution to the satellite data $\overline{T}_{\mbox{\small{sat}}}$. See Figure 1 for an illustration.

Such nonlinear satellite-fixed geodetic BVP [28] for the disturbing potential $\overline{T}$ can be formulated as follows

[TABLE]

Using the relation $|\vec{\xi}|=\frac{\vec{\xi}}{|\vec{\xi}|}\cdot\vec{\xi}$ for $\vec{\xi}\in\mathbb{R}^{3}\backslash\{0\}$ we get (7) in the form

[TABLE]

Letting

[TABLE]

be the actual gravity vector $\nabla G(\vec{x})$ unit direction, we can rewrite the nonlinear boundary condition (7) as follows

[TABLE]

Since the unit vector $\vec{v}(\vec{x})$ is unknown and depends on $\overline{T}(\vec{x})$, boundary condition (10) is still nonlinear. However, if we set $\overline{T}(\vec{x})=0$ in (9), which means that we approximate the unit direction $\vec{v}(\vec{x})$ of the actual gravity vector by the unit direction of the normal gravity vector equal to $\frac{\nabla U(\vec{x})}{|\nabla U(\vec{x})|}$, we get a linear(ized) boundary condition

[TABLE]

where $\overline{\gamma}(\vec{x})=\vec{V}(\vec{x})\cdot\nabla U(\vec{x})=\frac{\nabla U(\vec{x})}{|\nabla U(\vec{x})|}\cdot\nabla U(\vec{x})=|\nabla U(\vec{x})|$ is the so-called normal gravity. Since all quantities depending on $U$ are given analytically, the equation (11) represents a linear oblique derivative boundary condition. Together with equations (1a) and (1c), they are called the fixed gravimetric boundary value problem in the geodetic community [1, 24, 22, 23, 8, 16, 29] and give a basis for determining the Earth gravity field when gravity measurements are known on the Earth surface. When we denote $\overline{g}-\overline{\gamma}=g$ and consider the problem on a bounded domain outside the Earth, we end up with the oblique derivative BVP (1) (here, $\vec{V}$ has unit length and, as previously mentioned, the decomposition (2) is obtained after rescaling $g$).

Let us briefly mention some results that can be found in the literature regarding the numerical approximation of second order equations with oblique derivative boundary conditions. In [5, 4] authors deal with the finite volume method for the oblique derivative boundary value problem in 2D case. In [5] they consider the oblique BC in the form

[TABLE]

where $\alpha$ is a smooth function, $\overline{T}_{n}$ is a derivative in the normal direction and $(\alpha\overline{T})_{t}$ is a derivative in tangential direction. They develop a finite volume scheme based on the upwind principle, prove its convergence and obtain an error estimate of order $\sqrt{h}$. In [4], convergence results are established for the Poisson and a parabolic equation with oblique derivative boundary condition in which $\alpha$ is constant. The convergence results are not only obtained for the approximate finite volume solutions, but also for their discrete gradients. The error estimate of order $\sqrt{h}$ is obtained theoretically, but numerical experiments presented in these works indicate a first order rate of convergence. In [2] the authors present and analyse a 2D finite element method for the oblique derivative boundary problem, where the oblique derivative boundary is given by a graph of a real function. Finite volume methods for solving oblique derivative problems in 3D domains were suggested and numerically investigated in [26, 27, 25, 29]. These schemes are based either on upwind or central approximation of the oblique derivative. A numerical approximation of a nonlinear problem with eikonal-type boundary condition (4) was presented in [28]. This approximation is based on an iterative update of the oblique derivative condition.

In this paper we introduce and analyse novel numerical scheme for solving 3D oblique derivative boundary value problems for the Laplace equation. For the first time, there is presented a convergence analysis and error estimates for a finite volume scheme solving the oblique derivative problem in 3D. The model comes from the Earth gravity field modelling on real Earth topography but can be used in other applications as well. The presented numerical approach is general and covers various possible discretisations of the Laplace equation inside the domain and it treats in a robust and stable way the oblique derivative boundary condition. We present numerical tests showing convergence properties of the novel scheme and compare them to further alternative numerical treatments of the oblique derivative. We also present a local Earth gravity field modelling for the region of Slovakia where we compare obtained numerical results with GPS-leveling measurements.

The paper is organised as follows. In Section 2 we present a generic finite volume method for solving the oblique derivative problem for the Laplace equation and its error estimates. In Section 3 we specify approximation of inner and surface fluxes and present results of numerical computations. We also give alternative schemes for oblique derivative treatment and discuss their pros and cons. In Section 4 we present concluding remarks. Appendix A contains proof of error estimate for the generic scheme and Appendix B contains proof of coercivity and consistency of the suggested inner flux approximation.

2 Generic finite volume scheme

We describe here a generic finite volume approximation of (1). The discretisation is based on a recasting of the model to transform the oblique derivative into a normal derivative, handled as a Neumann boundary condition, and a boundary advection–reaction term along $\Gamma$. The method is “generic” in the sense that we do not impose any specific expression of the numerical fluxes, only broad assumptions that enable the convergence analysis of the method. Our approach and analysis therefore cover many possible choices of Finite Volume methods for discretising the Laplacian in the domain.

2.1 Mesh, space of unknowns and interpolant

Let $\mathfrak{T}$ be a partition of $\Omega$ into “generalised” polyhedral finite volumes $p$, the generalisation coming from the fact that the faces of the polyhedra could be curved (especially those lying on $\Gamma$). The mesh size is $h:=\max_{p\in\mathfrak{T}}{\rm diam}(p)$. We denote by $\mathfrak{S}$ the set of faces of the mesh, and by $\mathfrak{S}_{\rm int}$ the faces contained in $\Omega$. The boundary faces are assumed to be compatible with $\Gamma$ in the sense that each face in $\mathfrak{S}\backslash\mathfrak{S}_{\rm int}$ totally lies on $\Gamma$, or totally lies on the Dirichlet boundary $\partial\Omega\backslash\Gamma$. We let $\mathfrak{S}_{\Gamma}$ be the set of faces on $\Gamma$, and $\mathfrak{S}_{\rm Dir}=\mathfrak{S}\backslash(\mathfrak{S}_{\rm int}\cup\mathfrak{S}_{\Gamma})$ be the set of faces on $\partial\Omega\backslash\Gamma$.

For each cell $p\in\mathfrak{T}$ we take a point $\vec{x}_{p}\in p$ and we denote by $\mathfrak{S}(p)$ the set of faces of $p$, so that $\partial p=\cup_{\sigma\in\mathfrak{S}(p)}\sigma$. If $\sigma\in\mathfrak{S}(p)$, $\vec{n}_{p,\sigma}$ is the unit outer normal to $p$ on $\sigma$. Every face $\sigma$ in $\mathfrak{S}_{\Gamma}$ is a face of a unique finite volume $p$; the dependency of $p$ on $\sigma$ is not made explicit as there is no risk of confusion in the formulas. We assume that:

[TABLE]

Remark 1* (Assumption (13)).*

This assumption is not mandatory, and the design and analysis in the following sections could be adapted to meshes not satisfying (13) (see Remark 9); however, the method we consider in Section 3 naturally satisfies this property, which is why we assume it in our analysis.

For every $p\in\mathfrak{T}$ and $\sigma\in\mathfrak{S}(p)\cap\mathfrak{S}_{\rm int}$, we denote by $q_{p}(\sigma)$ the finite volume such that $\sigma=\overline{p}\cap\overline{q_{p}(\sigma)}$; here too, no risk of confusion arising we simply denote $q$ for $q_{p}(\sigma)$. We then set $d_{pq}=|\vec{x}_{p}-\vec{x}_{q}|$. A face $\sigma\in\mathfrak{S}(p)\cap\mathfrak{S}_{\rm Dir}$ on the Dirichlet boundary of a cell $p$ is sometimes considered as a “degenerate” cell, and also denoted by $q$; for such faces, we pick a point $\vec{x}_{q}\in\sigma$ and define again $d_{pq}=|\vec{x}_{p}-\vec{x}_{q}|$.

For each $\sigma\in\mathfrak{S}_{\Gamma}$ we take $\vec{x}_{\sigma}\in\sigma$, and we denote by $\mathfrak{E}(\sigma)$ the set of edges $e$ of $\sigma$. The set of all such edges is $\mathfrak{E}_{\Gamma}=\cup_{\sigma\in\mathfrak{S}_{\Gamma}}\mathfrak{E}(\sigma)$, and the edges that lie in the relative interior of $\Gamma$ are gathered in the set $\mathfrak{E}_{\Gamma\!,{\rm int}}$. For $e\in\mathfrak{E}(\sigma)$, $\vec{n}_{\sigma,e}$ is the unit normal outward to $\sigma$ on $e$ in the tangent space of $\Gamma$.

If $X$ is a control volume $p$, a face $\sigma$ or an edge $e$, $|X|$ denotes the Lebesgue measure of $X$ in the corresponding dimension of $X$ (dimension 3 for a control volume, 2 for a face, 1 for an edge).

Our space of approximation has unknowns in the finite volumes, on the Dirichlet faces (“degenerate cells”), and on each edge on $\Gamma$, with zero values for Dirichlet faces, and for edges that are not in the relative interior of $\Gamma$:

[TABLE]

Remark 2*.*

Introducing the zero-valued unknowns is of course not necessary, but will be useful to simplify some expressions.

The norm on $V_{h}$ is defined by

[TABLE]

where $d_{pe}^{\bot}$ is the orthogonal distance between $\vec{x}_{p}$ (which belongs to $\sigma$) and $e$. Remember that, in (14b), $p$ and $q$ are the two cells on each side of $\sigma$ if $\sigma\in\mathfrak{S}_{\rm int}$, and $q=\sigma$ if $\sigma\in\mathfrak{S}_{\rm Dir}$ (so that $\varphi_{q}=\varphi_{\sigma}=0$ in that case). The term $|\varphi|_{V_{h},\Omega}$ can thus be viewed as a discrete $H^{1}_{0}$-(semi)norm in $\Omega$ [13, Eq. (7.7f)], whilst $|\varphi|_{V_{h},\Gamma}$ plays the role of a discrete $H^{1}_{0}$-(semi)norm on the surface $\Gamma$. The presence of this boundary semi-norm, and its scaling by $h_{\Gamma}$, will be justified by the introduction of a small amount of diffusion on that surface to stabilise a centred approximation of an advective term on $\Gamma$ stemming from the oblique boundary condition (see (19a)). Notice that, in (14c), Assumption (13) was used to identify the unknown on a face $\sigma\in\mathfrak{S}_{\Gamma}$ with the value $\varphi_{p}$ corresponding to $p\in\mathfrak{T}$ such that $\sigma\in\mathfrak{S}(p)$.

The unknowns in the control volumes $p$ are destined to be approximations of the solution at $\vec{x}_{p}$, whereas those on the boundary edges approximate the average value of the solution on the corresponding edge. This leads to defining the following interpolant $I_{h}:C(\overline{\Omega})\to V_{h}$: for $\varphi\in C(\overline{\Omega})$ such that $\varphi=0$ on $\partial\Omega\backslash\Gamma$,

[TABLE]

The boundary condition $\varphi=0$ on $\partial\Omega\backslash\Gamma$ ensures that $\varphi_{\sigma}=0$ for all $\sigma\in\mathfrak{S}_{\rm Dir}$, and that $\varphi_{e}=0$ whenever $e\not\in\mathfrak{E}_{\Gamma\!,{\rm int}}$.

2.2 Prolegomena to the scheme

Integrating (1b) over a control volume $p\in\mathfrak{T}$, using Green’s theorem and introducing $\vec{W}$ defined in (2), it holds

[TABLE]

Denoting by $\overline{F}_{p,\sigma}(\overline{T})=-\iint_{\sigma}\nabla\overline{T}\cdot\vec{n}_{p,\sigma}\mathrm{d}\vec{x}$ the exact fluxes, we invoke the boundary condition (1b) to write

[TABLE]

The vector field $\vec{W}$ is tangential to $\Gamma$ and thus only the tangential gradient of $\overline{T}$ is involved in the quantity $\nabla\overline{T}\cdot\vec{W}$. We can therefore write $\nabla\overline{T}\cdot\vec{W}=\nabla_{\!\Gamma}\!\boldsymbol{\cdot}\!(\overline{T}\vec{W})-\overline{T}\nabla_{\!\Gamma}\!\boldsymbol{\cdot}\!\vec{W}$, where $\nabla_{\!\Gamma}\!\boldsymbol{\cdot}\!$ is the divergence operator on the manifold $\Gamma$. This leads, using the divergence theorem on each face $\sigma\in\mathfrak{S}(p)\cap\mathfrak{S}_{\Gamma}$, to

[TABLE]

Let us denote by $\left[\,\overline{T}\vec{W}\cdot\vec{n}\right]_{\sigma,e}=\int_{e}\overline{T}\vec{W}\cdot\vec{n}_{\sigma,e}$ the exact advection fluxes on the boundary, and by $\left[\,\overline{T}\nabla_{\!\Gamma}\!\boldsymbol{\cdot}\!\vec{W}\right]_{\sigma}=\iint_{\sigma}\overline{T}\nabla_{\!\Gamma}\!\boldsymbol{\cdot}\!\vec{W}\mathrm{d}\vec{x}$ the other contribution (akin to a reaction term) to the boundary term. This shows that the solution to (1) satisfies, for all $p\in\mathfrak{T}$,

[TABLE]

2.3 Scheme

The scheme for (1) is obtained discretising (16). As previously mentioned, we will assume generic properties on the diffusive numerical fluxes. The advective contribution to the boundary terms is discretised using a centred scheme, to which we add a small amount of (boundary) diffusion for stabilisation purposes. As discussed in Remark 16, the choice of a centred discretisation seems crucial to prove optimal error estimates.

Based on our choice of unknowns and interpolant (15), we make the following approximation, in which $T=((T_{p})_{p\in\mathfrak{T}},(T_{\sigma})_{\sigma\in\mathfrak{S}_{\rm Dir}},(T_{e})_{e\in\mathfrak{E}_{\Gamma}})$ is the sought approximation of $\overline{T}$:

[TABLE]

where $\left[\vec{W}\cdot\vec{n}\right]_{\sigma,e}=\int_{e}\vec{W}\cdot\vec{n}_{\sigma,e}$ and $\left[\nabla_{\!\Gamma}\!\boldsymbol{\cdot}\!\vec{W}\right]_{\sigma}=\int_{\sigma}\nabla_{\!\Gamma}\!\boldsymbol{\cdot}\!\vec{W}$. Here, we used Assumption (13) to utilise $T_{p}$ as approximate value of $\overline{T}$ on $\sigma\in\mathfrak{S}(p)\cap\mathfrak{S}_{\Gamma}$. The exact fluxes $\overline{F}_{p,\sigma}(\overline{T})$, for $p\in\mathfrak{T}$ and $\sigma\in\mathfrak{S}(p)\backslash\mathfrak{S}_{\Gamma}$, are discretised into numerical fluxes $\mathcal{F}^{\Omega}_{p,\sigma}(T)$ that satisfy the following conservativity condition: for all $\varphi\in V_{h}$ and all $\sigma\in\mathfrak{S}_{\rm int}$,

[TABLE]

We also select numerical diffusion fluxes $\mathcal{F}^{\Gamma}_{\sigma,e}(T)$ on the boundary, approximations of $-\int_{e}\nabla_{\Gamma}\overline{T}\cdot\vec{n}_{\sigma,e}$ for $\sigma\in\mathfrak{S}_{\Gamma}$ and $e\in\mathfrak{E}(\sigma)$.

The resulting finite volume scheme has the form: find $T=((T_{p})_{p\in\mathfrak{T}},(T_{\sigma})_{\sigma\in\mathfrak{S}_{\rm Dir}},(T_{e})_{e\in\mathfrak{E}_{\Gamma}})\in V_{h}$ such that:

[TABLE]

Remark 3* (Conservativity of the fluxes).*

Because they correspond to a cell-centred finite volume method, the inner fluxes $\mathcal{F}^{\Omega}_{p,\sigma}$ must satisfy by design the conservativity condition (18) on any vector $\varphi\in V_{h}$. On the contrary, the fluxes $\mathcal{F}^{\Gamma}_{\sigma,e}$ correspond to a cell- and edge-centred method and their conservativity is only imposed on the solution to the finite volume scheme (see Equation (19b)). See also [9, Sections 3.3.1 and 3.3.3] on this topic.

2.4 Error estimate

The following assumptions are made on the diffusive fluxes.

Assumption 4**.**

The numerical fluxes satisfy:

Coercivity: there is $\rho_{\Omega}>0$ and $\rho_{\Gamma}>0$ such that, for all $\varphi\in V_{h}$,

[TABLE]

and

[TABLE] 2. 2.

Consistency: there exist constants $C_{{\rm cons}}^{\Omega}$ and $C_{{\rm cons}}^{\Gamma}$ such that, for all $u\in C^{2}(\overline{\Omega})$ with $u=0$ on $\partial\Omega\backslash\Gamma$,

[TABLE]

and

[TABLE]

The error estimate will be established under the assumption that the following mesh regularity factor remains bounded above:

[TABLE]

where $d^{\bot}_{p,\sigma}$ is the orthogonal distance between $\vec{x}_{p}$ and $\sigma$.

Remark 5* (Interpretation of ${\rm reg}_{\mathfrak{T}}$).*

Bounding ${\rm reg}_{\mathfrak{T}}$ above imposes that each $\vec{x}_{p}$ must be “well within” its cell $p$, and that two neighbouring cells must have comparable diameters (which does not prevent local refinement, provided that it is done in layers of smoothly refined meshes).

Combining [13, Lemmas B.21 and B.31], we obtain the following discrete trace inequality: there is $C_{\rm tr}>0$ depending only on $\Omega$, $\Gamma$ and an upper bound of ${\rm reg}_{\mathfrak{T}}$ such that, for all $\varphi\in V_{h}$,

[TABLE]

In the rest of the paper, the notation $a\lesssim b$ means that $a\leq Cb$ for a constant $C$ that is independent of the quantities in $a$ and $b$, and of the mesh (but that may depend on $\Omega$, $\vec{W}$, $\Gamma$, $\rho_{\Omega}$, $\rho_{\Gamma}$, $C_{{\rm cons}}^{\Omega}$, $C_{{\rm cons}}^{\Gamma}$, $R$ and an upper bound of ${\rm reg}_{\mathfrak{T}}$). We can now state our main error estimate.

Theorem 6** (Error estimate).**

Under Assumption 4, suppose that $\vec{W}$ satisfies

[TABLE]

Assume that the solution $\overline{T}$ to (1) belongs to $C^{2}(\overline{\Omega})$, and let $T$ be the solution to the scheme (19). Then,

[TABLE]

Proof.

See Appendix A. ∎

Remark 7* (About Assumption (26)).*

Assumption (26a) imposes a relative smallness only of the positive part of $\nabla_{\!\Gamma}\!\boldsymbol{\cdot}\!\vec{W}$. In particular, this assumption holds if $\nabla_{\!\Gamma}\!\boldsymbol{\cdot}\!\vec{W}\leq 0$. Assumption (26b) shows how the user-defined parameter $R$ must be chosen to ensure the stability of the method.

Remark 8* (Regularity assumption on $\overline{T}$).*

In most situations, the $C^{2}$ regularity on $\overline{T}$ can be weakened to an $H^{2}$ regularity, upon additional technicalities that we do not address here to simplify the exposition. See, e.g., [13, Section 7.4] for lemmas useful for establishing consistency estimates under $H^{2}$-regularity of the function.

Remark 9* (Assumption (13)).*

In case Assumption (13) is not satisfied, that is the points $\vec{x}_{p}$ corresponding to cells that touch $\Gamma$ do not lie on $\Gamma$, the scheme has to be slightly modified the following way:

• •

Additional unknowns on the faces on $\Gamma$ are introduced, so that $V_{h}$ is changed into

[TABLE]

A point $\vec{x}_{\sigma}$ is chosen on each $\sigma\in\mathfrak{S}_{\Gamma}$ and the interpolant (15) is extended by setting, for these faces, $\varphi_{\sigma}=\varphi(\vec{x}_{\sigma})$.

• •

The seminorms $|{\cdot}|_{V_{h},\Omega}$ and $|{\cdot}|_{V_{h},\Gamma}$ are modified in the following way: in (14b) the sum is taken over $\sigma\in\mathfrak{S}$ with $\varphi_{q}=\varphi_{\sigma}$ whenever $\sigma\in\mathfrak{S}_{\Gamma}$; in (14c), $\varphi_{p}$ is replaced with $\varphi_{\sigma}$.

• •

Fluxes $\mathcal{F}^{\Omega}_{p,\sigma}$ are also considered for $\sigma\in\mathfrak{S}_{\Gamma}$ and the scheme consists in finding $T\in V_{h}$ solution to the conservativity equations (19b) and

[TABLE]

where, in (29), $p$ is the only cell that has $\sigma$ as face.

• •

The coercivity assumption (20) is changed into

[TABLE]

The analysis performed in Appendix A can then be adapted and leads to the same error estimate (27).

3 Numerical tests

The numerical tests presented here are obtained using internal fluxes $\mathcal{F}^{\Omega}_{p,\sigma}$ corresponding to a corrected version of the ones introduced in [29], and variants. For boundary fluxes $\mathcal{F}^{\Gamma}_{\sigma,e}$, used only for stabilisation purposes, we utilise the ones provided by the Hybrid Mimetic Mixed method [14].

3.1 Description of the scheme

3.1.1 Inner fluxes

We consider a structured, but not necessarily Cartesian, grid of points on $\Omega$. These points are called representative points, as this is where we will look for an approximation of the potential $\overline{T}$. The structured grid assumption means that the representative points can be denoted by $\vec{x}_{i,j,k}$, where $i\in\{0,\ldots,I+1\}$, $j\in\{0,\ldots,J+1\}$, $k\in\{0,\ldots,K+1\}$, and we assume that the extremal points (corresponding to $i=0$, $i=I+1$, $j=0$, $j=J+1$, $k=0$ or $k=K+1$) lie on $\partial\Omega$. We split the set of indices of these extremal points into $\mathcal{I}_{\Gamma}=\{(i,j,k)\,:\,k=0\}$ and its complement $\mathcal{I}_{D}$, and we assume that $\mathcal{I}_{\Gamma}$ corresponds to the points $\vec{x}_{i,j,k}\in\Gamma$, so that $\mathcal{I}_{D}$ is the set of indices for the points on the Dirichlet boundary $\partial\Omega\backslash\Gamma$. The points associated with two extremal indices lie on the edges of $\Omega$, whereas those with three extremal indices describe the corners of $\Omega$. Note that $\Omega$ is not necessarily a hexahedron since its “faces” may not be planar. We refer to Figs. 2 and 5 for illustrations.

For each $(i,j,k)\not\in\mathcal{I}_{D}$, a hexahedral finite volume is constructed around $\vec{x}_{i,j,k}$ using the following procedure. Note that points with indices in $\mathcal{I}_{D}$ are not associated with control volumes, as they lie on the Dirichlet boundary and they are therefore not associated with unknowns of the scheme.

• •

If $(i,j,k)\in\mathcal{I}_{\rm int}:=[2,I-1]\times[2,J-1]\times[2,K-1]$, then setting $A=\{(m,n,o)\in\{-1,0,1\}^{3}\,:\,|m|+|n|+|o|=3\}$ we define, for $(m,n,o)\in A$, the vertex $\vec{x}_{i,j,k}^{m,n,o}$ as an average of the eight neighbouring points in the grid, one of them being $\vec{x}_{i,j,k}$:

[TABLE]

where B(m,n,o)=\{(m,n,o),$$(m,n,0),$$(m,0,o),$$(m,0,0),$$(0,n,o),$$(0,n,0),$$(0,0,o),$$(0,0,0)\}. The finite volume around $\vec{x}_{i,j,k}$ is then the hexahedron (with possibly non-planar faces) defined by the vertices $\{\vec{x}_{i,j,k}^{m,n,o}\,:\,(m,n,o)\in A\}$. See Fig. 2 (left) for an illustration.

• •

If $(i,j,k)\in\mathcal{I}_{\Gamma}$, so that $k=0$, and $(i,j)\in[2,I-1]\times[2,J-1]$, we construct four vertices $\vec{x}_{i,j,k}^{m,n,1}$, for $(m,n,1)\in A$, as in (30). Four more vertices $\vec{x}_{i,j,k}^{m,n,0}$ are constructed by averaging the four neighbouring vertices on $\Gamma$:

[TABLE]

where $C(m,n)=\{(m,n),(m,0),(0,n),(0,0)\}$. The control volume associated with $\vec{x}_{i,j,0}$ is defined by the eight vertices thus constructed, and we notice that $\vec{x}_{i,j,0}$ lies on one of its faces (the one on $\Gamma$), so that (13) is satisfied.

• •

If $(i,j,k)\not\in(\mathcal{I}_{\rm int}\cup\mathcal{I}_{\Gamma})$, $\vec{x}_{i,j,k}$ is associated with a control volume touching the Dirichlet boundary and built from four vertices constructed as in (30) and four other vertices constructed in a similar way as in (31), using representative points on the Dirichlet boundary $\partial\Omega\backslash\Gamma$. See Fig. 2 (right).

• •

A similar construction is made for the remaining indices $(i,j,k)$, corresponding to control volumes with an edge along an edge of $\Omega$, or a vertex at one of the corners of $\Omega$; for example, the vertices of the control volumes lying on an edge of $\Omega$ are constructed as the average of two representative points $\vec{x}_{a,b,c}$ with two extremal indices. See Fig. 2 (right).

A generic finite volume is therefore identified by a triplet $(i,j,k)\not\in\mathcal{I}_{D}$. For simplicity and to relate more to the unstructured notations used in Section 2, we denote $(i,j,k)$ by $p$. The point $\vec{x}_{i,j,k}$ associated with $p$ is therefore denoted by $\vec{x}_{p}$. Any face $\sigma\in\mathfrak{S}(p)\backslash\mathfrak{S}_{\Gamma}$ can be associated with two representative points on each side: $\vec{x}_{p}$ itself, and $\vec{x}_{q}$ which might either be associated with a genuine control volume if $\sigma\in\mathfrak{S}_{\rm int}$, or with $\sigma$ itself (as a degenerate cell $q$) if $\sigma\in\mathfrak{S}_{\rm Dir}$. We write $\sigma=\sigma_{pq}$ and notice that $\sigma_{pq}$ may not be planar.

The four vertices $(\vec{x}_{i,j,k}^{m,n,o})_{(m,n,o)}$ of $\sigma_{pq}$ are ordered in a counterclockwise way, respective to the orientation compatible with the outer normal to $p$, and we denote them by $\vec{x}_{pq}^{\oplus}$, $\vec{x}_{pq}^{\boxplus}$, $\vec{x}_{pq}^{\ominus}$, $\vec{x}_{pq}^{\boxminus}$. For $\vec{x}_{pq}^{*}$ one of these vertices, we let $\mathcal{R}(\vec{x}_{pq}^{*})$ be the set of representative points involved in the construction of $\vec{x}_{pq}^{*}$; hence, if ${\rm Card}(\mathcal{R}(\vec{x}_{pq}^{*}))\in\{8,4,2,1\}$ is the cardinality of $\mathcal{R}(\vec{x}_{pq}^{*})$, we have

[TABLE]

We define four vectors related to the face $\sigma_{pq}$: the unit vector $\vec{s}_{pq}$ which points from $\vec{x}_{p}$ to $\vec{x}_{q}$

[TABLE]

two tangent vectors to the face

[TABLE]

and

[TABLE]

Due to the orientation chosen on $\sigma_{pq}$ and the ordering of the vertices of this face, if $\vec{n}_{pq}:\sigma_{pq}\to\mathbb{R}^{3}$ is the pointwise outer unit normal to $p$ on $\sigma_{pq}$ we have

[TABLE]

Since $(\vec{s}_{pq},\vec{t}_{pq}^{\ocircle},\vec{t}_{pq}^{\square})$ form a basis of $\mathbb{R}^{3}$, there are $\beta_{pq}>0$ and $(\alpha_{pq}^{\ocircle},\alpha_{pq}^{\square})\in\mathbb{R}^{2}$ such that

[TABLE]

The numerical fluxes are then given by: for $p\in\mathfrak{T}$, $\sigma\in\mathfrak{S}(p)\backslash\mathfrak{S}_{\Gamma}$, and $\varphi\in V_{h}$,

[TABLE]

where

• •

$d_{pq}=|\vec{x}_{p}-\vec{x}_{q}|$ (as in Section 2), $d_{pq}^{\ocircle}=|\vec{x}_{pq}^{\oplus}-\vec{x}_{pq}^{\ominus}|$ and $d_{pq}^{\boxempty}=|\vec{x}_{pq}^{\boxplus}-\vec{x}_{pq}^{\boxminus}|$, and

• •

for $*\in\{\oplus,\ominus,\boxplus,\boxminus\}$, each point $\vec{y}\in\mathcal{R}(\vec{x}_{pq}^{*})$ is associated with a genuine or degenerate cell $r$ (possibly an edge or corner of $\Omega$); we then let $\varphi_{\vec{y}}=\varphi_{r}$ (with $\varphi_{r}=0$ if $r$ is an edge or corner on $\partial\Omega$) and, following (32), the secondary unknown $\varphi_{pq}^{*}$ located at the vertex $\vec{x}_{pq}^{*}$ is defined by

[TABLE]

Remark 10* (Correction of the flux in [29]).*

In [29], a similar flux is defined with right-hand side multiplied by $\frac{|\widetilde{\vec{n}}_{pq}|}{|\sigma_{pq}|}$ in (36) – that is, $\beta_{pq}$ is multiplied by $\frac{|\sigma_{pq}|}{|\widetilde{\vec{n}}_{pq}|}$. The consistency analysis in Section B.2.2 shows that this choice of flux is only consistent if the faces are asymptotically flat (that is, $\frac{|\widetilde{\vec{n}}_{pq}|}{|\sigma_{pq}|}\to 1$ as the mesh size tends to zero). The flux we define above is consistent even if some faces remain non-flat as the mesh size tends to zero.

3.1.2 Surface fluxes on $\Gamma$

We use the fluxes of the Mixed Finite Volumes [12], which is the finite volume presentation of the Hybrid Mimetic Mixed method (see [14] and [13, Section 13.2.2]). Let $\varphi\in V_{h}$ and define, for $\sigma\in\mathfrak{S}_{\Gamma}$ and denoting as usual by $p\in\mathfrak{T}$ the unique control volume that contains $\sigma$ in its boundary,

[TABLE]

where $\overline{\vec{x}}_{e}=\frac{1}{|e|}\int_{e}\vec{x}$ is the centre of mass of $e$. Assuming that $\vec{n}_{\sigma,e}$ is constant along $e$ (but see Remark 11 below), the Stokes formula and the definition (15) of $I_{h}$ easily show that $\nabla_{\Gamma,\sigma}$ is a consistent approximation of the tangential gradient on $\Gamma$ in the sense that, if $\varphi\in C^{2}(\Gamma)$,

[TABLE]

As a consequence, $S_{p,e}$ can be seen as the remainder of a discrete first order Taylor expansion. The HMM fluxes are then defined by: for all $\varphi\in V_{h}$ and all $\sigma\in\mathfrak{S}_{\Gamma}$, the family $(\mathcal{F}^{\Gamma}_{\sigma,e}(\varphi))_{e\in\mathfrak{E}(\sigma)}$ is the unique solution to

[TABLE]

where we recall that $d_{pe}^{\bot}$ is the orthogonal distance between $\vec{x}_{p}$ and $e$.

Remark 11* (Curved edges).*

The definition (38) is consistent if the only curvature of the edges is due to the curvature of $\Gamma$, that is, the unit normal vector $\vec{n}_{\sigma,e}$ to $\sigma$ on $e$ along $\Gamma$ is constant. However, the HMM remains asymptotically consistent on faces with slightly curved edges [6], in the sense that

[TABLE]

3.1.3 Properties of the fluxes

In this section, we show that, upon some mesh regularity assumption (that can be checked in practice during implementation), the inner and surface fluxes described in Sections 3.1.1 and 3.1.2 are coercive and consistent. As a consequence, Theorem 6 applies to the numerical scheme (19) based on these fluxes, and the error estimate (27) holds for this scheme.

We first define three mesh regularity factors. The first two are required to establish the properties of the inner fluxes (see Appendix), whereas the third one is linked to the properties of the HMM fluxes

The first regularity factor is related to the faces not lying on $\Gamma$:

[TABLE]

where $\vec{d}_{pq}^{*}=(|\vec{x}_{r}-\vec{x}_{pq}^{*}|)_{r\in\mathcal{R}(\vec{x}_{pq}^{*})}$ and $d_{pq}^{*}=d_{pq}^{\ocircle}$ if $*\in\{\ominus,\oplus\}$, $d_{pq}^{*}=d_{pq}^{\boxempty}$ if $*\in\{\boxminus,\boxplus\}$. 2. 2.

The definition of the second regularity factor requires the introduction of a few notations associated to a pair $(\vec{x}_{pq}^{\star},r)$, where $\vec{x}_{pq}^{*}$ is a vertex of a face $\sigma_{pq}\in\mathfrak{S}\backslash\mathfrak{S}_{\Gamma}$ and $r=p$ or $q$. We refer to Figure 3 for an illustration of these notations.

• •

If $\vec{x}_{pq}^{*}\in\Omega$ is an internal vertex, we let $F_{r,pq}^{*}$ be the set of the two control volumes that are neighbours of $r$, have $\vec{x}_{pq}^{*}$ as vertex, but are neither $p$ or $q$. The two control volumes in $F_{r,pq}^{*}$ have two neighbours in common: $r$ itself, and another control volume that we denote by $e_{r,pq}^{*}$.

• •

If $\vec{x}_{pq}^{*}\in\Gamma$, we let $F_{r,pq}^{*}$ be the set made of the only control volume neighbour of $r$, that has $\vec{x}_{pq}^{*}$ as vertex, but that is not $p$ or $q$.

• •

If $\vec{x}_{pq}^{*}$ lies on the Dirichlet boundary $\partial\Omega\backslash\Gamma$, it is the vertex of a face $\sigma\in\mathfrak{S}(r)\cap\mathfrak{S}_{\rm Dir}$. We let $F_{r,pq}^{*}$ be the set made of this face, which is identified to the degenerate control volume $q$.

We then need to know, for a given face $\sigma_{ab}\in\mathfrak{S}\backslash\mathfrak{S}_{\Gamma}$, for which triplet $(p,q,*)$ we have, for some $r=p$ or $q$ and $f\in F^{*}_{r,pq}$, $\{e^{*}_{r,pq},f\}=\{a,b\}$ or $\{r,f\}=\{a,b\}$. These triplets are described by the following two sets

[TABLE]

The second regularity factor is then $\varrho_{\mathfrak{T},\Omega}$, assumed to be $>0$, such that

[TABLE]

where

[TABLE]

and

[TABLE] 3. 3.

The third regularity factor is

[TABLE]

Proposition 12** (Properties of the fluxes).**

The fluxes defined in Sections 3.1.1 and 3.1.2 satisfy Assumption 4 with $\varrho_{\Omega}=\varrho_{\mathfrak{T},\Omega}$, $C_{{\rm cons}}^{\Omega}$ depending only on an upper bound of ${\rm reg}_{\mathfrak{T}}+{\rm reg}_{\mathfrak{T},\Omega}$, and $\varrho_{\Gamma}>0$ and $C_{{\rm cons}}^{\Gamma}$ depending only on an upper bound of ${\rm reg}_{\mathfrak{T},\Gamma}$.

Proof.

See Appendix B. ∎

Remark 13* (About the regularity factors).*

Bounding ${\rm reg}_{\mathfrak{T},\Omega}$ above imposes the proximity of a vertex $\vec{x}_{pq}^{*}$ and the representative points $\mathcal{R}(\vec{x}_{pq}^{*})$ involved in its definition, as well as the non-degeneracy of the faces (whose diagonals $\vec{t}_{pq}^{\ocircle}$ and $\vec{t}_{pq}^{\boxempty}$ must have a minimal angle) and the transversality of the vector $\vec{x}_{p}\vec{x}_{q}$ and the face $\sigma_{pq}$. All these properties are natural given our construction of the control volumes.

The regularity factor ${\rm reg}_{\mathfrak{T},\Gamma}$ plays the same role, for the mesh $\mathfrak{S}_{\Gamma}$ of $\Gamma$, as the regularity factor ${\rm reg}_{\mathfrak{T}}$ for the mesh $\mathfrak{T}$ of $\Omega$. See Remark 5 for an interpretation of these terms.

Bounding $\varrho_{\mathfrak{T},\Omega}$ below imposes that faces that share a common vertex must have comparable measures and diagonal lengths (the terms $\frac{|\sigma_{pq}|}{|\sigma_{ab}|}$ and $\frac{d_{ab}}{d_{pq}^{\diamondsuit}}$ remain bounded), and that $\vec{s}_{pq}$ is “not too far” from the orthogonal direction to $\sigma_{pq}$ (so that, recalling (35), $\beta_{pq}$ remains close to $1$ while $\alpha_{pq}^{\ocircle}$ and $\alpha_{pq}^{\boxempty}$ remain small compared to $1$).

All these regularity factors, as well as ${\rm reg}_{\mathfrak{T}}$, are easy to numerically evaluate for a given mesh during the implementation. If, as the mesh is refined, these computed factors remain bounded above (for ${\rm reg}_{\mathfrak{T}}$, ${\rm reg}_{\mathfrak{T},\Omega}$ and ${\rm reg}_{\mathfrak{T},\Gamma}$) or below (for $\varrho_{\mathfrak{T},\Omega}$), then it ensures the robustness and accuracy of the numerical output since the error estimate (27) then holds. Note however that these conditions on the regularity factors are merely sufficient, not necessary; the scheme can still, in some cases, converge even if these factors do not remain properly bounded.

3.2 Alternative schemes

In the numerical tests, we will also present the results using two alternative schemes to the ones described above. The first alternative scheme is similar to (16) in a way that it approximate the oblique derivative as an advection equation on the boundary. It uses an upwind discretisation, instead of a numerically stabilized centred discretisation, for the convective term on the boundary. The second alternative approach is similar to the approximation of inner fluxes (36). It approximates the fluxes through a boundary face on $\Gamma$ by splitting the normal derivative into an oblique component, in the direction $\vec{V}$, and a tangential component to $\Gamma$. Similar splitting, but just on uniform rectangular, radial or spherical grids, was presented in [26, 27], where, however, additional points outside domain for treatment of normal derivative were introduced which is made possible with uniform structured grids.

Upwind scheme

The boundary advection term $\left[\overline{T}\vec{W}\cdot\vec{n}\right]_{\sigma,e}$ in (16) is here discretised using an upwind approach (which, contrary to (19), does not require the introduction of numerical diffusion for stabilisation). The resulting scheme has the form

[TABLE]

where the boundary advective numerical flux $\mathcal{F}^{\Gamma,{\rm\scriptsize adv}}_{\sigma,e}(T)$, which approximates $\left[\,\overline{T}\vec{W}\cdot\vec{n}\right]_{\sigma,e}$, is given by

[TABLE]

Remark 14* (Theoretical analysis).*

The theoretical analysis of this scheme can be conducted in a similar way as the scheme in Section 3.1, but leads to an $\mathcal{O}(\sqrt{h})$ theoretical convergence rate (see in particular Remark 16).

Splitting scheme

Here, the oblique derivative is not recast as a normal component and a boundary advective component. Instead, it is directly used together with a tangential approximation to reconstruct the normal fluxes. The resulting scheme has the form

[TABLE]

where the numerical normal flux $\mathcal{F}^{\Gamma}_{p,\sigma}(T)$, that approximates $\nabla\overline{T}\cdot\vec{n}$, is given by

[TABLE]

Here, the coefficients $\beta_{\sigma}$, $\alpha_{\sigma}^{\square}$ and $\alpha_{\sigma}^{\ocircle}$ are given by (35) with $\vec{s}_{pq}=\vec{V}(\vec{x}_{p})$. They therefore correspond to the decomposition of $\widetilde{\vec{n}}_{pq}$ on the basis $(\vec{V}(\vec{x}_{p}),\boldsymbol{t}_{\sigma}^{\square},\boldsymbol{t}_{\sigma}^{\ocircle})$. The equation (50) approximates $\nabla\overline{T}\cdot\vec{n}$ using the oblique derivative $\nabla\overline{T}\cdot\vec{V}=g$ (see (1b)) and a tangential component, reconstructed from the boundary values of $T$ using the same principles as in Section 3.1.1.

Remark 15* (Theoretical analysis).*

Given the close proximity of the approximation (50) with the discretisation (36), the ideas developed in Appendix B could be adapted to yield an error estimate for this scheme. However, when establishing the coercivity of the method, additional boundary terms would have to be accounted for in the regularity factor (44). The scale of these additional negative terms is proportional to how tangential the oblique vector is, and the method fails to be coercive for problems with an oblique field that is too tangential to the boundary of the domain.

3.3 Results

To test the proposed methods we present three sets of numerical experiments. For all experiments, the exact solution is chosen to be $\overline{T}(\boldsymbol{x})=\frac{1}{\boldsymbol{x}-\boldsymbol{x}_{0}}$, where $\boldsymbol{x}_{0}=(-0.3,-0.2,-0.1)$. The regularity parameters (24), (42), (44), and (47) and the Experimental Order of Convergence (EOC) are presented.

3.3.1 Cubic domain and non-uniform mesh

The computational domain $\Omega$ for the first set of experiments is a cube with unit edge length. The boundary $\Gamma$, on which the oblique boundary condition is prescribed, corresponds to the bottom face of the cube. The mesh is a non-uniform one obtained constructing first a uniform grid with distance between representative points equal to $h_{u}$, and then moving each point by a random vector $\boldsymbol{r}$ with components in $(-0.15h_{u},0.15h_{u})$. Points on $\partial\Omega$ are only moved in a direction tangential to the boundary. Experiments with different oblique vector fields are presented, and the regularity parameters are presented in Table 1. We notice that all parameters remain in a range that makes Theorem 6 and Proposition 12 applicable.

The first experiment, whose results are presented in Table 2, shows the convergence of the method for a constant vector field $\vec{V}(\vec{x})=(-1,-1,-1)$. The method (19) displays a first order convergence in $L^{2}$ and energy norms, which confirms the theoretical prediction of Theorem 6 and Proposition 12. The absence of super-convergence in $L^{2}$ norm is not surprising, as specific Finite Volume methods are only known to super-converge under certain geometric conditions, and to fail to super-converge in some cases [15]. The rates of convergence for the upwind method (48) are around $1$ in $L^{2}$ norm but tend to $1/2$ in $V_{h,\Omega}$ norm, which is expected (see Remark 14). The splitting method (49) shows the best convergence rates: above second order in $L^{2}$ norm, and above first order in $V_{h,\Omega}$.

The second experiment considers the non-constant vector field $\vec{V}(x,y,z)=(x,y,-1)$. In this case the surface divergence of $\vec{W}(x,y,z)=(x,y,0)$ (see (2)) is $\nabla_{\!\Gamma}\!\boldsymbol{\cdot}\!\vec{W}(\vec{x})=2$. The tests show similar orders of convergence, albeit slightly reduced, as in the experiment with a constant vector field; see Table 3. The slight degradation could stem from the fact that the assumption (26a) is not fully satisfied on these meshes and with this vector field, or that the asymptotic rate has not been achieved at these mesh sizes.

In the third experiment on the cube we consider a divergence free rotational vector field $\vec{V}(x,y,z)=(-x,z,-1)$. The results in Table 4 show that, here again, the schemes behave in a similar way as with the other two vector fields.

3.3.2 Tesseroid domain with non-planar $\Gamma$

The next experiments are run on a computational domain with a non-planar boundary $\Gamma$. The discrete computational domain then does not exactly match $\Omega$, and vertices of the boundary faces do not have to lie on the boundary $\Gamma$. Moreover, in this construction, the tangent space to $\Gamma$ is not well defined everywhere so the co-normal $\vec{n}_{\sigma,e}$ is not well defined either in the Eq. (19). In this case we approximate the normal vector $\vec{n}_{\sigma,e}$ by the normalised version of $\left(\frac{\vec{N}_{p}+\vec{N}_{q}}{2}\right)\times\vec{e}$, where the vector $\vec{N}_{p}$ is a normal to the face $\sigma$, the vector $\vec{N}_{q}$ is normal to the neighbouring face on the other side of $e$, and the vector $\vec{e}$ is a tangent vector to the edge $e$, chosen such that $\vec{n}_{\sigma,e}$ is an outward normal to $\sigma$.

The experiments are performed on a non-uniform mesh of the tesseroid

[TABLE]

See Fig. 5 for an illustration. The oblique boundary condition is prescribed on the non-planar face corresponding to $r=1$:

[TABLE]

The regularity parameters of the considered meshes are presented in Table 5. As can be seen there, the regularity factor $\varrho_{\mathfrak{T},\Omega}$ seems to degenerate as the mesh size is reduced, indicating that the condition that ensure the coercivity of the inner fluxes (see Proposition 12) might not hold – which does not necessarily mean that the scheme actually fails to be coercive or to converge, since this is only a sufficient condition.

The tests present the convergence of the methods for a non-constant vector field $\boldsymbol{V}(\boldsymbol{x})=(0.3,0.2,0.1)-\boldsymbol{x}$. The results presented in Table 6 are similar to the ones obtained in Section 3.3.1, with perhaps slightly better rates of convergence across the board. In any case, the apparent decay of the regularity factor $\varrho_{\mathfrak{T},\Omega}$ does not seem to negatively impact the convergence of the schemes.

3.3.3 Spherical section domain with perturbed bottom $\Gamma$

This series of experiments is performed on a section of a spherical domain, with a perturbed bottom boundary $\Gamma$ (see Fig. 6):

[TABLE]

The regularity parameters for the considered meshes are presented in Table 7. The coercivity constant $\varrho_{\mathfrak{T},\Omega}$ is worse as in the tesseroid case, as it becomes negative. However, once again, since a lower bound on this constant is only a sufficient condition for the theoretical analysis, this does not mean that the schemes fail to converge, as the numerical results will show. We take the non-constant vector field $\boldsymbol{V}(\boldsymbol{x})=(0.3,0.2,0.1)-\boldsymbol{x}$. Table 8 shows that all three schemes behave in a similar way as in the previous tests of Sections 3.3.1 and 3.3.2. This indicates that our coercivity analysis (based on $\varrho_{\mathfrak{T},\Omega}$) is actually a bit too conservative regarding the robustness range of the discretisations.

3.3.4 Cubic domain, almost tangential vector field $\vec{V}$

In this final series of numerical experiments with an analytical solution, we show the advantage of the proposed scheme (19) and of the upwind scheme (48) over the splitting scheme (49). The computational domain $\Omega$ is the unit cube, with $\Gamma$ being its bottom. We take the vector field $\vec{V}=(11.4301,0,-1)$, which corresponds to the outer normal on $\Gamma$ rotated by $85\^{A}\textdegree$; this vector field is therefore almost tangential to the boundary. The tests are run on uniform meshes.

Table 9 presents the EOC for the proposed scheme and upwind scheme. The rates are sometimes degraded compared to the previous tests, but there is a clear convergence.

For the splitting scheme (49), the fact that $\vec{V}$ is almost tangential leads to very large values of $\frac{\alpha_{\sigma}^{*}}{\beta_{\sigma}}$ in (50). As a consequence, the negative coefficients in the coercivity factor are too large to be controlled by the positive coefficients; the scheme really becomes non-coercive and unstable, and the BiCGStab algorithm used to solve the system fails. This breakdown of a numerical method is probably the worst situation that one wants to avoid in practice, which indicates that in severely oblique situations the proposed new methods (19) and (48) should be preferred, despite yielding sometimes reduced rates of convergence.

3.4 Local gravity field modelling

In this section we present local gravity field modelling over Slovakia using terrestrial gravity data. The goal of this experiment is to compute a disturbing potential using presented FVM schemes with oblique BC from terrestrial measurements and Dirichlet BCs obtained from satellite based model. Then we transform obtained potential to quasi-geoidal heights and compare them with real measurements. On the upper and side boundaries, the GO_CONS_GCF_2_DIR_R5 model [7] was used and on the bottom boundary we used the surface gravity disturbances obtained from the available regular grid of gravity anomalies, with the resolution $20^{\prime\prime}\times 30^{\prime\prime}$, that was compiled from original gravimetric measurements [19]. The gravity anomalies were transformed into the gravity disturbances by official digital vertical reference model DVRM (www.geoportal.sk).

The domain was bounded by $\langle 16{}^{\circ},23{}^{\circ}\rangle$ meridians and $\langle 47{}^{\circ},50.5{}^{\circ}\rangle$ parallels. The side boundaries were chosen sufficiently far from the area of Slovakia in order to mitigate an influence of the prescribed Dirichlet BC generated from the satellite-only geopotential model. For more details about this influence see [16]. The heights were interpolated from SRTM30 PLUS model and the upper boundary is in the height of 240km above the reference ellipsoid.

Three experiments with the grid density $841\times 631\times 301$ were performed using the FVM schemes (19), (48) and (49). The accuracy of the simulations was tested using GNSS-leveling. From the available dataset of 61 GNSS-leveling benchmarks, three evident outliers were removed. Hence, we tested the obtained local quasi-geoid model at 58 points. The results are summarised in Table 10 and, for the method (19), they are visualized in Fig. 7. We see a comparable precision of all the methods in this experiment. With this grid resolution the standard deviation of residuals between numerical results and measurements for all schemes is around 5cm.

4 Conclusion

We developed a framework for designing and analysing Finite Volume schemes for the Laplace equation with oblique boundary conditions. This framework, which can easily be extended to more general second order differential equations, consists in splitting the boundary condition into a normal and a tangential component, the later being handled as an advection term along the boundary of the domain; to ensure optimal convergence rates, this advection term is discretised using a centered scheme, with added numerical diffusion for stability purposes. The convergence analysis was carried out under usual coercivity and consistency assumptions on the numerical fluxes, and therefore applies to a range of possible FV discretisations. This analysis establishes first-order rates of convergence in a discrete energy ($H^{1}$) norm.

We then constructed specific fluxes, in the case where the computational domain is discretised using generalised hexahedra, and we identified geometrical conditions, easy to check during simulations, that ensure their coercivity and consistency. Two alternative discretisations of the oblique boundary conditions were also presented: the first one uses an upwind FV discretisation of the boundary advection, the second is not based on a FV discretisation on the boundary, but rather on splitting the outer normal to the boundary into its oblique component, and a tangential component discretised using finite differences and the specific geometry of the mesh.

We then provided extensive numerical tests, designed to assess the accuracy and robustness of the method, for various choices of the computational domain, and of the oblique vector field defining the boundary conditions. These tests confirmed, for all three schemes, the theoretical first-order rate of convergence in energy norm. In some tests, the energy rate of convergence is actually apparently higher than the theoretical one (but the asymptotic convergence rate might not have been attained at the considered mesh sizes). The second variant, based on a splitting of the outer normal, seems to present the best accuracy in our initial tests, when the velocity field is not too tangential to the boundary. For a nearly tangential velocity field, this splitting scheme breaks down as the numerical solver fails to find a solution to it. On the contrary, the other two variants remain robust and convergent in this extreme situation, albeit with a reduced accuracy. All three schemes were used to compute a quasi-geoidal height in the region of Slovakia. For this test, all methods give results with comparable quality.

Appendix A Proof of Theorem 6

The proof hinges on the 3rd Strang lemma of [9]. Let us first recast the scheme (19) under a variational formulation. Take $\varphi\in V_{h}$, multiply (19a) by $\varphi_{p}$ and summing over $p\in\mathfrak{T}$ to get

[TABLE]

Using the conservativity of the fluxes $\mathcal{F}^{\Omega}_{p,\sigma}$ (see (18)), $\mathcal{F}^{\Gamma}_{\sigma,e}(T)$ (see (19b)) and $T_{e}\left[\vec{W}\cdot\vec{n}\right]_{\sigma,e}$ (by definition of $\left[\vec{W}\cdot\vec{n}\right]_{\sigma,e}$), and the zero value of $T_{e}$ if $e$ is a boundary edge in $\Gamma$, we gather the sums in the left-hand side by faces and edges as in [9, Proofs of Theorem 27 and 33] to find

[TABLE]

The solution $T\in V_{h}$ to the scheme thus satisfies $a_{h}(T,\varphi)=\ell_{h}(\varphi)$ for all $\varphi\in V_{h}$, with $a_{h}(T,\varphi)$ (resp. $\ell_{h}$) the bilinear form (resp. linear form) in the left-hand side (resp. right-hand side) of (51). Owing to the 3rd Strang lemma [9, Theorem 10], the estimate (27) follows if we establish the coercivity and consistency properties:

[TABLE]

and, letting $\mathcal{E}_{h}(\overline{T};\varphi):=\ell_{h}(\varphi)-a_{h}(\overline{T},\varphi)$ be the consistency error,

[TABLE]

A.1 Coercivity

The coercivity properties (20) and (21) show that

[TABLE]

Simple algebraic identities show that

[TABLE]

By conservativity of $\left[\vec{W}\cdot\vec{n}\right]_{\sigma,e}$ and zero value of $\varphi_{e}$ on boundary edges of $\Gamma$,

[TABLE]

Hence, since $\sum_{e\in\mathfrak{E}(\sigma)}\left[\vec{W}\cdot\vec{n}\right]_{\sigma,e}=\left[\nabla_{\!\Gamma}\!\boldsymbol{\cdot}\!\vec{W}\right]_{\sigma}$,

[TABLE]

Using $\left[\nabla_{\!\Gamma}\!\boldsymbol{\cdot}\!\vec{W}\right]_{\sigma}\leq\|(\nabla_{\!\Gamma}\!\boldsymbol{\cdot}\!\vec{W})^{+}\|_{C(\Gamma)}|\sigma|$ and the trace inequality (25), we write

[TABLE]

Plugging this into (55) and noticing, since $d_{pe}^{\bot}\leq h_{\Gamma}$, that

[TABLE]

we obtain

[TABLE]

Coming back to (54), we infer that

[TABLE]

Owing to Assumption 26, this proves (52).

A.2 Consistency

Using (16) and recalling that $\ell_{h}(\varphi)$ is defined by the right-hand side of (51), we write

[TABLE]

where we have used the conservativity of the fluxes to gather the sums by faces in the second equality. Subtracting $a_{h}(I_{h}\overline{T},\varphi)$ (given by the left-hand side of (51) with $T$ replaced by $I_{h}\overline{T}$), we can split the consistency error into four terms:

[TABLE]

with, setting $I_{h}\overline{T}=((\overline{T}_{p})_{p\in\mathfrak{T}},(\overline{T}_{\sigma})_{\sigma\in\mathfrak{S}_{\rm Dir}},(\overline{T}_{e})_{e\in\mathfrak{E}_{\Gamma}})$,

[TABLE]

We now estimate each of these terms.

Term $\mathcal{T}_{c,1}$.

Introducing $\sqrt{d_{pq}/|\sigma|}$ and using a Cauchy–Schwarz inequality, the consistency property (22) yields

[TABLE]

Let $D_{p\sigma}$ be the convex hull of $\vec{x}_{p}$ and $\sigma$. If $\sigma$ is flat, by [13, Lemma B.2] we have $|D_{p\sigma}|=d^{\bot}_{p,\sigma}|\sigma|/3$ and thus, by definition of ${\rm reg}_{\mathfrak{T}}$ (which implies $d_{pq}\leq{\rm diam}(p)+{\rm diam}(q)\lesssim{\rm diam}(p)\lesssim d_{p\sigma}^{\bot}$),

[TABLE]

This final estimate also holds in case of non-flat $\sigma$, as can be seen approximating $\sigma$ by piecewise flat surfaces. Hence,

[TABLE]

Term $\mathcal{T}_{c,2}$.

We first estimate the consistency of the fluxes involved in this term. Using the definition of the interpolant (15) we have $\int_{e}(\overline{T}-\overline{T}_{e})=0$ and thus

[TABLE]

where the conclusion follows from a mean value theorem on $\overline{T}$ and $\vec{W}$. Applying a Cauchy–Schwarz inequality and using (61), we infer

[TABLE]

In a similar way as in the last equalities in (59), $D_{pe}$ being the convex hull of $\vec{x}_{p}$ and $e$ we have

[TABLE]

Since $|\varphi|_{V_{h},\Gamma}\leq h_{\Gamma}^{-\frac{1}{2}}\|\varphi\|_{V_{h}}$, we conclude that

[TABLE]

Remark 16* (Centred discretisation of the advective term).*

The approximation in the first term of (17) corresponds to a centred discretisation of the advection term $\nabla_{\!\Gamma}\!\boldsymbol{\cdot}\!(\overline{T}\vec{W})$ on $\Gamma$. To stabilise this centred discretisation and ensure the coercivity of the scheme, we have to add the artificial diffusion through the terms $Rh_{\Gamma}\mathcal{F}^{\Gamma}(T)$ in (19a). A standard option to avoid adding numerical diffusion is to directly use an upwind discretisation of the advective term, as in (48). In this case, since stability would not require to introduce artificial diffusion, we would only consider cell unknowns in $V_{h}$ (and not introduce edge unknowns), and we would take $\|\cdot\|_{V_{h}}=|\cdot|_{V_{h},\Omega}$. The resulting scheme would be (48). Such a choice, however, would prevent us from introducing $\left[\vec{W}\cdot\vec{n}\right]_{\sigma,e}$ in (61) and the resulting estimate would be in $\mathcal{O}(h_{\Gamma})$ instead of $\mathcal{O}(h_{\Gamma}^{2})$. Carrying on as in (62) but with $\|\cdot\|_{V_{h},\Omega}$ instead of the absent $|\cdot|_{V_{h},\Gamma}$, with the natural assumption that $|e|h_{\Gamma}\lesssim|\sigma|$, we would arrive at (64) with $h_{\Gamma}^{\frac{1}{2}}$ instead of $h_{\Gamma}^{\frac{3}{2}}$. The final consistency estimate, and thus error estimate, would then be in $\mathcal{O}(h^{\frac{1}{2}})$ instead of $\mathcal{O}(h)$.

Term $\mathcal{T}_{c,3}$.

Notice first that

[TABLE]

Hence, using a Cauchy–Schwarz inequality and the trace inequality (25),

[TABLE]

Term $\mathcal{T}_{c,4}$.

Introducing the exact surface fluxes $\overline{F}_{\sigma,e}(\overline{T})=-\int_{e}\nabla\overline{T}\cdot\vec{n}_{\sigma,e}$ on $\Gamma$, we write

[TABLE]

A Cauchy–Schwarz inequality, the consistency property (23), and (63) show that

[TABLE]

For $\mathcal{T}_{c,4,2}$, we use the conservativity of the fluxes $\overline{F}_{\sigma,e}(\overline{T})$, the fact that $\varphi_{e}=0$ for edges on the boundary of $\Gamma$, $\sum_{e\in\mathfrak{E}(\sigma)}\overline{F}_{\sigma,e}(\overline{T})=-\iint_{\sigma}\Delta_{\Gamma}\overline{T}$, and the trace inequality (25) to write

[TABLE]

Combined with (66) and recalling that $\mathcal{T}_{4,c}=Rh_{\Gamma}\mathcal{T}_{4,c,1}+Rh_{\Gamma}\mathcal{T}_{4,c,2}$ this shows that

[TABLE]

Gathering (60), (64), (65) and (67) in (58), we infer that (53) holds, which concludes the proof of Theorem 6.

Appendix B Proof of Proposition 12

B.1 Boundary fluxes

The coercivity of the HMM fluxes result from the construction of the method as a Gradient Discretisation Method, see [13, Chapter 13] – we note that this coercivity is purely algebraic, and not impacted by the curvature of the faces $\sigma\in\mathfrak{S}_{\Gamma}$ or of their edges. In the case of flat faces and edges, the consistency of the fluxes (41) is a consequence of (40), see [13, Chapter 13] or [9, Example 31].

Consider now a curved face $\sigma\in\mathfrak{S}_{\Gamma}$, and assume that all edges of $\sigma$ are “only curved along $\sigma$” in the sense that $\vec{n}_{\sigma,e}$ is constant over $e$ for all $e\in\mathfrak{E}(\sigma)$ (see Remark 11 otherwise). Because $\Gamma$ is a smooth surface, $\overline{\vec{x}}_{e}-\vec{x}_{p}$ is asymptotically close to the tangent plane to $\Gamma$ at any point of $\sigma$ (that is, the projection of $\overline{\vec{x}}_{e}-\vec{x}_{p}$ in any normal direction to $\sigma$ has length $\mathcal{O}(h_{\Gamma}^{2})$). Taylor expansions at any point of $\sigma$ and the consistency property (40) thus give a constant $C$ independent of the mesh such that, for $\varphi\in C^{2}(\Gamma)$,

[TABLE]

Using this estimate and (40), the arguments developed in [13, Chapter 13] can then easily be adapted and yield (23).

B.2 Inner fluxes

B.2.1 Coercivity

Without any loss of generality we can select the vertex labels $\vec{x}_{pq}^{\oplus}$ and $\vec{x}_{pq}^{\ominus}$ (resp. $\vec{x}_{pq}^{\boxplus}$ and $\vec{x}_{pq}^{\boxminus}$) such that $\alpha_{pq}^{\ocircle}\geq 0$ (resp. $\alpha_{pq}^{\square}\geq 0$). Recalling the definition (36) of the fluxes, we use the zero value on the Dirichlet boundary and the Young inequality $xy\geq-\frac{1}{2}x^{2}-\frac{1}{2}y^{2}$ to write

[TABLE]

In order to establish (20), we now need to find a lower bound of this quantity in terms of sums of $(\varphi_{a}-\varphi_{b})^{2}$ for $(a,b)$ pairs of neighbouring control volumes. The first stage is to recast $\varphi_{pq}^{\oplus}-\varphi_{pq}^{\ominus}$ and $\varphi_{pq}^{\boxplus}-\varphi_{pq}^{\boxminus}$ as combinations of differences of $\varphi$ on neighbouring control volumes. Without loss of generality, we consider $\varphi_{pq}^{\oplus}-\varphi_{pq}^{\ominus}$. We have to deal with three cases, depending if the corresponding vertices are both internal, if one lies on $\Gamma$, or if one lies on the Dirichlet boundary $\partial\Omega\backslash\Gamma$.

Case 1: Internal vertices. We assume here that $\vec{x}_{pq}^{\oplus}$ and $\vec{x}_{pq}^{\ominus}$ are both in $\Omega$. Let $\vec{x}_{pq}^{*}$ denote any one of these two vertices. Recalling the definitions in Section 3.1.3 (see also Figure 3) of $F_{r,pq}^{*}$ and $e^{*}_{r,pq}$, we see that the set $\cup_{r=p,q}(F_{r,pq}^{*}\cup\{r,e_{r,pq}^{*}\})$ is made of the eight control volumes around $\vec{x}_{pq}^{*}$ whose unknowns are involved in the definition (37) of $\varphi_{pq}^{*}$. Hence, we can decompose $\varphi_{pq}^{\oplus}-\varphi_{pq}^{\ominus}$ as

[TABLE]

Inside the sum in the right-hand side, each cell unknown appears with the coefficients represented in Figure B.2.1 (left). Our goal is to gather these terms together in order to write $\varphi_{pq}^{\oplus}-\varphi_{pq}^{\ominus}$ as a combination of terms $\varphi_{b}-\varphi_{a}$ with $a$ and $b$ neighbouring control volumes. This is done by splitting the coefficients in order to associate (parts of) each cell unknown with a neighbouring cell unknown, as in Figure B.2.1 (right).