# Determination of a high spatial resolution geopotential model using atomic clock comparisons

Guillaume Lion, Isabelle Panet, P. Wolf, C. Guerlin, S. Bize, P. Delva

► To cite this version:

Guillaume Lion, Isabelle Panet, P. Wolf, C. Guerlin, S. Bize, et al. Determination of a high spatial resolution geopotential model using atomic clock comparisons. Journal of Geodesy, 2017, 10.1007/s00190-016-0986-6. hal-01436820

HAL Id: hal-01436820

https://hal.sorbonne-universite.fr/hal-01436820

Submitted on 16 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ORIGINAL ARTICLE

# Determination of a high spatial resolution geopotential model using atomic clock comparisons

G. Lion<fnref n="1" />,<fnref n="2" /> ID · I. Panet<fnref n="2" /> · P. Wolf<fnref n="1" /> · C. Guerlin<fnref n="1" />,<fnref n="3" /> · S. Bize<fnref n="1" /> · P. Delva<fnref n="1" />

Received: 26 July 2016 / Accepted: 10 December 2016
© The Author(s) 2016. This article is published with open access at Springerlink.com

**Abstract** Recent technological advances in optical atomic clocks are opening new perspectives for the direct determination of geopotential differences between any two points at a centimeter-level accuracy in geoid height. However, so far detailed quantitative estimates of the possible improvement in geoid determination when adding such clock measurements to existing data are lacking. We present a first step in that direction with the aim and hope of triggering further work and efforts in this emerging field of chronometric geodesy and geophysics. We specifically focus on evaluating the contribution of this new kind of direct measurements in determining the geopotential at high spatial resolution (≈10 km). We studied two test areas, both located in France and corresponding to a middle (Massif Central) and high (Alps) mountainous terrain. These regions are interesting because the gravitational field strength varies greatly from place to place at high spatial resolution due to the complex topography. Our method consists in first generating a synthetic high-resolution geopotential map, then drawing synthetic measurement data (gravimetry and clock data) from it, and finally reconstructing the geopotential map from that data using least squares collocation. The quality of the recon-

structed map is then assessed by comparing it to the original one used to generate the data. We show that adding only a few clock data points (less than 1% of the gravimetry data) reduces the bias significantly and improves the standard deviation by a factor 3. The effect of the data coverage and data quality on the results is investigated, and the trade-off between the measurement noise level and the number of data points is discussed.

**Keywords** Chronometric geodesy · High spatial resolution · Geopotential · Gravity field · Atomic clock · Least squares collocation (LSC) · Stationary covariance function

## 1 Introduction

Chronometry is the science of the measurement of time. As the time flow of clocks depends on the surrounding gravity field through the relativistic gravitational redshift predicted by Einstein (Landau and Lifshitz 1975), chronometric geodesy considers the use of clocks to directly determine Earth's gravitational potential differences. Instead of using state-of-the-art Earth's gravitational field models to predict frequency shifts between distant clocks (Pavlis and Weiss (2003), ITOC project<fnref n="1" />), the principle is to reverse the problem and ask ourselves whether the comparison of frequency shifts between distant clocks can improve our knowledge of Earth's gravity and geoid (Bjerhammar 1985; Mai 2013; Petit et al. 2014; Shen et al. 2016; Kopeikin et al. 2016). For example, two clocks with an accuracy of 10⁻¹⁸ in terms of relative frequency shift would detect a 1-cm geoid height variation between them, corresponding to a geopotential variation ΔW

<figure><img src="image_3.png" /></figure>

G. Lion
Guillaume.Lion@obspm.fr

I. Panet
isabelle.panet@ensg.eu

<footnote n="1">¹ LNE-SYRTE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, 61 avenue de l Observatoire, 75014 Paris, France</footnote>

<footnote n="2">² LASTIG LAREG, IGN, ENSG, Univ Paris Diderot, Sorbonne Paris Cité, 35 rue Hélène Brion, 75013 Paris, France</footnote>

<footnote n="3">³ Laboratoire Kastler Brossel, ENS-PSL Research University, CNRS, UPMC-Sorbonne Universités, Collège de France, 24 rue Lhomond, 75005 Paris, France</footnote>

<footnote n="1">¹ http://projects.npl.co.uk/.</footnote>
of about $0.1 \text{ m}^2\text{s}^{-2}$ (for more details, see, e.g., Delva and Lodewyck 2013; Mai 2013; Petit et al. 2014).

Until recently, the performances of optical clocks had not been sufficient to make applications in practice for the determination of Earth's gravity potential. However, ongoing quick developments of optical clocks are opening these possibilities. Chou et al. (2010) demonstrated the ability of the new generation of atomic clocks, based on optical transitions, to sense geoid height differences with a 30-cm level of accuracy. To date, the best of these instruments reach a stability of $1.6 \times 10^{-18}$ (NIST, RIKEN + Univ. Tokyo, Hinkley et al. 2013) after 7 hours of integration time. More recently, an accuracy of $2.1 \times 10^{-18}$ (JILA, Nicholson et al. 2015) has been obtained, equivalent to geopotential differences of $0.2 \text{ m}^2\text{s}^{-2}$, or 2 cm on the geoid. Recently, Takano et al. (2016) demonstrated the feasibility of cm-level chronometric geodesy. By connecting clocks separated by 15 km with a long telecom fiber, they found that the height difference between the distant clocks determined by the chronometric leveling (see Vermeer 1983) was in agreement with the classical leveling measurement within the clocks uncertainty of 5 cm. Other related work using optical fiber or coaxial cable time-frequency transfer can be found in (Shen 2013; Shen and Shen 2015).

Such results stress the question of what can we learn about Earth's gravity and mass sources using clocks that we cannot easily derive from existing gravimetric data. Recent studies address this question; for example, Bondarescu et al. (2012) discussed the value and future applicability of chronometric geodesy for direct geoid mapping on continents and joint gravity potential surveying to determine subsurface density anomalies. They find that a geoid perturbation caused by a 1.5-km radius sphere with 20 percent density anomaly buried at 2 km depth in the Earth's crust is already detectable by atomic clocks with present-day accuracy. They also investigate other applications, for earthquake prediction and volcanic eruptions (Bondarescu et al. 2015b), or to monitor vertical surface motion changes due to magnetic, post-seismic, or tidal deformations (Bondarescu et al. 2015a, c).

Here we will consider the “static” or “long-term” component of Earth’s gravity. Our knowledge of Earth’s gravitational field is usually expressed through geopotential grids and models that integrate all available observations, globally or over an area of interest. These models are, however, not based on direct observations with the potential itself, which has to be reconstructed or extrapolated by integrating measurements of its derivatives. Yet, this quantity is needed in itself, like using a high-resolution geoid as a reference for height on land and dynamic topography over the oceans (Rummel and Teunissen 1988; Rummel 2002, 2012; Sansò and Venuti 2002; Zhang et al. 2008; Sansò and Sideris 2013; Marti 2015).

The potential is reconstructed with a centimetric accuracy at resolutions of the order of 100 km from GRACE and GOCE satellite data (Pail et al. 2011; Bruinsma et al. 2014) and integrated from near-surface gravimetry for the shorter spatial scales. As a result, the standard deviation (rms) of differences between geoid heights obtained from a global high-resolution model as EGM2008, and from a combination of GPS/leveling data, reaches up to 10 cm in areas well covered in surface data (Gruber 2009). The uneven distribution of surface gravity data, especially in transitional zones (coasts, borders between different countries) and with important gaps in areas difficult to access, indeed limits the accuracy of the reconstruction when aiming at a centimeter-level of precision. This is an important issue, as large gravity and geoid variations over a range of spatial scales are found in mountainous regions, and because a high accuracy on altitudes determination is crucial in coastal zones. Airborne gravity surveys are thus realized in such regions (Johnson 2009; Douch et al. 2015); local clock-based geopotential determination could be another way to overcome these limitations.

In this context, here, we investigate to what extent clocks could contribute to fill the gap between the satellite and near-surface gravity spectral and spatial coverages in order to improve our knowledge of the geopotential and gravity field at all wavelengths. By nature, potential data are smoother and more sensitive to mass sources at large scales than gravity data, which are strongly influenced by local effects. Thus, they could naturally complement existing networks in sparsely covered places and even also contribute to point out possible systematic patterns of errors in the less recent gravity data sets. We address the question through test case examples of high-resolution geopotential reconstructions in areas with different characteristics, leading to different variabilities of the gravity field. We consider the Massif Central in France, marked by smooth, moderate altitude mountains and volcanic plateaus, and an Alps–Mediterranean zone, comprising high reliefs and a land/sea transition.

Throughout this work, we will treat clock measurements as direct determinations of the disturbing potential $T$ (see below and Sect. 3 for details). We implicitly assume that the actual measurements are the potential differences between the clock location and some reference clock(s) within the area of interest. These measurements are obtained by comparing the two clocks over distances of up to a few 100 km. Currently two methods are available for such comparisons, fiber links (Lisdat et al. 2016) and free space optical links (Deschênes et al. 2016). The free space optical links are most promising for the applications considered here, but are presently still limited to short (few km) distances. However, projects for extending these methods based on airborne or satellite relays are on the way, but still require some effort in technology development.
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Scheme of the numerical approach used to evaluate the contribution of atomic clocks to determine the geopotential</figcaption></figure>

The paper is organized as follows. In Sect. 2, we briefly summarize the method schematically. In Sect. 3, we describe the regions of interest and the construction of the high-resolution synthetic data sets used in our tests. In Sect. 4, we present the methodology to assess the contribution of new clock data in the potential recovery, in addition to ground gravity measurements. Numerical results are shown in Sect. 5. We finally discuss in Sect. 6 the influence of different parameters like the data noise level and coverage.

## 2 Method

The rapid progress of optical clocks performances opens new perspectives for their use in geodesy and geophysics. While they were until recently built only as stationary laboratory devices, several transportable optical clocks are currently under construction or test (see, e.g., Bongs 2015; Origlia et al. 2016; Vogt et al. 2016). The technological step toward state-of-the-art transportable optical clocks is likely to take place within the next decade. In parallel, in order to assess the capabilities of this upcoming technology, we chose an approach based on numerical simulation in order to investigate whether atomic clocks can improve the determination of the geopotential. Based on the consideration that ground optical clocks are more sensitive to the longer wavelengths of the gravitational field around them than gravity data, our method is adapted to the determination of the geopotential at regional scales. In Fig. 1a scheme of the method used in this paper is shown:

1. In the first step, we generate a high spatial resolution grid of the gravity disturbance $\delta g$ and the disturbing potential $\tilde{T}$, considered as our reference solutions. This is done using a state-of-the-art geopotential model (EIGEN-6C4) and by removing low and high frequencies. It is described in details in Sect. 3.

2. In the second step, we generate synthetic measurements $\delta g$ and $\tilde{T}$ from a realistic spatial distribution, and then we add generated random noise representative of the measurement noise. This is described in details in Sect. 4.

3. In a third step, we estimate the disturbing potential $\tilde{T}$ from the synthetic measurements $\delta g$ and/or $\tilde{T}$ on a regular grid thanks to least square collocation (LSC) method. Interpolating spatial data are realized by making an assumption on the a priori gravity field regularity on the target area, as described in Sect. 5. This prior is expressed by the covariance function of the gravity potential and its derivatives. It allows to predict the disturbing potential on the output grid from the observations using the signal correlations between the data points and with the estimated potential.

4. Finally, we evaluate the potential recovery quality for different data distribution sets, noise levels, and types of data, by comparing the statistics of the residuals $\delta$ between the estimated values $\tilde{T}$ and the reference model $\tilde{T}$.

Let us underline that in this work, we use synthetic potential data while a network of clocks would give access to potential differences between the clocks. We indeed assume that the clocks-based potential differences have been connected to one or a few reference points, without introducing additional biases larger than the assumed clock uncertainties. Note that these reference points are absolute potential points determined by other methods (GNSS/geoid for example).

In this differential method, significant residuals $\delta$ (higher than the machine precision) can have several origins, depending on the parameters of the simulation that can be varied:

1. The modeled instrumental noise added to the reference model at step 2. This noise can be changed in order to determine, for instance, whether it is better to reduce gravimetry noise by one order of magnitude, rather than using clock measurements.
2. The data distribution chosen in step 2. This is useful to check for instance the effect of the number of clock measurements on the residuals or to find an optimal coverage for the clock measurements.

3. The potential estimation error, due to the intrinsic imperfection of the covariance model chosen for the geopotential. In our case, this is due to the low-frequency content of the covariance function chosen for the least square collocation method (see Sect. 5).

All these sources of errors are somewhat entangled with one another, such that a careful analysis must be done when varying the parameters of the simulation. This is discussed in details in Sect. 6.

# 3 Regions of interest and synthetic gravity field reference models

## 3.1 Gravity data and distribution

Our study focuses on two different areas in France. The first region is the Massif Central located between 43° to 47°N and 1° to 5°E and consists of plateaus and low mountain range, see Fig. 2. The second target area, much more hilly and mountainous, is the French Alps with a portion of the Mediterranean Sea located at the limit of different countries and bounded by 42° to 47°N and 4.5° to 9°E, see Fig. 3. Topography is obtained from the 30-m digital elevation model over France by IGN, completed with Smith and Sandwell (1997) bathymetry and SRTM data.

Available surface gravity data in these areas, from the BGI (International Gravimetric Bureau), are shown in Figs. 2b–3b. Note that the BGI gravity data values are not used in this study, but only their spatial distribution in order to generate realistic distribution in the synthetic tests. In these figures, it is shown that the gravity data are sparsely distributed: The plain is densely surveyed while the mountainous regions are poorly covered because they are mostly inaccessible by the conventional gravity survey. The range of free-air gravity anomalies (see Moritz 1980; Sansò and Sideris 2013) which are quite large reflects the complex structure of the gravity field in these regions, which means that the gravitational field strength varies greatly from place to place at high resolution. The scarcity of gravity data in the hilly regions is thus a major limitation in deriving accurate high-resolution geopotential model.

## 3.2 High-resolution synthetic data

Here, we present the way to simulate our synthetic gravity disturbances $\delta g$ and disturbing potentials $\bar{T}$ by subtracting

<figure><img src="image_5.png" /><figcaption>Fig. 2 Topography and gravity data distribution in the Alps–Mediterranean area. a Topography. b Terrestrial and marine free-air gravity anomalies</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3 Topography and gravity data distribution in the Massif Central area. a Topography. b Terrestrial and marine free-air gravity anomalies</figcaption></figure>

the gravity field long and short wavelengths influence of a high-resolution global geopotential model.

The generation of the synthetic data $\delta g$ and $\bar{T}$ at the Earth's topographic surface was carried out, in ellipsoidal approximation, with the FORTRAN program GEOPOT<fnref n="2" /> (Smith 1998) of the National Geodetic Survey (NGS). This program allows to compute gravity field-related quantities at given locations using a geopotential model and additional information such as parameters of the ellipsoidal normal field, tide system. The ellipsoidal normal field is defined by the parameters of the geodetic reference system GRS80 (Moritz 1984). As input, we used the static global gravity field model EIGEN-6C4 (Förste et al. 2014). It is a combined model up

<figure><img src="image_7.png" /></figure>

<footnote n="2">² http://www.ngs.noaa.gov/GEOID/RESEARCH Software/ research_software.html.</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 4** High-pass filter based on a Poisson wavelet $\Phi$ at order $m = 3$. The cutoff is $n_{\text{cut}} = 100$ and the wavelet scale is 0.03</figcaption></figure>

to degree and order (d/o) 2190 containing satellite, altimetry, terrestrial gravity, and elevation data. By using the spherical harmonics (SH) coefficients up to d/o 2000, it allows us to map gravity variations down to 10 km resolution. Thus, these synthetic data do not represent the full geoid signal. The choice is motivated by the fact that at a centimeter-level of accuracy, we expect large benefit from clocks at wavelengths $\ge 10$ km.

Our objective is to study how clocks can advance knowledge of the geoid beyond the resolution of the satellites. In a first step, as illustrated in Fig. 4, the long wavelengths of the gravity field covered by the satellites and longer than the extent of the local area are completely removed up to the degree $n_{\text{cut}} = 100$ (200 km resolution). This data reduction is necessary for the determination of the local covariance function in order to have centered data, or close to zero, as detailed in Knudsen (1987, 1988). Between degree 101 and 583, the gravity field is progressively filtered using 3 Poisson wavelets spectra (Holschneider et al. 2003), while its full content is preserved above degree 583. In this way, we realize a smooth transition between the wavelengths covered by the satellites and those constrained from the surface data.

To subtract the terrain effects included in EIGEN-6C4, we used the topographic potential model dV_ELL_RET2012 (Claessens and Hirt 2013) truncated at d/o 2000. Complete up to d/o 2160, this model provides in ellipsoidal approximation the gravitational attraction due to the topographic masses anywhere on the Earth's surface. The results of this data reduction yields to the reference fields $\delta g$ and $\bar{T}$ for both regions, shown in Figs. 5 and 6.

Figures 5 and 6 show the different characteristics of the residual field in these two regions. The residual anomalies have smaller amplitudes in the Massif Central area when compared to the Alps. In addition, the presence of high moun-

<figure><img src="image_9.png" /><figcaption>**Fig. 5** Synthetic reference fields of gravity disturbances $\delta g$ and disturbing potential $\bar{T}$ in the Massif Central area. Anomalies are computed at the Earth's topographic surface from the EIGEN-6C4 model up to d/o 2000 after removal of the low and high frequencies of the gravity field</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 6** Synthetic reference fields of gravity disturbances $\delta g$ and disturbing potential $\bar{T}$ in the Alps–Mediterranean area. Anomalies are computed at the Earth’s topographic surface from the EIGEN-6C4 model up to d/o 2000 after removal of the low and high frequencies of the gravity field</figcaption></figure>

tains on part of the latter zone results in an important spatial heterogeneity of the residual gravity anomalies, with large signals also at intermediate resolutions.

## 4 Data set selection and synthetic noise

Gravimetric location points selection. Our goal is to reproduce a realistic spatial distribution of the gravity points. The BGI gravity data sets contain hundreds of thousands points for the target regions (see Figs. 2b–3b). In order to reduce the size of the problem and make it numerically more tractable,
<figure><img src="image_11.png" /><figcaption>**Fig. 7** Distribution of the gravity and clock data used in the synthetic tests. **a** Massif Central: 4374 gravity data and 33 potential data, **b** Alps: 4959 gravity data and 32 potential data</figcaption></figure>

we build a distribution with no more than several thousand
points from the original one.

Starting from the spatial distribution of the BGI gravity
data sets, a grid $\delta g$ of $N$ cells is built with a regular step
of about 6.5 km. Each cell contains $n_i$ points with $i =$
{1, 2, ..., $N$}. These $n_i$ points are replaced by one point
which location is given by the geometric barycenter of the $n_i$
points, in the case that $n_i > 0$. If $n_i = 0$, then there is no
point in the cell $i$. Figure 7 show the new distributions of
gravimetric data for the Massif Central and the Alps regions;
they have, respectively, 4374 and 4959 location points. These
new spatial distributions reflect the initial BGI gravity data
distribution but are be more homogeneous. They will be used
in what follows.

Chronometric location points selection. We choose to put
clock measurements only where existing land gravity data
are located. Indeed, these data mainly follow the roads and
valleys which could be accessible for a clock comparison.
Then, we use a simple geometric approach in order to put
clock measurements in regions where the gravity data cover-
age is poor. Since the potential varies smoothly compared to
the gravity field, a clock measurement is affected by masses
at a larger distance than in the case of a gravimetric measure-
ment. For that reason, a clock point will be able to constrain
longer wavelengths of the geopotential than a gravimetric
point. This is particularly interesting in areas poorly sur-
veyed by gravity measurement networks. Finally, in order
to avoid having clocks too close to each other, we define a
minimal distance *d* between them. We chose *d* greater than
the correlation length of the gravity covariance function (in
this work *λ* ~ 20 km, see Table 1).

Here we give more details about our algorithm to select
the clock locations:

1. First, we initialize the clock locations on the nodes of a regular grid $\overline{T}$ with a fixed interval $d$. This grid is included in the target region at a setback distance of about 30 km from each edge (outside possible boundary effects).

2. Secondly, we change the positions of each clock point to the position of the nearest gravity point from the grid $\delta g$, located in cell $i$ (see the previous paragraph); in cell $i$ are located $n_i$ points of the initial BGI gravity data distribution.

3. Finally, we remove all the clock points located in cells where $n_i > n_{max}$. This is a simple way to keep only the clock points located in areas with few gravimetric measurements.

This method allows to simulate different realistic clock mea-
surement coverages by changing the values of *d* and *n*<sub>max</sub>.
The number of clock measurements increases when the dis-
tance *d* decreases or when the threshold *n*<sub>max</sub> increases and
vice versa. It is also possible to obtain different spatial dis-
tributions but the same number of clock measurements for
different sets of *d* and *n*<sub>max</sub>.

In Fig. 7, we propose an example of clock coverage
used hereafter for both target regions with 32 and 33 clock
locations, respectively, in the Massif Central and the Alps,
corresponding to ~0.7% of the gravity data coverage. For
the chosen distributions, the value of *d* is about 60 km
and *n*<sub>max</sub> = 15.

Synthetic measurements simulation For each data point,
the synthetic values of δg and T are computed by applying
the data reduction presented in Sect. 3.2. It is important to
note that the location points of the simulated data T are not
necessarily at the same place than the estimated data T̄.

A Gaussian white noise model is used to simulate the
instrumental noise of the measurements. We chose, for the
main tests in the next section, a standard deviation σ<sub>δg</sub> =
1 mGal for the gravity data and σ<sub>T</sub> = 0.1 m<sup>2</sup>/s<sup>2</sup> for the
potential data. In terms of geoid height, the latter noise level
is equivalent to 1 cm. Other tests with different noise levels
are discussed in Sect. 6.

**5 Numerical results**

In this section, we present our numerical results showing the
contribution of clock data in regional recovery of the geoapo-
tential from realistic data points distribution in the Massif
Central and the Alps. The reconstruction of the disturbing
potential is realized from the synthetic measurements δg
and T, and by applying the least squares collocation (LSC)
method.

Planar Least Squares Collocation. The LSC method,
described in Moritz (1972, 1980), is a suitable tool in geodesy
to combine heterogeneous data sets in gravity field modeling.
Assuming that the measured values are linear functionals of the disturbing potential $T$, this approach allows us to estimate any gravity field parameter based on $T$ from many types of observables.

Consider $\mathbf{I} = [\mathbf{I}_T, \mathbf{I}_{\delta g}] = l_k$ a data vector composed by $p$ data $T$ and $q$ data $\delta g$, affected by measurement errors $\varepsilon_k$, with $k = \{1, 2, \dots, p+q\}$. The estimation of the disturbing potential $\tilde{T}_P$ at point $P$ from the data $\mathbf{I}$ can be performed with the relation

$$ \tilde{T}_P = \mathbf{C}_{T_P, l} \cdot \mathbf{C}_{n, n}^{-1} \cdot \mathbf{I} \qquad (5.1) $$

$$ \mathbf{C}_{n, n} = \mathbf{C}_{l, l} + \omega \mathbf{C}_{\varepsilon, \varepsilon} \qquad (5.2) $$

with $\mathbf{C}_{l,l}$ the covariance matrix of the measurement vector $\mathbf{I}$, $\mathbf{C}_{\varepsilon,\varepsilon}$ the covariance matrix of the noise, $\mathbf{C}_{T_P,l}$ the cross-covariance matrix between the estimated signal $T_P$ and the data $\mathbf{I}$, and $\omega$ the Tikhonov regularization factor (Neyman 1979), also called weight factor.

In practice, the data $\mathbf{I}$ are synthesized as described in Sects. 3 and 4. Therefore, the measurement noise is known to be a Gaussian white noise. Noise and signal (errorless part of $l_k$) are assumed to be uncorrelated, and the covariance matrix of the noise can be written as

$$ \mathbf{C}_{\varepsilon, \varepsilon} = \begin{bmatrix} \mathbf{I}_p \cdot \sigma_T^2 & 0 \\ 0 & \mathbf{I}_q \cdot \sigma_{\delta g}^2 \end{bmatrix} \qquad (5.3) $$

with $\mathbf{I}_n$ the identity matrix of size $n$.

Because $\mathbf{C}_{l,l}$ can be very ill-conditioned, the matrix (5.3) plays an important role in its regularization before inversion, since positive constant values are added to the elements of its main diagonal. To avoid any iterative process to find an optimum value of $\omega$ in case where this matrix $\mathbf{C}_{l,l}$ is not definite positive, we chose to fix the weight factor $\omega = 1$ and to apply a singular value decomposition (SVD) to pseudo-inverse the matrix. As shown in (Rummel et al. 1979), these two approaches are similar.

Estimation of the covariance function. Implementation of the collocation method requires to compute the covariance matrices $\mathbf{C}_{T_P,l}$ and $\mathbf{C}_{l,l}$. This step has been carried out using a logarithmic spatial covariance function from (Forsberg 1987), see “A Covariance function.” This stationary and isotropic model is well adapted to our analysis. Indeed, it provides the auto-covariances (ACF) and cross-covariances (CCF) of the disturbing potential $T$ and its derivatives in 3 dimensions with simple closed-form expressions.

The spatial correlations of the gravity field are analyzed with the program GPFIT (Forsberg and Tscherning 2008). The variance $C_0$ is directly computed from the gravity data on the target area, and the parameters $\alpha$ and $\beta$ (see “A Covariance function”) are estimated by fitting the a priori covariance function to the empirical ACF of the gravity disturbances $\delta g$.

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 8** Empirical and best fitting covariance function of the ACF of $\delta g$. Values of the parameters are given in Table 1. **a** Massif Central, **b** Alps-Mediterranean</figcaption></figure>

Results of the optimal regression analysis for both regions are given in Fig. 8 and Table 1. The estimated covariance models reflect the different characteristics of the gravity signals in the two areas and the data sampling, which is less dense in high relief areas. Finally, the gravity anomaly covariances show similar correlation lengths, with a larger variance for the case of the Alps; their shapes, however, slightly differ, with a broader spectral coverage for the Alps.

Knowing the parameter values of the covariance model, we can now estimate the potential anywhere on the Earth’s surface.

Contribution of clocks. The contribution of clock data in the potential recovery is evaluated by comparing the residuals of two solutions to the reference potential on a regular grid interval of 10 km. The first solution corresponds to the errors between the estimated potential model computed solely from gravity data and the potential reference model, while the second solution uses combined gravimetric and clock data. To
<table><caption>**Table 1** Estimation of the auto-covariance function parameters on the gravity data δg using the logarithmic model from Forsberg (1987) with, μ the mean, C₀ the variance, α and β, respectively, a shallow and a compensating depth parameter</caption><thead><tr><th>Area</th><th>Nb data</th><th>μ (mGal)</th><th>C₀ (mGal²)</th><th>α (km)</th><th>β (km)</th><th>λ [km]</th></tr></thead><tbody><tr><td>Massif Central</td><td>4374</td><td>0.41</td><td>63.4</td><td>24</td><td>15</td><td>21</td></tr><tr><td>Alps-Med.</td><td>4959</td><td>1.15</td><td>352.5</td><td>6</td><td>47</td><td>18</td></tr></tbody></table>

Here, λ is the correlation length defined as the distance at which the covariance is half of the variance

<figure><img src="image_14.png" /><figcaption>**Fig. 9** Accuracy of the disturbing potential *T* reconstruction on a regular 10-km step grid in Massif Central, obtained by comparing the reference model and the reconstructed one. In **a**, the estimation is realized</figcaption></figure>

from the 4374 gravimetric data δg only and in **b** by adding 33 potential data *T* to the gravity data. **a** Without clock data, **b** With clock data

avoid boundary effects in the estimated potential recovery, a grid edge cutoff of 30 km has been removed in the solutions.

For the Massif Central region, the disturbing potential is estimated with a bias μT ≈ 0.041 m²s⁻² (4.1 mm) and a rms σT ≈ 0.25 m²s⁻² (2.5 cm) using only the 4374 gravimetric data, see Fig. 9a. When we now reconstruct *T* by adding the 33 potential measurements to the gravimetric measurements, the bias is improved by one order of magnitude (μT ≈ −0.002 m²s⁻² or −0.2 mm) and the standard deviation by a factor 3 (σT ≈ 0.07 m²s⁻² or 7 mm), see Fig. 9b.

For the Alps, Fig. 10, the potential is estimated with a bias μT ≈ 0.23 m²s⁻² (2.3 cm) and a standard deviation

σT ≈ 0.39 m²s⁻² (3.9 cm) using only the 4959 gravimetric data. When adding the 32 potential measurements, we note that the bias is improved by a factor 4 (μT ≈ −0.069 m²s⁻² or −6.9 mm) and the standard deviation by a factor 2 (σT ≈ 0.18 m²s⁻² or 1.8 cm).

It can be noticed that the residuals in both areas differ. This results from the covariance function that is less well modeled when the data survey has large spatial gaps. It should also be stressed that a trend appears in the reconstructed potential with respect to the original one when no clock data are added in both regions. This effect is discussed in Sect. 6.
<figure><img src="image_15.png" /><figcaption>**Fig. 10** Accuracy of the disturbing potential *T* reconstruction on a regular 10-km step grid in Massif Central, obtained by comparing the reference model and the reconstructed one. In **a**, the estimation is real-</figcaption></figure>

ized from the 4959 gravimetric data $\delta g$ only and in **b** by adding 32 potential data *T* to the gravity data. **a** Without clock data, **b** With clock data

# 6 Discussion

Effect of the number of clock measurements.

Figure 11 shows the influence of the number of clock data in the potential recovery, and therefore, of their spatial distribution density. We vary the number and distribution of clock data by changing the mesh grid size *d*, which represents the minimum distance between clock data points (see Sect. 4). The particular cases shown in detail in Sect. 5 are included. We characterize the performance of the potential reconstruction by the standard deviation and mean of the differences between the original potential on the regular grid and the reconstructed one. When increasing the density of the clock network, the standard deviation of the differences tends toward the centimeter-level, for the Massif Central case, and the bias can be reduced by up to 2 orders of magnitude. Note that we have not optimized the clock locations such as to maximize the improvement in potential recovery. The chosen locations are simply based on a minimum distance

and a maximum coverage of gravity data (c.f. Sect. 4). An optimization of clock locations would likely lead to further improvement, but is beyond the scope of this work and will be the subject of future studies.

Moreover, the results indicate that it is not necessary to have a large number of clock data to improve the reconstruction of the potential. We can see that only a few tens of clock data, i.e., less than 1% of the gravity data coverage, are sufficient to obtain centimeter-level standard deviations and large improvements in the bias. When continuing to increase the number of clock data, the standard deviation curve seems to flatten at the cm-level.

Effect of the number of gravity measurements.

We have performed numerical tests in order to study the influence of the density of gravity measurements on the reconstructed disturbing potential, with or without clocks. We take the case of the Massif Central region and set up simulations where the clock coverage is fixed (either no clocks, or 38 clocks at fixed locations where we also have gravity data).
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 11** Performance of the potential reconstruction (expressed by the standard deviations and mean of differences between the original potential on the regular grid and the reconstructed one) wrt the number of clocks. In green, number of clock data in terms of percentage of $\delta g$ data. **a** Massif Central area, **b** Alps area</figcaption></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 12** Effect of the number of gravity data combined with 38 clock data on the disturbing potential recovery in the Massif Central region. Panel **a**: absolute value of the mean of the residuals of $T$; panel **b**: the rms. The noise of the measurements is 1 mGal for $\delta g$ and $0.1 \text{ m}^2\text{s}^{-2}$ for $T$. Note that for each coverage of gravity data, a new covariance model is fitted on the empirical covariance model</figcaption></figure>

Then, we progressively increase the spatial resolution of the gravity data, from 91 to 6889 points, and evaluate as before the quality of the potential reconstruction with or without clocks. Here, in contrast with the tests presented in the previous section, the gravity points are randomly generated from a complete 5-km step grid. Figure 12 shows the results of these tests. If we compare the rms values between configurations where we add clocks or not, we observe that the behavior of the results is globally similar and improved with clocks. The interpolation error due to a too low resolution of the gravity data with respect to scales of the field variations predominates when we have less than ~1500 gravity measurements, leading to large rms values even with clocks. Above this number, the large-scale reconstruction errors significantly contribute to the rms of residuals, explaining that the rms

further decreases only when clocks are added. Looking at the bias between the reconstructed and original potential, we can see that it is poorly dependent on the number of gravity data in the tests without clocks. It probably reflects the fact that these data are more sensitive to the smaller scale components of the gravity potential. When we add clocks, the improvement on the bias is always important, which is consistent with the fact that the higher sensitivity of clocks to the longer wavelengths of the field reduces significantly the trend from the modeling error.

Covariance function consistency.

In Figs. 10a and 11a, a trend appears in the residuals, but disappears when gravimetric and clock data are combined. This is due to the fact that the covariance function
<table><caption>**Table 2** Noise level effect on the disturbing potential recovery in the Massif Central region</caption><thead><tr><th rowspan="3">σ<sub>T</sub></th><th colspan="8">σ<sub>δg</sub></th></tr><tr><th colspan="2">5 mGal</th><th colspan="2">1 mGal</th><th colspan="2">0.1 mGal</th><th colspan="2">0.01 mGal</th></tr><tr><th>μ</th><th>σ</th><th>μ</th><th>σ</th><th>μ</th><th>σ</th><th>μ</th><th>σ</th></tr></thead><tbody><tr><td>No clock</td><td>2.2 × 10<sup>-1</sup></td><td>3.7 × 10<sup>-1</sup></td><td><b>4.1 × 10<sup>-2</sup></b></td><td><b>2.5 × 10<sup>-1</sup></b></td><td>1.5 × 10<sup>-1</sup></td><td>1.7 × 10<sup>-1</sup></td><td>2.6 × 10<sup>-1</sup></td><td>1.8 × 10<sup>-1</sup></td></tr><tr><td>1 m<sup>2</sup>s<sup>-2</sup></td><td>-4.4 × 10<sup>-3</sup></td><td>2.8 × 10<sup>-1</sup></td><td>-1.8 × 10<sup>-4</sup></td><td>1.7 × 10<sup>-1</sup></td><td>-1.1 × 10<sup>-2</sup></td><td>1.6 × 10<sup>-1</sup></td><td>-2.0 × 10<sup>-2</sup></td><td>1.7 × 10<sup>-1</sup></td></tr><tr><td>0.1 m<sup>2</sup>s<sup>-2</sup></td><td>-1.4 × 10<sup>-2</sup></td><td>2.0 × 10<sup>-1</sup></td><td><b>-2.4 × 10<sup>-3</sup></b></td><td><b>7.3 × 10<sup>-2</sup></b></td><td>-6.7 × 10<sup>-3</sup></td><td>5.2 × 10<sup>-2</sup></td><td>-1.1 × 10<sup>-3</sup></td><td>4.8 × 10<sup>-2</sup></td></tr></tbody></table>

In bold: results presented in Sect. 5. Values are given in m²s⁻²

<table><caption>**Table 3** Noise level effect on the disturbing potential recovery in the Alps region</caption><thead><tr><th rowspan="3">σ<sub>T</sub></th><th colspan="8">σ<sub>δg</sub></th></tr><tr><th colspan="2">10 mGal</th><th colspan="2">1 mGal</th><th colspan="2">0.1 mGal</th><th colspan="2">0.01 mGal</th></tr><tr><th>μ</th><th>σ</th><th>μ</th><th>σ</th><th>μ</th><th>σ</th><th>μ</th><th>σ</th></tr></thead><tbody><tr><td>No clock</td><td>5.8 × 10<sup>-1</sup></td><td>6.6 × 10<sup>-1</sup></td><td><b>2.2 × 10<sup>-1</sup></b></td><td><b>3.9 × 10<sup>-1</sup></b></td><td>2.1 × 10<sup>-1</sup></td><td>4.2 × 10<sup>-1</sup></td><td>2.1 × 10<sup>-1</sup></td><td>4.2 × 10<sup>-1</sup></td></tr><tr><td>1 m<sup>2</sup>s<sup>-2</sup></td><td>1.8 × 10<sup>-1</sup></td><td>6.2 × 10<sup>-1</sup></td><td>1.4 × 10<sup>-1</sup></td><td>3.4 × 10<sup>-1</sup></td><td>1.2 × 10<sup>-1</sup></td><td>3.3 × 10<sup>-1</sup></td><td>1.2 × 10<sup>-1</sup></td><td>3.3 × 10<sup>-1</sup></td></tr><tr><td>0.1 m<sup>2</sup>s<sup>-2</sup></td><td>2.0 × 10<sup>-1</sup></td><td>5.6 × 10<sup>-1</sup></td><td><b>6.8 × 10<sup>-2</sup></b></td><td><b>1.7 × 10<sup>-1</sup></b></td><td>4.7 × 10<sup>-2</sup></td><td>1.5 × 10<sup>-1</sup></td><td>1.7 × 10<sup>-2</sup></td><td>1.6 × 10<sup>-1</sup></td></tr></tbody></table>

In bold: results presented in Sect. 5. Values are given in m²s⁻²

the data. As previously, numerical simulations are performed with and without adding clocks, and the results are presented in the first column of Tables 2 and 3. We can see that decreasing the accuracy of the gravimetric measurements increases the residuals as compared to the previous solutions. This is due to the fact that the signal-to-noise ratio decreases, degrading the covariance function modeling. However, our previous conclusions on the benefit of clocks remain the same, even in the presence of significant signals at the shortest spatial scales.

## 7 Conclusions

Optical clocks provide a tool to measure directly the potential differences and determine the geopotential at high spatial resolution. We have shown that the recovery of the potential from gravity and clock data with the LSC method can improve the determination of geopotential at high spatial resolution, beyond what is available from satellites. Compared to a solution that does not use the clock data, the standard deviation of the disturbing potential reconstruction can be improved by a factor 3, and the bias can be reduced by up to 2 orders of magnitude with only a few tens of clock data. This demonstrates the benefit of this new potential geodetic observable, which could be put in practice in the medium term when the first transportable optical clocks and appropriate time transfer methods will be developed (see Bongs 2015; Lisdat et al. 2016; Deschênes et al. 2016; Vogt et al. 2016). Since clocks are sensitive to low frequencies of the gravity field, this method is particularly well adapted in

hilly and mountainous regions for which the gravity coverage is more sparsely distributed, allowing to fill areas not covered by the classical geodetic observables (gravimetric measurements). Additionally, adding new observables helps to reduce the modeling errors, e.g., coming from a mismatch between the covariance function used and the real gravity field.

In the same way, GPS and leveling data have been used, in combination with gravity data, to derive high-resolution gravimetric geoids (Kotsakis and Sideris 1999; Duquenne 1999; Denker et al. 2000; Duquenne et al. 2005; Nahavandi and Soltanpour 2006). Using clocks is, however, different from performing GPS and leveling measurements. They provide an information of similar nature as the gravity data, in contrast with these geometric observations. The latter are affected by different sources of errors (e.g., Duquenne 1998; Marti et al. 2001) and quite expensive in the case of leveling campaigns. We can expect that clocks could help identify and reduce errors in the gravity and GPS/leveling through their joint analysis for geopotential determination. Beyond the application considered in this work, the clocks can also contribute to the unification of height systems realizations (Shen et al. 2011, 2016; Denker 2013; Kopeikin et al. 2016; Takano et al. 2016), connecting distant points to a high-resolution reference potential network.

To our knowledge, this is the first detailed quantitative study of the improvement in field determination that can be expected from chronometric geodesy observables. It provides first estimates and paves the way for future more detailed and in depth works in this promising new field.
To overcome some limitations in the a priori model, as discussed in the previous section, we intend in a forthcoming work to investigate in more details the imperfections of the covariance function model. Moreover, as the gravity field is in reality non-stationary in mountainous areas or near the coast, some numerical tests with non-stationary covariance functions will be conducted. Another promising source of improvement could be the optimization of the positioning of the clock data. For example, the correlation lengths and the variations of the gravity field could be used as constraints. A genetic algorithm could also be considered to solve this location problem. Finally, it will be interesting to focus on the improvement of the potential recovery quality by combining other types of observables such as leveling data and gradiometric measurements. As knowledge of the geopotential provides access to height differences, this could be a way to estimate errors of the GNSS technique for the vertical positioning or contribute to regional height systems unification.

**Acknowledgements** We thank René Forsberg for providing us the FORTRAN code of the logarithmic covariance function. We gratefully acknowledge financial support from Labex FIRST-TF, ERC AdOC (Grant No. 617553 and EMRP ITOC (EMRP is jointly funded by the EMRP participating countries within EURAMET and the European Union). We thank Olivier Jamet and Matthias Holschneider for discussions about the collocation method. We thank Gwendoline Pajot-Métivier for discussions on high-resolution gravity signals. We thank anonymous reviewers and the associate editor for their useful comments on this manuscript.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

## Appendix: Covariance function

Let us consider two points P and Q with the Cartesian coordinates $(x_P, y_P, z_P)$ and $(x_Q, y_Q, z_Q)$, respectively. To compute the ACF and CCF of the disturbing potential T and its derivatives, Forsberg (1987) proposed a planar attenuated logarithm covariance model with upward continuation that can be expressed in the generic form

$$C(x, y, z_1 + z_2) = S \sum_{i=0}^{3} \lambda_i K(x, y, z_i) \quad (7.1)$$

with

$$x = x_Q - x_P, \quad y = y_Q - y_P \quad (7.2a)$$

$$z_i = z_P + z_Q + \alpha_i \quad (7.2b)$$

$$\alpha_i = \alpha + i\beta \quad (7.2c)$$

$$\lambda_i = \{1, -3, 3, -1\} \quad (7.2d)$$

$$S = C_0 \log^{-1} \left( \frac{\alpha_1^3 \alpha_3}{\alpha_0 \alpha_2^3} \right) \quad (7.2e)$$

This model is characterized by three parameters: $C_0$ the variance of the gravity disturbance $\delta g$ and two scale factors acting as high and low-frequency attenuators: $\alpha$ the shallow depth parameter and $\beta$ the compensating depth, respectively. The function $K = K(x, y, z_i)$ is logarithmic function modeling the covariances between the gravity field quantities. For example, by putting $r_i = \sqrt{d^2 + \alpha_i^2}$ and $d = \sqrt{x^2 + y^2}$, the ACF of $\delta g$ and T can be evaluated, respectively, with

$$K = -\log(\alpha_i + r_i) \quad (7.3)$$

$$K = \frac{3}{4} z_i r_i + \left( \frac{r_i^2}{4} - \frac{3}{4} z_i^2 \right) \log(z_i + r_i) \quad (7.4)$$

and the CCF between T and $\delta g$ with

$$K = r_i - z_i \log(z_i + r_i) \quad (7.5)$$

## References

* Bjerhammar A (1985) On a relativistic geodesy. Bull Déd 59(3):207–220. doi:10.1007/BF02520327

* Bondarescu R, Bondarescu M, Hetényi G, Boschi L, Jetzer P, Balakrishna J (2012) Geophysical applicability of atomic clocks: direct continental geoid mapping. Geophys J Int 191(1):78–82. doi:10.1111/j.1365-246X.2012.05636.x

* Bondarescu M, Bondarescu R, Jetzer P, Lundgren A (2015a) The potential of continuous, local atomic clock measurements for earthquake prediction and volcanology. In: European Physical Journal Web of Conferences, European Physical Journal Web of Conferences, vol 95, p 4009, doi:10.1051/epjconf/20159504009, arXiv:1506.02853

* Bondarescu R, Schärer A, Jetzer P, Angélil R, Saha P, Lundgren A (2015b) Testing general relativity and alternative theories of gravity with space-based atomic clocks and atom interferometers. In: European Physical Journal Web of Conferences, European Physical Journal Web of Conferences, vol 95, p 2002, doi:10.1051/epjconf/20159502002, arXiv:1412.2045

* Bondarescu R, Schärer A, Lundgren A, Hetényi G, Houlié N, Jetzer P, Bondarescu M (2015c) Ground-based optical atomic clocks as a tool to monitor vertical surface motion. Geophys J Int 202:1770–1774. doi:10.1093/gji/ggv246. arXiv:1506.02457

* Bongs K, Singh Y, Smith L, He W, Kock O, Świerad D, Hughes J, Schiller S, Alighanbari S, Origlia S, Vogt S, Sterr U, Lisdat C, Targat RL, Lodewyck J, Holleville D, Venon B, Bize S, Barwood GP, Gill P, Hill IR, Ovchinnikov YB, Poli N, Tino GM, Stuhler J, Kaenders W (2015) Development of a strontium optical lattice clock for the SOC mission on the ISS. C R Phys 16(5):553–564. doi:10.1016/j.crhy.2015.03.009

* Bruinsma SL, Förste C, Abrikosov O, Lemoine JM, Marty JC, Mulet S, Rio MH, Bonvalot S (2014) Esa’s satellite-only gravity field model via the direct approach based on all geoid data. Geophys Res Lett 41(21):7508–7514. doi:10.1002/2014GL062045L062045

* Castaldo R, Fedi M, Florio G (2014) Multiscale estimation of excess mass from gravity data. Geophys J Int p ggu082
Chou CW, Hume DB, Rosenband T, Wineland DJ (2010) Optical clocks and relativity. *Science* 329(5999):1630–1633. doi:10.1126/ science.1192720

Claessens SJ, Hirt C (2013) Ellipsoidal topographic potential: new solutions for spectral forward gravity modeling of topography with respect to a reference ellipsoid. *J Geophys Res Solid Earth* 118(11):5991–6002. doi:10.1002/2013JB010457B010457

Delva P, Lodewyck J (2013) Atomic clocks: new prospects in metrology and geodesy. *Acta Futura*, Issue 7, p 67–78 7:67–78, arXiv:1308.6766

Denker H (2013) Regional gravity field modeling: theory and practical results. Springer, Berlin. doi:10.1007/978-3-642-28000-9_5

Denker H, Torge W, Wenzel G, Ihde J, Schirmer U (2000) Investi- gation of different methods for the combination of gravity and gps/levelling data. In: Geodesy Beyond 2000, Springer, Berlin. pp 137–142

Deschênes JD, Sinclair LC, Giorgetta FR, Swann WC, Baumann E, Bergeron H, Cermak M, Coddington I, Newbury NR (2016) Syn- chronization of distant optical clocks at the femtosecond level. *Phys Rev X* 6(021):016. doi:10.1103/PhysRevX.6.021016

Douch K, Panet I, Pajot-Métivier G, Christophe B, Foulon B, Lequentrec-Lalancette MF, Diament M (2015) Error analysis of a new planar electrostatic gravity gradiometer for airborne surveys. *J Geod* 89:1217–1231. doi:10.1007/s00190-015-0847-8

Duquenne H (1998) Qgf98, a new solution for the quasigeoid in France. In: *Proceeding of the Second Continental Workshop on the Geoid in Europe. Reports of the Finnish Geodetic Institute*, vol 98, pp 251–255

Duquenne H (1999) Comparison and combination of a gravimetric quasigeoid with a levelled gps data set by statistical analysis. *Phys Chem Earth Part A Solid Earth Geod* 24(1):79–83. doi:10.1016/ S1464-1895(98)00014-3

Duquenne H, Everaerts M, Lambot P (2005) Merging a gravimetric model of the geoid with GPS/levelling data: an example in Bel- gium. Springer, Berlin. doi:10.1007/3-540-26932-0_23

Forsberg R (1987) A new covariance model for inertial gravimetry and gradiometry. *J Geophys Res* 92:1305–1310. doi:10.1029/ JB092iB02p01305

Forsberg R, Tscherning CC (2008) An overview manual for the GRAV- SOFT. University of Copenhagen, Denmark

Förste C, Bruinsma S, Abrikosov O, Flechtner F, Marty JC, Lemoine JM, Dahle C, Neumayer H, Barthelmes F, König R, Biancale R (2014) EIGEN-6C4 - The latest combined global gravity field model including GOCE data up to degree and order 1949 of GFZ Potsdam and GRGS Toulouse. In: EGU General Assembly Con- ference Abstracts, EGU General Assembly Conference Abstracts, vol 16, p 3707

Gruber T (2009) Evaluation of the egm2008 gravity field by means of gps-levelling and sea surface topography solutions. External quality evaluation reports of EGM08, Newton’s Bulletin 4, Bureau Gravimétrique International (BGI) / International Geoid Service (IGeS)

Hinkley N, Sherman JA, Phillips NB, Schioppo M, Lemke ND, Beloy K, Pizzocaro M, Oates CW, Ludlow AD (2013) An atomic clock with 10–18 instability. *Science* 341(6151):1215–1218. doi:10.1126/ science.1240420

Holschneider M, Chambodut A, Mandea M (2003) From global to regional analysis of the magnetic field on the sphere using wavelet frames. *Phys Earth Planet Inter* 135(2–3):107–124. doi:10.1016/ S0031-9201(02)00210-8

Johnson B (2009) Noaa project to measure gravity aims to improve coastal monitoring. *Science* 325(5939):378–378. doi:10.1126/ science.325_378

Knudsen P (1988) Determination of local empirical covariance func- tions from residual terrain reduced altimeter data. Tech. rep, DTIC Document

Knudsen P (1987) Estimation and modelling of the local empirical covariance function using gravity and satellite altimeter data. Bull *Géod* 61(2):145–160. doi:10.1007/BF02521264

Kopeikin SM, Kanushin VF, Karpik AP, Tolstikov AS, Gienko EG, Goldobin DN, Kosarev NS, Ganagina IG, Mazurova EM, Karaush AA, Hanikova EA (2016) Chronometric measurement of orthome- tric height differences by means of atomic clocks. *Gravit Cosmol* 22(3):234–244. doi:10.1134/S0202289316030099

Kotsakis C, Sideris MG (1999) On the adjustment of combined gps/levelling/geoid networks. *J Geod* 73(8):412–421

Landau L, Lifshitz EM (1975) *The Classical Theory of Fields. No. vol. 2* in Course of theoretical physics, Butterworth-Heinemann

Lisdat C, Grosche G, Quintin N, Shi C, Raupach SMF, Grebing C, Nicolodi D, Stefani F, Al-Masoudi A, Dörscher S, Häfner S, Robyr JL, Chiodo N, Bilicki S, Bookjans E, Koczwara A, Koke S, Kuhl A, Wiotte F, Meynadier F, Camisard E, Abgrall M, Lours M, Leg- ero T, Schnatz H, Sterr U, Denker H, Chardonnet C, Le Coq Y, Santarelli G, Amy-Klein A, Le Targat R, Lodewyck J, Lopez O, Pottie PE (2016) A clock network for geodesy and fundamen- tal science. *Nature Communications* 7:12,443 EP -, doi:10.1038/ ncomms12443

Mai E (2013) Time, atomic clocks, and relativistic geodesy. Deutsche Geodätische Kommission bei der Bayerischen Akademie der Wis- senschaften, Reihe A, *Theoretische Geodäsie*, Beck

Marti U (2015) Gravity, Geoid and Height Systems: Proceedings of the IAG Symposium GGHS2012, October 9-12, 2012, Venice, Italy. International Association of Geodesy Symposia, Springer, Berlin https://books.google.fr/books?id=2f8qBgAAQBAJ

Marti U, Schlatter A, Brockmann E (2001) Combining levelling with gps measurements and geoid information

Moritz H (1972) Advanced Least-squares Methods. Ohio State Univer- sity, Department of Geodetic Science, Ohio State University

Moritz H (1980) Advanced physical geodesy

Moritz H (1984) Geodetic reference system 1980. Bulletin géodésique 58(3):388–398. doi:10.1007/BF02519014

Nahavandchi H, Soltanpour A (2006) Improved determination of heights using a conversion surface by combining gravimet- ric quasi-geoid/geoid and gps-levelling height differences. *Studia Geophysica et Geodaetica* 50(2):165–180. doi:10.1007/ s11200-006-0010-3

Neyman YM (1979) The variational method of physical geodesy. Bul- letin géodésique. Nedra Publishers, Moscow

Nicholson TL, Campbell SL, Hutson RB, Marti GE, Bloom BJ, McNally RL, Zhang W, Barrett MD, Safronova MS, Strouse GF, Tew WL, Ye J (2015) Systematic evaluation of an atomic clock at $2 \times 10^{-18}$ total uncertainty. *Nat Commun* 6:6896. doi:10.1038/ ncomms7896. arXiv:1412.8261

Origlia S, Schiller S, Pramod MS, Smith L, Singh Y, He W, Viswam S, Świerad D, Hughes J, Bongs K, Sterr U, Lisdat C, Vogt S, Bize S, Lodewyck J, Le Targat R, Holleville D, Venon B, Gill P, Barwood G, Hill IR, Ovchinnikov Y, Kulosa A, Ertmer W, Rasel EM, Stuhler J, Kaenders W, SOC2 consortium contributors t (2016) Develop- ment of a strontium optical lattice clock for the SOC mission on the ISS. arXiv:1603.06062

Pail R, Bruinsma S, Migliaccio F, Förste C, Goiginger H, Schuh WD, Höck E, Reguzzoni M, Brockmann JM, Abrikosov O, Veichert s M, Fecher T, Mayrhofer R, Krasbutter I, Sansò F, Tsch- erning CC (2011) First goce gravity field models derived by three different approaches. *J Geod* 85(11):819–843. doi:10.1007/ s00190-011-0467-x

Pavlis NK, Weiss MA (2003) The relativistic redshift with $3 \times 10^{-17}$ uncertainty at nist, boulder, colorado, usa. *Metrologia* 40(2):66

Petit G, Wolf P, Delva P (2014) Atomic time, clocks, and clock com- parisons in relativistic spacetime: a review. In: Kopeikin SM (ed) *Frontiers in Relativistic Celestial Mechanics -, vol 2.*, Applications
and Experiments De Gruyter Studies in Mathematical Physics, De Gruyter, pp 249–279

Rummel R (2002) Global Unification of Height Systems and GOCE. Springer, Berlin. pp 13–20 doi:10.1007/978-3-662-04827-6_3

Rummel R, Schwarz KP, Gerstl M (1979) Least squares collocation and regularization. Bull Geod 53:343–361. doi:10.1007/BF02522276

Rummel R (2012) Height unification using GOCE. J Geod Sci 2:355–362. doi:10.2478/v10156-011-0047-2

Rummel R, Teunissen P (1988) Height datum definition, height datum connection and the role of the geodetic boundary value problem. Bull Géod 62(4):477–498. doi:10.1007/BF02520239

Sansò F, Sideris M (2013) Geoid Determination: Theory and Methods. Lecture Notes in Earth System Sciences, Springer, Berlin

Sansò F, Venuti G (2002) The height datum/geodetic datum problem. Geophys J Int 149(3):768–775. doi:10.1093/gji/149.3.768

Shen WB (2013) Orthometric height determination based upon optical clocks and fiber frequency transfer technique. In: 2013 Saudi International Electronics, Communications and Photonics Conference, pp 1–4. doi:10.1109/SIECPC.2013.6550987

Shen W, Ning J, Liu J, Li J, Chao D et al (2011) Determination of the geopotential and orthometric height based on frequency shift equation. Nat Sci 3(05):388

Shen Z, Shen WB, Zhang S (2016) Formulation of geopotential difference determination using optical-atomic clocks onboard satellites and on ground based on Doppler cancellation system. Geophys J Int. doi:10.1093/gji/ggw198

Shen Z, Shen W (2015) Geopotential difference determination using optic-atomic clocks via coaxial cable time transfer technique and a synthetic test. Geodesy and Geodynamics 6(5):344–350. doi:10.1016/j.geog.2015.05.012, http://www.sciencedirect.com/science/article/pii/S1674984715000816

Smith DA (1998) There is no such thing as “the” egm96 geoid: subtle points on the use of a global geopotential model. IGeS Bull 8:17–28

Smith WHF, Sandwell DT (1997) Global sea floor topography from satellite altimetry and ship depth soundings. Science 277(5334):1956–1962. doi:10.1126/science.277.5334.1956

Takano T, Takamoto M, Ushijima I, Ohmae N, Akatsuka T, Yamaguchi A, Kuroishi Y, Munekane H, Miyahara B, Katori H (2016) Geopotential measurements with synchronously linked optical lattice clocks. Nat Photon 10(10):662–666. doi:10.1038/nphoton.2016.159 letter

Tscherning CC (1976) Covariance expressions for second and lower order derivatives of the anomalous potential. Tech. Rep. 225, Ohio State University. Department of Geodetic Science

Tscherning CC, Rapp RH (1974) Closed covariance expressions for gravity anomalies, geoid undulations, and deflections of the vertical implied by anomaly degree variance models. Tech. Rep. 208, Ohio State University. Department of Geodetic Science

Tziavos IN, Vergos GS, Grigoriadis VN (2009) Investigation of topographic reductions and aliasing effects on gravity and the geoid over Greece based on various digital terrain models. Surv Geophys 31(1):23. doi:10.1007/s10712-009-9085-z

Vermeer M (1983) Chronometric Levelling. Reports of the Finnish Geodetic Institute, Geodeettinen Laitos, Geodetiska Institutet

Vogt S, Häfner S, Grotti J, Koller S, Al-Masoudi A, Sterr U, Christian L (2016) A transportable optical lattice clock. Journal of Physics: Conference Series 723(1):012,020, http://stacks.iop.org/1742-6596/723/i=1/a=012020

Yale MM, Sandwell DT, Herring AT (1998) What are the limitations of satellite altimetry? The Lead Edge 17(1):73–76

Zhang L, Li F, Chen W, Zhang C (2008) Height datum unification between Shenzhen and Hong Kong using the solution of the linearized fixed-gravimetric boundary value problem. J Geod 83(5):411. doi:10.1007/s00190-008-0234-9