Exploring potential applications of optical lattice clocks in a plate subduction zone

Optical clocks improved the accuracy of state-of-the-art cesium clocks by more than two orders of magnitude and enabled frequency comparison with a fractional uncertainty of one part in 1018. Gravitational redshift of two such clocks allows determining their height difference with an uncertainty of 1 cm. In Europe, chronometric leveling has been extensively conducted for unifying the height reference systems. Temporal response of the leveling, which affords monitoring a cm height variation within hours of averaging time, may offer new opportunities to explore seismology and volcanology. Superb stability of optical lattice clocks will be best used for such applications. This article outlines the prospects of chronometric leveling in Japan. Combining optical lattice clocks with an existing observation network of GNSS, crustal deformations may be monitored with unprecedented accuracy in the future.


Introduction
Atomic clocks are used in various applications, which are indispensable to support modern society, such as Global Navigation Satellite Systems (GNSSs), high-speed communications, and basic science. In the International System of Units (SI), one second is defined by cesium atomic clocks with an uncertainty of approximately 5 × 10 -16 , which corresponds to an uncertainty of 1 s in 60 million years.
Recent progress of optical atomic clocks has improved this uncertainty by more than 100 times and achieved fractional uncertainties of low 10 -18 (Nicholson et al. 2015;Ushijima et al. 2015;Brewer et al. 2019;McGrew et al. 2018). Such unprecedented accuracy of optical clocks allows not only testing fundamental laws of physics (Safronova et al., 2018), but also finds useful applications, including chronometric leveling (Delva et al. 2019) and the redefinition of the SI second (Riehle et al. 2018).
According to the theory of general relativity, clocks tick slower for deeper gravitational potential, known as gravitational redshift. Clocks with an uncertainty of 1 × 10 −18 enable us to measure height difference of a centimeter on the ground; clocks serve as novel quantum sensors for investigating gravitational potential or height difference (Müller et al. 2018). Such clock-based geodesy is referred to as chronometric leveling and opens new research fields (Mehlstäubler et al. 2018). For example, continental-scale fiber links of optical clocks (Lisdat et al. 2016) are expected to unify height reference systems across countries, which improves accuracy and cost for geodetic leveling (Denker et al. 2017). A more precise determination of the geoid using fiber-linked optical clocks in combination with traditional geodetic methods has also been studied (Delva et al. 2019). Theoretical frameworks for accurately evaluating relativistic effects on geodetic applications have also been developed, which are commensurate with improved observational precisions. It has been revealed that the geodetic measurements of the vertical displacements of the order of a few millimeters would require post-Newtonian effects to be considered in the description of the geoid (see, Sect. 7.3 of Kopeikin (2019)).
In contrast to these approaches regarding height reference systems in Europe, we have a distinct geophysical feature; the Japanese Island Arc located in a plate boundary has suffered from earthquakes and volcanic activities. As a result, geodesy in Japan has focused more on investigating minute temporal as well as spatial variations in the crust rather than resolving international inconsistency, as exemplified by a dense observation network for GNSS (GEONET: Tsuji and Hatanaka 2018) operated by the Geospatial Information Authority of Japan (GSI). To discuss the applicability of chronometric leveling to seismological and volcanological geodesy, stability and accuracy of atomic clocks play a key role.
This paper aims at extending the geodetic use of optical clocks to applications in geologically unstable regions in Japan. We first give an overview on the use of atomic clocks for the height and geoid determinations (Sect. 2) with an emphasis on the stability aspect of optical lattice clocks. A progress of chronometric leveling in Japan is mentioned therein. Next, we explore applications of optical lattice clocks to earthquake and volcano studies (Sects. 3-6). We emphasize that gravitational potential measurements can be used to observe crustal vertical motions. Ellipsoidal height observed by GNSS and potential change obtained by chronometric leveling are fundamentally different physical quantities since the latter reflects not only geometric changes but also mass changes. However, when considering some tectonic applications, the potential changes are dominated by the geometric effects, by which both quantities become approximately equal to each other. The performance of crustal deformation monitoring will be improved by networked optical clocks that allow accessing height uncertainty of 1 cm with an averaging time shorter than that required for GNSS.

Determination of the geoid
By definition, the Earth's gravitational potential consists of the potential due to the total mass of the Earth and the centrifugal potential due to the Earth's rotation. The geoid usually refers to the static equipotential surface, which best fits the long-term mean sea level.
There are three methods to determine the geoid, and each one has its own spatial resolution. The most traditional method is to use gravity data obtained on the ground and the ocean. The geodetic boundary value problem (GBVP) is solved with these data (Hofmann-Wellenhof and Moritz 2006), which is suitable to determine the local geoid with spatial resolutions ranging from 1 to 2 km to several tens of km.
The GNSS-leveling method takes the difference between the ellipsoidal height obtained by GNSS and the orthometric height obtained by leveling. This method allows us to directly measure the geoid height at a site (Hofmann-Wellenhof and Moritz 2006) without degrading the spatial resolution, but the measurement accuracy for leveling deteriorates for longer distances. (See Sect. 4.) Satellite gravity measurements are able to determine the geoid with spatial variations typically greater than 100 km because the potential change generated by the fine distribution of the topographic and other masses is significantly reduced at the satellite's altitude. Recent satellite-based geoid models have uncertainties of approximately 1.5 cm with a spatial resolution of 100 km (e.g., Delva et al. 2019).
These data with different spatial resolutions are often integrated, and the least-squares method constructs a geoid model that best explains all the employed data (Delva et al. 2019).

Chronometric approach
In contrast to the above methods based on Newton's gravitational theory, chronometric leveling relies on the relativistic effect. The idea of using clocks to determine the gravitational potential or height dates back to the 1980s (Bjerhammar 1985;Delva and Lodewyck 2013). In practice, we assume two distant clocks linked by a phase-stabilized optical fiber (Ma et al. 1994). Figure 1 shows a schematic for chronometric leveling. The ratio of the time rates and the gravitational potential is related by where , t, W , and c denote the proper time at the measurement site, the geocentric coordinate time related to an inertial system, the gravitational potential, and the speed of light, respectively (Müller et al. 2018). Due to this nature, the beat note of the clocks Δ = 2 − 1 measures the gravitational potential difference ΔW = W 2 − W 1 between two sites as Δ ∕ 1 = ΔW∕c 2 . Assuming that the gravitational acceleration at the two sites g 1,2 is equal to the average surface gravity g = 9.8m∕s 2 , ΔW can be approximated as gΔh where Δh = h 2 − h 1 . C o n s e q u e n t l y , w e o b t a i n Δ ∕ 1 ≈ 1.1 × 10 −18 Δhcm −1 and Δh = c 2 g Δ 1 ≈ 9.1 × 10 17 Δ 1 (cm) (see, Hofmann-Wellenhof and Moritz (2006) for the rigorous definition of g 1,2 and h 1,2 which are used in the geodetic determination of the potential difference). The fractional uncertainty for the measurement of the frequency difference is given by the quadratic sum, ≈ √ 2 1 + 2 2 + 2 link , where 1 , 2 and link denote fractional uncertainties of clock 1, clock 2, and optical frequency transfer, respectively. The corresponding height uncertainty is then given by h ≈ 9.1 × 10 17 cm. The magnitude of typically improves as 1∕ √ for an averaging time up to hours or days, which is termed as stability, and finally reaches floors determined by the clocks' uncertainties. For example, an uncertainty of ≈ 5 × 10 −16 achieved by cesium clocks corresponds to a height uncertainty of approximately h ≈ 5 m , which is inferior to conventional non-relativistic leveling. However, optical clocks with an uncertainty ≈ 1 × 10 −18 , which corresponds to a height uncertainty of 1 cm, will bring new avenue to chronometric leveling. Figure 2 shows the uncertainty of frequency comparison as a function of time, where the right axis shows the corresponding height uncertainty. Typical optical lattice clocks improve their stabilities with averaging time in the range of OLC = 1 × 10 −15 ( ∕s) −1∕2 up to 5 × 10 −17 ( ∕s) −1∕2 (Oelker Fig. 1 A schematic for chronometric leveling. On site 1 and 2, the gravitational potential is given by W 1 and W 2 , and the clock frequency with gravitational redshift is given by 1 and 2 , respectively. The clock signal is transferred by a phase-stabilized optical fiber (red line). The beat note of the clocks Δ = 2 − 1 measures gravitational potential difference ΔW = W 2 − W 1 between two sites as Δ ∕ 1 = ΔW∕c 2 . Assuming the gravitational acceleration g = g 1 = g 2 = 9.8m∕s 2 , the height difference Δh = h 2 − h 1 is given by Δh = c 2 g Δ 1 ≈ 9.1 × 10 17 Δ 1 (cm)  Takano et al. (2016). (b). future optical clocks with no-dead-time operation will improve the stability as −1 (c). operational magic condition may improve the clock uncertainty to low 10 -19 , and (d). an example of the uncertainty of the GNSS (precise point positioning) (e.g., Choy et al. 2017). (A-C) represent events which are expected to be (more precisely) detected by using the fiberlinked optical lattice clocks. (A) seasonal crustal movements including slow slip events (SSEs) and the secular deformation, (B) SSEs and (C) solid Earth and ocean tides, SSEs which are not yet discovered, possible slow slip prior to great earthquakes and volcanic eruptions with shorter durations than 1 day 93 Page 4 of 14 et al. 2019), which is mainly limited by the Dick effect caused by a periodic interrogation of the clock transition (Santarelli et al. 1998). Optical frequency transfer negligibly contributes to the overall uncertainties for > 100 s as it improves with link ∼ 1 × 10 −16 ( ∕s) −3∕2 (Newbury et al. 2007) or faster for fiber length of 100 km. As an example, we illustrate a clock stability = 6 × 10 −16 ( ∕s) −1∕2 (reddashed line) measured for a 15-km-long clock link (Takano et al. 2016).
As the long-term stability of the frequency transfer can be as low as low as 10 −21 (Akatsuka et al. 2020), atomic clocks essentially determine the uncertainty floor. Elimination of the blackbody radiation shift (BBR) improved clocks' uncertainties down to low 10 −18 (Nicholson et al. 2015;Ushijima et al. 2015;McGrew et al. 2018). At this level of uncertainties, evaluation of higher order light shifts becomes crucial. An operational magic condition is proposed  to cancel out higher order light shift to low 10 −19 , which is under investigation experimentally (Brown et al. 2017;Ushijima et al. 2018). In the future, no-dead-time operation of clocks (Dick et al. 1990;Schioppo et al. 2016;Katori 2021), which improves clocks' stability as −1 , will be used to access low 10 −19 uncertainty in an averaging time of less than an hour (dashed-green line).
Unlike geodetic leveling, the measurement uncertainty does not pile up for longer distances in chronometric leveling, as the uncertainty is guaranteed by that of atomic clocks. Moreover, the measurement time to achieve cm uncertainties can be less than a day and is nearly independent of distance for distance up to 100-1000 km, where optical frequency transfer uncertainty link can be less than that of clocks. These features are in striking contrast with geodetic leveling, requiring measurement time of about 2 km/ day. Gravitational potential differences obtained by chronometric leveling dramatically improve height measurements needed to establish a height reference system as well as provide us with a new observable for the GBVP to determine the geoid (Delva et al. 2019).

Progress of chronometric leveling
The national height reference systems adopted in European countries joined by land have discrepancies exceeding 10 cm in some places (e.g., Gruber et al. 2012). Consequently, one of the main tasks of the International Association of Geodesy (IAG) is to establish a globally unified height system. The IAG has formed a study group including geodetic applications of atomic clocks and an international project to realize chronometric leveling has been carried out in collaboration between physicists and geodesists (Margolis et al. 2013; http:// proje cts. npl. co. uk/ itoc/; https:// quge. iag-aig. org/). The earlier results are summarized in Delva et al. (2019).
In Japan, a high-resolution geoid model with an uncertainty of approximately 2 cm has already been constructed, thanks to the dense GNSS observation and leveling networks for crustal deformation monitoring (Miyahara et al. 2014). Nevertheless, chronometric leveling by optical lattice clocks allows constraining uncertainties in the GNSSleveling method, which accumulate as the measurement distance increases, further reducing the uncertainty of the long-wavelength geoid (Kuroishi 2017).
The first technical verification of chronometric leveling at cm level was reported for optical lattice clocks connected between Riken (Wako city, Saitama) and the University of Tokyo (Bunkyo-ku, Tokyo) with an uncertainty of 5 × 10 −18 (Takano et al. 2016). The gravitational redshift of the clocks was compared with the gravitational potential difference determined by a geodetic survey. Optical lattice clocks between Paris and Braunschweig demonstrated a continental scale frequency link with uncertainty 5 × 10 −17 of over 1000 km (Lisdat et al. 2016). Transportable optical lattice clocks (Koller et al. 2017) with an uncertainty of 7 × 10 −17 offer a new tool for chronometric leveling (Grotti et al. 2018). A pair of transportable optical lattice clocks demonstrated chronometric leveling with an uncertainty of 5 cm for a height difference of approximately 450 m at Tokyo SKYTREE (Takamoto et al. 2020).

Tectonic temporal variations in the Earth's gravitational field
The geoid model which is used as a height reference represents only the static component. However, the geoid as the actual equipotential surface changes in time. The main cause to disturb the gravitational potential at the surface is tides (Voigt et al. 2016). If the distance between two measurement sites for chronometric leveling is short, tidal effects are generally negligible because tides have relatively large wavelengths. However, the uncertainty caused by neglecting temporally varying tidal effects can reach 2 cm in terms of equivalent height difference, even when chronometric leveling is performed over a Japanese Island Arc scale (= 1000 km) (Kurosihi, 2017), which can be detected by optical lattice clocks (Fig. 2C). Such effects due to the solid Earth and ocean tides need to be accurately estimated, since it is common to remove tidal effects in geodetic applications including the determination of the static geoid. Examples of geophysical phenomena that cause temporal changes in the gravity field other than tides include variations in atmospheric pressure, continental water, non-tidal ocean bottom pressure (OBP), present-day ice mass changes and long-term viscoelastic responses of the mantle to deglaciation (= postglacial rebound) and great earthquakes (e.g., Wouters et al. 2014;Rodell et al. 2018) (Note that some of these effects are represented in terms of the equivalent water height anomaly (Wahr et al. 1998) instead of the geoid height anomaly in these papers). Furthermore, irregularities in the Earth's rotation cause not only classical but also relativistic changes in the gravity field (Fateev et al. 2015;Müller et al. 2018). These phenomena induce variations in the geoid with different space-time scales, but generally, the amplitudes are smaller than that for tides by an order of magnitude.
Now consider the extent that earthquakes and volcanic activities change the geometry of the masses and the gravity field over time. Here, it should be noted that the observed potential variations at the surface reflect not only the mass redistributions but also the height (or more general geometric) changes at which the clock is placed. The latter can be approximated by Eq. (1), where Δh now represents the temporal change in the relative height difference between two measurement sites.
Here we use the 2011 earthquake off the Pacific coast of Tohoku as an example of a great earthquake (M9). Figure 3 shows the geoid height change due to the earthquakeinduced mass redistribution as estimated by the method of Tanaka et al. (2015a). A decrease and an increase of 2 cm occurred at the Oshika Peninsula and off the southeast, respectively. At the Oshika Peninsula, subsidence by approximately 1 m was also observed by GNSS (Fig. 4). Assuming that the potential variation at this site is measured relative to a sufficiently far place where the geoid height change due to the earthquake is negligible, 98% of the measured potential change is due to the contribution of height change.
Next, as an example of a large eruption, consider the 1914 event of the Sakurajima volcano. Applying the simplest point-mass model and assuming that an eruptive volume of 1.2 km 3 (Yokoyama 2013) was lost from the magma chamber located at a 10-km depth, the estimated maximum decrease in the geoid is approximately 2 mm. A subsidence of approximately 80 cm was observed by a leveling survey, which was conducted before and after the eruption along the route surrounding Sakurajima volcano (Yokoyama 2013). For a rough estimate, it is adequate to consider that the subsidence and the geoid change were measured at the same location. Subsequently, the geoid height change due to mass redistribution is only a few percent of the contribution from the height change.
The geoid height changes due to a mass redistribution caused by an M9 earthquake or a large eruption are, at Reference site is the Fukue station in Nagasaki Prefecture, which is outside the figure.
(Lower) Time-series plot for the relative height change near the Mizusawa VLBI Observatory most, a few cm, but the apparent potential changes originating from height changes at the measurement site are much larger. The contribution of the mass redistribution caused by earthquakes with M ≤ 8 is less than 1 mm. That caused by secular vertical crustal motion due to the steady-state plate subduction is much smaller than 1 mm.
It is known that groundwater disturbances can affect terrestrial gravity measurements at 1-10 μGal (10-100 nm/s 2 ) levels (Crossley et al. 2013). However, the corresponding effects on the gravitational potential are generally below 0.1 mm (see, Sect. 2.7 and Fig. 3 of Mehlstäubler et al. 2018).

Measuring vertical crustal motions
Chronometric leveling mainly measures the apparent effect of the height change rather than the potential differences associated with the mass redistribution; therefore fiberlinked clocks can serve as an altimeter (e.g., Bondarescu et al. 2015). In this section, chronometric leveling using optical lattice clocks, geodetic leveling, and GNSS are compared as means to observe height changes.
The maximum allowable uncertainty of the first-order geodetic leveling over the measurement distance S is controlled with 2.5 √ (S∕km) mm in Japan. It follows that chronometric leveling with an uncertainty of 1 cm can be superior to geodetic leveling when the measurement involves distances longer than approximately 20 km. In addition, chronometric leveling with fiber-linked optical lattice clocks can greatly shorten the measurement time, as already mentioned.
Can optical lattice clocks measure the height changes better than GNSS? Before answering this question, consider the typical rate for the vertical crustal motion on the geodetic time scale. The rates derived from geodetic leveling and GNSS in the Japanese Islands on the whole are distributed within ± 1 cm/year (Murakami and Ozawa 2004). These rates are an order of magnitude larger than the average rates in most areas in Europe, except for Northern Europe where postglacial rebounds are remarkable (Altamimi et al. 2016). In addition, the 2011 M9 Tohoku earthquake triggered an ongoing, slow crustal uplift with rates exceeding 1 cm/year over eastern Japan (Fig. 4 (Lower)). The mechanism of this long-lasting deformation is considered to be a viscoelastic relaxation in the mantle, which gradually releases the earthquake-induced stress (e.g., Wang et al. 2012).
To discuss GNSS positioning uncertainty, Fig. 5a shows an example of the daily averaged coordinates in southwest (SW) Japan obtained by the GEONET. The static, doubledifference method was used for data analysis. The coordinates deviate 2-3 mm in the EW (red) and NS (green) components and approximately 10 mm in the UD component (blue). These deviations are mainly due to atmospheric noise. The deviations in the two horizontal components are smaller and mm-level variations are traceable. In the NS component, a slow slip event with a duration of one month (Heki and Kataoka 2008;Obara and Kato 2016), which is discussed later, is visible.
Generally, uncertainties for the GNSS-derived coordinates become larger for shorter measurement time.  (Zumberge et al. 1997) was used in the static precise point positioning mode with the International GNSS Service (IGS) Ultra-Rapid product. The other processing parameters were set to default. The standard deviation and the systematic bias with respect to the case for the 24-h average at the station #0751, where the above slow slip was detected, are 21 mm and 1.3 mm, respectively. When taking the difference from the station #0500, common errors (due to satellite orbits and atmospheric noise, for instance) which cannot be completely modeled are reduced and the standard deviation becomes smaller (16 mm). However, a systematic bias of + 7.2 mm remains with respect to the case of the 24-h average. The sum of these two uncertainties is larger than the uncertainty of 1 cm which fiber-linked optical lattice clocks can achieve within an hour (Fig. 2 b).
To extract rapid crustal vertical motion with shorter durations than a day, tidal effects must be removed. The correction for tides can be done by determining amplitude and phase of each tidal constituent using observed time-series data (e.g., with BAYTAP (Tamura et al. 1991)) or applying tidal prediction models such as GOTIC2 (Matsumoto et al. 2001). In the former approach based on observation data, a precision of the correction would depend on S/N (tidal signal/ observation noise). Considering a 100-km baseline, for example, tidal effects are several mm on both the height difference (Hatanaka et al. 2001) and the potential difference (in terms of height) (Kuroishi 2017). By the higher stability of the chronometric leveling for a diurnal band, less noise is expected for the chronometric approach than for the GNSS height measurement (Fig. 2 b). In the latter approach based on tidal models, Hatanaka et al. (2001) reported that the ocean tide model of GOTIC2 reproduces the GNSSderived vertical displacement for the M2 constituent with a mean uncertainty of 1 mm over the Japanese Islands. To what extent the same ocean tide model could reproduce the corresponding tidal gravitational potential change is not yet known. Chronometric leveling allows the potential change at such a local scale to be observed, which can be used to evaluate an uncertainty of the ocean tide model.

Atmospheric noise in GNSS data processing
It should be noted that the above uncertainty for GNSS positioning is evaluated after modeling and removing the ionospheric and tropospheric delays for the microwaves from satellites (Fig. 6). In general, the delays, which can be converted into the corresponding distance by multiplying by the speed of light, can exceed 5 m for the ionosphere and 2 m for the troposphere (e.g., Seeber (2003)). The effects of these delays can be reduced to a few cm using multiple frequency observations and inferring the atmospheric model parameters. A least-squares method is applied to solve the delays and the 3D coordinates of an observation site. Satellite positions move with time and the GNSS receiver receives satellite signals from different elevations and azimuth angles. Consequently, for longer measurement times (e.g., 24 h), most of the atmospheric delays cancel for the azimuth due to the symmetry of the satellite positions, resulting in smaller deviations in the horizontal components. However, the atmospheric delays do not completely cancel for the elevation because the receiver cannot see the satellites on the other side of the Earth. The remaining atmospheric delays can affect the obtained vertical coordinate due to the correlation between the associated parameters, which results in a larger deviation in the vertical component. To overcome this difficulty, a number of methods to correct for atmospheric delays have been proposed, including those based on improved parameterizations of the atmospheric effects and numerical weather models (Boehm and Schuh 2013). Another issue is the difficulty to evaluate uncertainty of GNSS-derived shorter-term coordinate changes using independent methods. Collocations with other space geodetic techniques have been elaborated, but they have usually lower temporal resolutions. As a result, uncertainty of GNSS coordinate time series is often evaluated with respect to the coordinate determined from a longer time-span GNSS observation data. The statistical average of atmospheric noise for a given time span will not always become zero. Hence, the precision does not necessarily improve for a higher sampling frequency in the GNSS observation.
The uncertainty in the vertical coordinate and the atmospheric delays can affect the horizontal coordinates. Figure 7d shows the daily displacements in the east-west (EW) and the vertical (UD) directions for 3 years obtained at a GEONET station in SW Japan near the Nankai Trough. In summer, the larger scattering of data is observed not only in the vertical coordinate (UD) but also in the horizontal coordinate (EW). Moreover, seasonal variations with amplitude of ~ 5 mm and ~ 10 mm are visible in the EW and the UD components, respectively, superimposed on the westward and downward secular motions caused by the inter-plate coupling. These seasonal variations can represent real crustal deformations due to atmosphere, ocean, groundwater and snow and/ or apparent deformations due to insufficient modeling of atmospheric noise (Munekane et al. 2008), and separation of real and apparent variations is not an easy task. We will discuss this issue in Sect. 6.2 in more detail.
6 Monitoring of crustal deformation by combining optical lattice clocks with GNSS

Separation of atmospheric noise from the vertical coordinate
Chronometric leveling based on fiber-linked optical lattice clocks is independent of the atmospheric noise that appear in the GNSS data analysis. Moreover, chronometric leveling allows measuring height changes with significantly shorter averaging time than GNSS. Incorporating optical lattice The blue dots in a and b show the daily coordinates in the vertical and the toward-trench (azimuth 120 • ) directions, respectively, obtained at a GNSS station #1189. The GIPSY X software was used in the analysis (static PPP with the IGS Final product). The other processing parameters were set to default. The location of the station is represented with the black square in (d).The red curves denote the deformation due to surface loads, derived from the GRACE gravity field data (Munekane et al. 2008). The color in c and d represent a surface displacement in the vertical and the toward-trench (azimuth 120 • ) directions, respectively, due to a seasonal slow slip event that is expected to occur during autumn to winter. The small boxes denote the locations of the faults with a uniform slip of 4 cm on the plate boundary. The dotted line denotes the Nankai trough clocks into a GNSS network greatly improves the monitoring performance for vertical crustal deformation. Furthermore, side-by-side observations of GNSS and optical lattice clocks can be used to verify atmospheric correction models adopted in the GNSS data analysis. The height changes that occur on short (< 1 day) or seasonal time scales obtained by chronometric leveling can be used as the true (known) value in the analysis, and the atmospheric model parameters can be inferred independently from the vertical coordinate. It should be noted that the effects of mass changes need to be carefully estimated and removed from the potential change observed by chronometric leveling, to compare with the geometric GNSS ellipsoidal height. After such treatments, it will be worth confirming whether the use of the obtained atmospheric parameters can reduce deviations in the timeseries horizontal coordinate data.

Detection of slow slip events
Fiber-linked optical lattice clocks may be used to detect slow slip events (SSEs). An SSE is caused by a slow fault slip on the plate interface. SSEs have been found to occur in the proximity of the source region of great earthquakes, and it has been pointed out that SSEs could trigger a great earthquake by providing a last push to the fault rupture (Obara and Kato 2016).
Using GNSS data, Nishimura (2014) detected 223 shortterm slow slip events in Japan over 14 years. Of these, 93 events were reported to contain large uncertainties. In the geophysical inversion for detecting SSEs, the contribution of GNSS vertical displacement data is much smaller than that of the horizontal displacement data, because the former has larger uncertainty. However, the amount of the vertical signal due to an SSE is comparable with that of the horizontal one [an example can be seen in Fig. 7c, d]. Chronometric leveling using fiber-linked optical lattice clocks will improve identification of SSEs which have been searched by GNSS, by making the most of information on crustal vertical deformation (Fig. 2B).
Furthermore, we expect fiber-linked optical clocks allow detection of SSEs with even shorter durations than a few days, which have not been discovered so far (Fig. 2C). The presence of such SSEs is expected by a scaling law for slow earthquakes (Ide et al. 2007). Finding those events will validate the scaling law and lead us to understand how a smaller rupture grows up to a larger earthquake. Sato (1977) reported that a slow coastal uplift with a rate of approximately 3 mm/hour initiated 1.5 days before the M w 8.3 Nankaido earthquake in 1946 at a tide gauge station. Linde and Sacks (2002) interpreted that this uplift was caused by a slow slip on the deep extension of the megathrust fault. The potential change expected for this slow slip is ~ 3 mm/hour in terms of height. By deploying a reference clock at a sufficiently far location from the coastal area (e.g., Kyoto), chronometric leveling can find the slow slip earlier than GNSS (Fig. 2C).
The fiber-linked optical lattice clocks may be used to identify SSEs on a seasonal time scale (Fig. 2A). It has been known that large earthquakes in the Nankai trough tend to occur in autumn and winter (e.g., Ohtake and Nakahra 1999), which may suggest the presence of SSEs that accelerate plate subduction during autumn to winter. Figure 7a, b shows the daily displacements observed at a GEONET station #1189, which is located near the Nankai trough, in the vertical and the towardtrench (azimuth 120 • ) directions, respectively. Linear trends are subtracted from the time-series plots. Assuming that a deeper part of the source region of an anticipated megathrust earthquake slips, we estimated a surface displacement caused by an SSE, using the method of Tanaka et al. (2015a). A uniform slip of 4 cm is given at depths of 31-43 km on the plate boundary (small boxes in Fig. 7d), and the slip is assumed to occur during autumn to winter so that it agrees with the above seismicity of large earthquakes. The moment magnitude is 7.0. We see from Fig. 7b that a toward-trench (i.e., positive) motion occurs from September to January, which agrees with the positive change expected by the SSE (Fig. 7d). However, the observed vertical motion from September to January is negative (Fig. 7a), which does not agree with uplift expected by the SSE (Fig. 7c).
On seasonal bands, variations in atmospheric pressure, terrestrial water and OBP also cause crustal deformations (Munekane et al. 2008;Schröder et al. 2021). Such deformations due to surface loads are dominated by vertical motion. Munekane et al. (2008) estimated annual vertical deformation due to surface loads, based on the variation in the gravity field observed with GRACE satellites. The GRACE-derived annual vertical deformation has amplitude of 3-4 mm and shows the maximum subsidence in the middle of October in SW Japan (Fig. 2b of their paper). We assumed that amplitude of the corresponding horizontal displacement is 1/3 of the vertical displacement (e.g., Chanard et al. 2018). This annual deformation is superimposed on Fig. 7a, b. We see that surface loading can explain the vertical component, but not the horizontal component.
The observed horizontal and vertical displacements suggest the SSE and surface loading, respectively. This contradiction may be resolved if the observed displacement is caused by a superposition of the SSE and surface loading. In such a case, we are able to separate these two effects using a network of optical clocks and GNSS as follows.
Observed temporal changes of the potential and height difference can be represented as (2) Δ obs = Δ tect + Δ load where Δ tect and Δh tect denote the effects due to tectonic events, including SSE, and Δ load and Δh load the effects of surface loading. In order to discuss seasonal changes, tidal effects are already reduced in these variables by applying tidal model corrections. Based on the discussion in Sect. 3, the potential changes due to mass redistributions are small for the tectonic effects, so where g denotes the average surface gravity. Suppose that surface loading is caused by a change in surface density L( , ) , which is distributed on a sphere with radius R (= Earth's mean radius), where and denote colatitude and longitude, respectively. This mass redistribution changes the surface potential by V( , ) = ∑ n 4 GR 2n+1 L n ( , ) , where L n represents spherical harmonics of L with degree n and G the universal gravitational constant (chapter 1 of Hofmann-Wellenhof and Moritz (2006); Sect. 6.5 of Dehant and Mathews (2015)). Using V n ≡ 4 GR 2n+1 L n , the terms representing surface loading in Eqs. (2) and (3) can be written as where k n and h n denote the load Love numbers for the potential change and the vertical displacement (Farrell 1972). Multiplying Eq. (3) with g and subtracting the result from Eq. (2) and using Eqs. (4-6), we obtain a new observable which depends only the effects of surface loading: The change in surface density L( , ) or L n ( , ) due to surface loads such as atmospheric, oceanic and hydrological loads can be estimated using GRACE data on the assumption that mass redistributions occur within a thin surface layer (Eq. (14) of Wahr et al. 1998). The spatial resolution of GRACE data for the monthly gravity field variation is ∼ 200 km, corresponding to a spherical harmonic degree of ∼ 100. This means that, using GRACE data, the potential change V n in Eq. (7) can also be inferred up to such a spatial scale. Annual surface loading that occurs on the Japanese Islands was estimated with GRACE data by Munekane et al. (2008), for instance.
If the observation on the left-hand side of Eq. (7) agrees with the right-hand side of Eq. (7) which is inferred from GRACE data, it follows that the left-hand side does not contain significant shorter-wavelength (< 200 km) surface loading. Then, because Δ load can be quantified using GRACE data and Eq. (5), we obtain the tectonic height change from Eq.
(2), based on chronometric leveling. If the observation on the left-hand side of Eq. (7) does not agree with the right-hand side of Eq. (7) which is inferred from GRACE data, it follows that the left-hand side contains shorter-wavelength surface loading that cannot be observed with a spatial resolution of GRACE data. In this case, V n for higher degrees in Eq. (5) need to be evaluated without using GRACE data. One possibility is to use a denser optical clock network to improve the space-time resolution for the detection of surface loading (Schröder et al. 2021). Another possibility is to use fine resolution OBP data, calculated by a dynamic ocean model assimilating real data (e.g., Wu et al. 2003;Usui et al. 2006;Forget et al. 2015;Tanaka et al. 2015b), although precision to represent OBP is generally limited because of the sparse distributions of real observations of OBP. Finally, subtracting the evaluated Δ load from Eq. (2), we could estimate the chronometric height change due to tectonic effects.

Improvement of volcano gravimetry
Monitoring vertical crustal motion is important for gravity observations to detect magma movement in volcanic areas. When a gravity measurement site moves vertically, an apparent gravity change arises ( ∼ −2 Gal/cm). As an example, Okubo (2020) estimates a gravity signal caused by magma movement in Asama volcano in Japan. The typical signal amounts to the apparent change due to a vertical motion of a few cm (Figs. 11 and 12 of the paper). Therefore, the vertical motion must be observed and the apparent change should be removed as precisely as possible. On August 15, 2015, a rapid dilatation of the Sakurajima volcano occurred in approximately 12 h (e.g., Okubo et al. 2017). This is attributed to magma intrusion. It is difficult to measure such a fast vertical motion at the gravity measurement site by GNSS with an uncertainty of 1 cm. However, this limitation can be overcome by using fiber-linked optical clocks (Fig. 2c). A detailed discussion for the detection of magma transport using optical clocks can be found in Bondarescu et al. (2015).

Detection of secular deformation
The rates of the long-term vertical crustal motion as observed in geological methods are on the order of 0.1 mm/ year (Goldfinger et al. 2013). When using the continuous GNSS data, an approximately 20-year observation span is necessary to detect the vertical motion with an uncertainty of 0.1 mm/year. Fiber-linked optical clocks can be used for such a long-term trend over tens of years as the longterm stability is guaranteed by atomic clocks. To elucidate how the land was formed, an accuracy of 0.1 mm/year is sometimes necessary to rigorously determine whether a site is uplifting or subsiding (Fig. 2a).

Summary and future plans
International collaborative studies to unify height reference systems using optical atomic clocks have been actively conducted in Europe. Besides such purely geodetic applications, we have discussed their potential applications in seismological and volcanological studies in Japan. When considering such tectonic applications, the gravitational potential difference measured by clocks is mainly given by the vertical displacement of clocks. The effect of the mass redistribution due to earthquakes and volcanic activities is much smaller than that of the vertical motion (at most a few percent). Thus, fiber-linked optical clocks can be used as an extremely accurate altimeter.
Continuous GNSS monitoring provides crustal deformation data that are indispensable as basic information for seismic and volcanic risk evaluations. However, it suffers from the atmospheric delay, which is inherent to space geodetic techniques. It is not easy to reduce the height uncertainty to a few mm in a 24 h average even after applying a correction model for the atmospheric delay. Uncertainties arising from incomplete corrections may prevail in the horizontal coordinates. They preclude us from extracting seasonal-scale crustal deformations.
To overcome this limitation, a terrestrial observation of height changes will be desirable, where signals do not propagate in the atmosphere. Fiber-linked optical clocks are best used for this purpose and allow determining height changes faster than GNSS. By combining such clock network with GNSS network, precise crustal deformation monitoring will be possible. Such concomitant use enables us to conduct studies to confirm atmospheric delay correction models, to quantify shorter-wavelength surface loading and to detect slow slip events and magma movement (Fig. 8).
As the first step to confirm the validity of fiber-linked optical clocks for crustal deformation monitoring, we plan to compare time-resolved chronometric leveling measured with a relatively short fiber link around Tokyo area (Akatsuka et al. 2020) with those obtained by GNSS. As a next step, the fiber link will be extended to Mizusawa VLBI (Very Long Baseline Interferometry) Observatory of the National Astronomical Observatory of Japan, located in the Tohoku region (Fig. 4). Using this baseline with approximately 400 km, we plan to observe the ongoing post-seismic deformation, which was triggered by the 2011 M9 earthquake, to validate long-distance chronometric leveling (Fig. 4 (Lower)). Such chronometric leveling will also reveal short-term and seasonal potential changes. Based on the observed potential changes, an uncertainty of atmospheric delay models used in GNSS analyses and the effects due to tides and surface loads can be investigated.
We expect that presenting observational data obtained with state-of-the-art optical lattice clocks deployed in a plate subduction zone and proposing models to interpret the temporal variations in the gravitational potential would contribute to the seismic and volcanic hazards mitigation studies in the future.
Author contributions YT and HK designed research. YT performed GNSS data analyses and calculated gravitational potential changes due to geophysical phenomena.

Data availability
The GNSS datasets analyzed during the current study are available from the web site of the GSI (http:// datah ouse1. gsi. go. jp/ terras/ terras_ engli sh. html). All data generated during the current study reported in this article are available on reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated Fig. 8 Schematic representation of the crustal deformation monitoring using fiber-linked optical lattice clocks. Vertical motion is masked by atmospheric noise in conventional GNSS observation. By combining chronometric leveling with GNSS observation, we can detect vertical deformation more quickly and precisely. Improved crustal deformation data with significantly reduced atmospheric noise allows the fault slip on the plate subduction interface to be inferred more clearly Page 13 of 14 93 otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.