Estimation of temporal and spatial variation of sound speed in ocean from GNSS-A measurements for observation using moored buoy

Using Global Navigation Satellite System–Acoustic (GNSS-A) technique, we have been developing observation system on a moored buoy for continuous monitoring of seafloor crustal deformation. The sound speed structure near a warm current has heterogeneity, which is the main cause of a seafloor positioning error. Assuming a sloping structure, previous studies proposed sound speed model to reduce positioning error. We examined the validity of the model by comparing the estimated structure with the actual structure measured at multiple points around our observation site. The result shows that the gradient parameter estimated from GNSS-A data acquired by vessel is appropriate. The numerical examination indicates that modeling error caused by the misinterpretation of the depth of gradient layer occurs, and it can be suppressed by performing acoustic ranging at the point near the centroid of units. From the calculation of estimation error of sound speed variation, the predicted acoustic ranging error observed using the moored buoy staying near the centroid is 9.0 cm or below. Therefore, seafloor displacement can be detected with centimeter class via moored buoy in the basin of a warm current.


Introduction
originally presented observation method for seafloor crustal deformation by combining Global Navigation Satellite System (GNSS) positioning with acoustic ranging. It is difficult to estimate the absolute position on the seafloor using GNSS only because the electromagnetic wave propagates poorly in seawater. GNSS-Acoustic (GNSS-A) technique uses a vessel or buoy as a relay point and estimates the relative position on the seafloor by acoustic ranging. Although the estimation accuracy for stationary displacement has been enhanced by the previous studies (Fujita et al. 2006;Honsho et al. 2019;Ikuta et al. 2008;Kido et al. 2008;Yasuda et al. 2017;Yokota et al. 2019), it is still difficult to detect that for non-stationary displacement or separate with high-accuracy coseismic from postseismic one. Yokota and Ishikawa (2019b) conducted an observation at frequent intervals and were successful to detect the displacement due to the slow slip event. From the results of this study, continuous monitoring is expected to increase the estimation accuracy. Recently, a buoy has been used as the continuous observation platform. Kido et al. (2018) developed the GNSS-A observation system on the moored buoy and had been operating it for almost 1 year. Kato et al. (2005) has utilized a buoy, named GNSS buoy, to observe tsunami. GNSS buoy was successful to detect tsunamis generated by the 2001 Peru earthquake, the 2004 off the Kii peninsula earthquake in Japan, and so on. Our research group then conducted experiments farther from the coast and extended the function of the GNSS buoy by introducing GNSS-A observation equipment. We previously reported the development of a comprehensive observation system for tsunami, ionospheric total electron content, precipitable water, and seafloor crustal deformation (Kato et al. 2018).
The moored buoy we used in this research is owned by Kochi Prefecture, Japan, which is located approximately 30 km off the coast of cape Ashizuri (Fig. 1). We installed three acoustic units on the seafloor around the buoy in 2017 and finished installing observation equipment, such as GNSS antenna and receiver, gyro sensor for attitude measuring, acoustic ranging equipment, computers for data processing, and satellite communication unit for data transfer, on the deck in 2018. The measurements of GNSS and gyro are acquired with 1 Hz sampling. The acoustic ranging at 3 min intervals for units one by one in sequence is continuously operated.
As shown by Yokota et al. (2019), from 2014 to 2016, water temperature at the surface layer around this region had gradient toward the southeast. The temperature gradient was generated by warm current called Kuroshio that flows northward in the Pacific Ocean along the eastern coast of Japan. Yokota end Ishikawa (2019a) also pointed out that sound speed structure has strong gradient when the observation point is located at the edge of the Kuroshio. Kido (2007) presented sound speed model considering horizontal gradient because the heterogeneity in sound speed structure is the main cause of positioning error for GNSS-A technique. Assuming a sloping structure for sound speed to estimate the variation from reference sound speed profile, Yasuda et al. (2017), Honsho et al. (2019), and Yokota et al. (2019) developed an analytical method. They eventually demonstrated that the estimation for the time series of unit position was improved.
In our case, we needed to use the measured in the past for analysis because our buoy was not equipped with the function of measuring sound speed. Therefore, the estimation accuracy for sound speed variation has more direct effects on the positioning accuracy. Accordingly, this study aimed to evaluate estimation accuracy for its temporal and spatial variations. We then measured sound speed at multiple points around seafloor units in parallel with acoustic ranging and obtained its three-dimensional structure. We applied the sloping structure model to GNSS-A data acquired using vessel and obtained seafloor unit position and sound speed variations. Then, we compared the estimated horizontal gradient with the measured structure. We also calculated the anticipated ranging error observed using buoy to examine the capability of the buoy to detect the seafloor displacement. First, we presented sound speed model and analytical method considering its sloping structure. Then, we showed the analytical results of array shape and horizontal gradient of sound speed by using GNSS-A data acquired using vessel. We discussed the estimated acoustic ranging error observed using buoy. The modeling error caused by misinterpretation of sound speed structure is also discussed.

GNSS-A observation using vessel
We conducted observation using research vessel Yugemaru from 2017 to 2019 to determine the initial position of seafloor units and survey the sound speed structure at the site. Figure 2a shows the track maps for six epochs. First, the vessel ran clockwise and counterclockwise around the unit array on the circular route in either case. The vessel also navigated along straight orbits passing both the vicinity of the center of the circular route and the point right above each seafloor unit. The acoustic ranging was performed at intervals of 8 s for units one by one in sequence for about 5-6 h for each epoch. During acoustic ranging, the onboard GNSS receiver recorded satellite signal and the gyrocompass measured attitude of the vessel. The position of onboard transducer can be calculated from the measured attitude and the position of GNSS antenna placed on the vessel estimated using kinematic positioning. We also measured conductivity, temperature, and depth (CTD) of sea water in parallel with acoustic ranging.

Sound speed structure obtained from CTD measurements
We conducted underway CTD observations at the points as indicated by the star symbols on the circular route in Fig. 2a. We performed at multiple points to understand sound speed structure at our observation site. Figure 2b shows sound speed profiles at 1 m intervals of water depth for each epoch, which is calculated using the equation proposed by Del Grosso (1974). The similar variation patterns according to the depth were measured The horizontal structure of sound speed has been generally considered to be sloping in the region near the Kuroshio basin (Kido 2007;Yasuda et al. 2017;Yokota et al. 2019;Honsho et al. 2019). We applied this model to detect horizontal gradient of sound speed from CTD measurements. Then, sound speed at the point (x, y, z) is modeled by where V 0 (z) (m/s) is sound speed at the origin and δV x and δV y (m/s/m) are the constant gradients of east and south components, respectively. Using the least squares method, the parameters V 0 (z), δV x (z), and δV y (z) at 1 m intervals of the depth from multiple profiles (Fig. 2b) can be obtained. Figure 2c shows the perturbation from V 0 (z) calculated by where R (m) is radius of circular route and ϕ (radian) is the direction of gradient axis from the true north. We considered that if the estimation errors of δV x (z) and δV y (z) are extremely large, the horizontal structure at a particular depth is not sloping. Then, the perturbation was not painted at the depth where the radius of error ellipse for ΔV(z, ϕ) with 95% confidence probability was higher than its estimation value.
The results in 2017 and 2019 showed that a strong gradient mainly from southeast to south is generated at the depth 200-400 m. Contrarily, the structures in 2018 have weak gradient or complicated gradient field. Yokota and Ishikawa (2019a) categorized a structure in accordance with the relative position between observation site and the Kuroshio basin. The structure in shallow depth has weak gradient in the area inside the Kuroshio. On the edge of the Kuroshio, there is strong thick gradient layer. The unstable disturbance is likely to be generated because there is no strong current flow in the area outside the Kuroshio. According to this categorization, our observation site was considered to be on the edge for the epochs in 2017 and 2019 and outside the Kuroshio for the epochs in 2018.

Model definition for travel time of acoustic signal
Travel time between onboard transducer x s = (x s , y s , z s ) and seafloor unit x u = (x u , y u , z u ) is defined as the integration of slowness S, which is an inverse of sound speed, along the propagation path: where L = |x s − x u | is path length (m). To express the equation of travel time simply, we used S instead of sound speed. As shown in Eq. (1), several previous works considered that horizontal sound speed structure is sloping. We also assumed that horizontal structure consists of N layers, the thickness of each layer is H (m), and the gradient is uniform within each layer. Then, the slowness at the depth in nth layer at the time t is modeled by where S 0 (z) is slowness (s/m) at the origin and∇S n = ( δS x,n , δS y, n ) are the east and north components of horizontal gradient (s/m/m) in nth layer, respectively. According to Ikuta et al. (2008), the temporal variation of slowness a(t) is represented by a cubic B-spline curve: where α i is control points and B i,4 (t) is basis function. We set 3 h for the knot interval in this study to detect long-term variation caused by the ocean tide. Assuming that acoustic signal propagates linearly as described in Fig. 3a, we obtained dl = L/D dz (D = z s − z u ). By substituting Eq. (4) into (3), the theoretical travel time can be written as where the average sound speed at the origin is expressed as and the contribution of horizontal gradient to travel time is By substituting the coordinate of the point on the path (8), the east component of ΔT grad is given by Kinugasa et al. Progress in Earth and Planetary Science (2020) The contributions of gradient to travel time are determined by its strength and thickness and depth of gradient layer. As suggested in Eq. (9), the parameters of strength δS x,n and thickness H are inseparable. Therefore, the previous works generally set H to 500-1000 m (Yasuda et al. (2017), Honsho et al. (2019)). Since the depth of our observation site is below 800 m ( Fig. 1), this study assumed that strength and direction of gradient are constant from sea bottom to sea surface (Fig. 3b). Then, we obtained by substituting N = 1, H = D, and ∇S n = ∇S = (δS x , δS y ) into Eqs. (8) and (9). ΔT grad indicates a sum of product of gradient value and horizontal distance from the origin to the middle of path. Accordingly, theoretical travel time model applied in this study is obtained from Eqs. (6) and (10).

Analytical process for detecting seafloor displacement
We need to determine the shape of seafloor unit array in advance for detecting precisely seafloor displacement observed using buoy. Then, first, we determine the array shape using data acquired using vessel, and second, we estimate the displacement of array from the initial position with assumptions that every single unit moves in the same direction and also at the same speed due to the crustal deformation. The theoretical travel time from the transducer of the vessel to the unit i for kth shot at time t k is given by The position of unit i for mth epoch is expressed as x u,i + δx m , where x u,i is its initial position and δx m is the vector of array displacement. To determine the array shape, the parameters to be estimated are x u,i , δx m , and temporal variation a(t k ) and horizontal gradient ∇S m from the slowness profile S 0,m . The average of their profiles was applied for S 0,m , since this analysis is performed without the information of horizontal sound speed structure obtained from multiple profiles (Fig. 2c).
As shown in Fig. 2c, the horizontal structure changes day by day. Although the time variation of ∇S should be considered, we treat ∇S as constant value during a certain time in accordance with Honsho et al. (2019). We assumed that ∇S is constant during observation for single epoch since the observation using vessel takes approximately 6 h for this site. For the observation using buoy, we were planning to divide observation period by a prescribed time interval and to consider ∇S constant within interval.
The theoretical travel time from the buoy is also given by T cal (x u,i + δx p , x s,k , S 0 , a, ∇S q , t k ). δx p and ∇S q are those of the pth and qth divisions when observation period is divided by a prescribed time interval for each. For the determination of array displacement, the unknown parameters are δx p , a(t k ), and ∇S q . We need to use the profile measured by vessel in the past because our observation system on the buoy has no function for CTD measuring. Both analyses use the least squares method and determine the unknown parameters that minimize the sum of residuals between observed and theoretical travel time.

Array shape of seafloor units
First, we examined the effect of considering horizontal gradient of sound speed for array shape determination. For this purpose, we compared array shape determined using sound speed model considering and ignoring gradient. Here, the model considering gradient (Eq. (4)) is called constant gradient model, and model ignoring gradient is called stratification model. By conducting analysis for each single epoch, a set of six shapes was obtained. Figure 4 compares the obtained array shapes for each sound speed model. The centroid position of each shape is arranged to the origin. For both horizontal and vertical shapes, the dispersion was suppressed by considering gradient. Table 1 also compares averaged distance between units with the standard deviation. The maximum dispersion of horizontal one between unit 1 and unit 3 was decreased from approximately 10 to 7 cm by considering gradient. The dispersions of vertical one were also decreased. We also showed the positional shift of unit array between two models for each epoch in Fig. 5. Violet and blue vectors in Fig. 5a denote the estimated horizontal gradient and the centroid shift, respectively. Two arrows  point in almost the same direction. Since violet arrow points in the direction of increasing sound speed, the estimated unit position is shifted toward faster sound speed by ignoring gradient. Moreover, the shift amount seems to be proportional to the quantity of gradient parameter. Vertical array shift in Fig. 5b suggests that the depth of unit is estimated shallower than the true position for the unit where the sound speed is faster. Second, we conducted analysis to determine common array shape for six epochs by using constant gradient model. Table 2   north, longitude 133.20850 east, and at ellipsoidal height − 756.15 m). From the calculated major axial radius of error ellipse with 95% confidence probability, the unit position was determined with 1.5 cm accuracy.

Horizontal sound speed structure
We also obtained horizontal gradient from the analysis to determine common array shape in the above section. The estimation results are summarized in the left column in Table 3. The magnitude of horizontal gradient ΔV (m/s/km) was converted from ∇S by ΔV ≈ j∇Sj=S 2 0 Â 10 3 . The estimated direction of gradient axis ϕ is also displayed. The estimated ϕs of epochs (1) and (2) in 2017 are southeast, which corresponds with the measured structure (Fig. 2c). The estimated ΔV value of epoch (3) is about half those of epoch (1). The effect of gradient on the observed travel time was canceled because the measured structure of epoch (3) is complicated. The estimated ΔV is smallest for epoch (4) because the structure has weak gradient. Although the structure of epoch (6) has large gradient at the depth 200-300 m, ΔV is smaller than those of epochs (1) and (2). It suggests that the effect of gradient is suppressed because the thickness of gradient layer is thin.

Estimation error of sound speed variations
A comparison with actual structure demonstrated in the previous section that horizontal gradient estimated from GNSS-A data with sound speed profile obtained in parallel with acoustic ranging was appropriate. Although our buoy has no function for measuring sound speed, the profile measured in the past needs to be used. Accordingly, we examined the effect of substituting the profile obtained on other days on the estimation of sound speed variations. We applied the profiles measured on 6 June 2017, and 5 June 2019, to the calculation for other epochs and obtained temporal variation and horizontal gradient from those reference profiles. Table 3 compares the horizontal gradient estimated using different profiles. The root-mean-square (rms) errors of residual between observed and calculated travel time are also shown. The rms error was 0.04-0.06 ms in the case of using profile obtained in each epoch. We then considered estimation using other profile to be validated when the rms error is less than 0.06 ms. Based on this index, the profile of epoch (1) should be selected for epochs (2)-(4) and that of epoch (6) should be selected for epoch (5). As shown in Fig. 2b, observed variation patterns along the depth for two epochs in 2019 were quite different from others. Thus, applying the profile markedly different from the actual pattern has adverse effect on the estimation. Even if the suitable one was selected as the reference, the estimation errors for ΔV and ϕ are maximum of 0.03 m/s/km and 3°, respectively. Figure 6 also compares the estimated temporal variation of averaged sound speed (1=aðtÞS 0 ). As mentioned above, the profile of epoch (6) should not be selected as the reference for epoch (1)-(4). However, regardless of the used profile, the temporal variation patterns were successfully identified. The offset in the figures represents the average of the residual of estimation values between used profiles. The results of all epochs indicate that using suitable profile can suppress the offset and the maximum value is 0.11 m/s.

Discussion
Separability between array displacement and sound speed gradient Honsho et al. (2019) discussed the requirements to detect array displacement separately from horizontal gradient. In two-dimensional example, they theoretically demonstrated that the measurements obtained from the point survey at the center of an equilateral triangle array have little information to separate them. The buoy we applied is in a similar situation as shown in Fig. 7 showing its movement for a month observed through GNSS. The position of buoy had been stable within a square of approximately 400 m × 400 m. Therefore, we also considered the separability with threedimensional example for observation using buoy.  According to Eq. (11), the observed travel time for unit i for kth shot is given by where L i,k = |x s,k − (x u,i + δx)|. The quantity of travel time changes caused by the array displacement δx becomes where L 0 i;k ¼ jx s;k −x u;i j. The temporal variation of S 0 is ignored here for brevity. The contribution of horizontal gradient to travel time is The position of onboard transducer x s,k is fixed at the origin if the buoy stops at the centroid of unit array. Assuming that the array shape is equilateral triangle and all units are at the same depth, the path length L and the incident angle of acoustic signal θ are independent of unit and shot (Fig. 8c). Then, the quantities of travel time changes in Eqs. (13) and (14) are rewritten as and where the directions of array displacement and slowness gradient are γ and ϕ S and x u,i is expressed by azimuth ψ i (Figs. 8a, b, d). From these equations, the separation of two quantities is theoretically impossible when the buoy stops at the centroid.
Our buoy stays stable but actually moves 50-100 m a day. To show the separability in the actual condition, we performed the numerical calculation assuming that the buoy moves around the centroid. Synthetic travel time for 1 year was generated from Eq. (12) by setting array displacement (δx, δy) = (−44.2, 35.7) mm/year and horizontal gradient (|∇S| = 0.2 ns/m/m, ϕ S = γ = 309°). We examined the effect of the movement range of x s,k by giving the random value within squares of 10 m × 10 m and 50 m × 50 m according to the uniform distribution. Figure 9 compares the results of δx and ∇S estimated once a week and every 12 h, respectively. A comparison between two ranges of movement shows that the estimation accuracy becomes higher with wider movement range. Therefore, the actual condition for our buoy meets the requirement for determining array displacement correctly.

Acoustic ranging error caused by estimation error of sound speed
We examined the effect of substituting profile obtained on other days for the estimation of sound speed variation (Table 3 and Fig. 6). To examine the capability of buoy to detect seafloor displacement, we discussed the estimated acoustic ranging error for observation using buoy in this section. From the positional relationship between buoy and seafloor units (Fig. 7), the average of observed travel time can be calculated and summarized in Table 4.
First, we consider ranging error due to the misinterpretation of sound speed average. As indicated in Fig. 6, the estimation error of averaged sound speed ϵ V was maximum of 0.11 m/s. Since the maximum of T obs is 0.82 s for unit 3 (Table 4), the maximum ranging error given by ϵ V T obs is 9.0 cm. As indicated in Table 3, estimation error of horizontal gradient ϵ ΔV was less than 0.03 m/s/km. ϵ V caused by ϵ ΔV is also given by ϵ V = x grad ϵ ΔV , where x grad is horizontal distance from the origin to the middle of path along the gradient axis. Since the maximum of x grad is approximately 521 m for unit 1 (Table 4), the ranging error due to the estimation error of horizontal gradient becomes 1.0 cm as a maximum.
We demonstrated that acoustic ranging from buoy without CTD measuring parallelly will probably be performed with under 10 cm accuracy. Consequently, proposed analytical method can achieve seafloor positioning with accuracy of several centimeters class, which is sufficient for detecting seafloor displacement.

Validity of proposed sound speed model
The quantities of contribution of horizontal gradient to travel time is determined by its strength, thickness, and depth of gradient layer, as well as signal route passing through the layer, as indicated in Eq. (8). Particularly, the depth of layer determines the ratio of the contributions between x u and x s into the change of travel time.
Assuming that sound speed structure has constant gradient from seafloor to sea surface, the ratio is fixed to 1: 1, as shown by the coefficients of x u and x s in Eq. (10). However, the actual structures are not so simple (Fig.  2c). Then, modeling error can occur, which can cause seafloor positioning error.
In this section, we demonstrated the contribution of horizontal gradient to observed travel time by numerical calculation and verify the validity of the assumption for sound speed model. We examined five types of structure shown in the left column in Fig. 10. The contribution of gradient ΔT grad was calculated for each structure by using Eq. (8) assuming that acoustic ranging is performed from the point in the displayed area. The distribution of ΔT grad for each unit is described by the color scale.
Type A is the structure having constant gradient from sea surface to seafloor, which is the same as the introduced model in this study. The distribution of ΔT grad suggests that if the gradient is ignored, the estimated position for every single unit can be shifted into the direction of gradient axis. Additionally, the unit where sound speed is faster can be estimated to be shallower. This result completely corresponds to that in Fig. 5.
The structure of type B has gradient layer in the middle depth and the thickness of layer is set to 100 m. The distribution of ΔT grad is nearly the same as type A. Accordingly, even if the thickness of actual gradient layer is thinner than the assumed, the effect of gradient can be removed when the middle depth of actual layer corresponds to that of the assumed.
Type C has gradient layer in the depth of 200-300 m. In this case, the signal propagates through the layer near from the ranging position x s . Then, the contribution ratio of x u becomes small. Contrarily, if the structure has gradient layer near the seafloor, ΔT grad is almost independent of x s . Type D simply reproduces the measured structure of epoch (3), which consists of two layers. The distribution of ΔT grad shows that major effect of gradient is cancelled because the direction of gradient axis is opposite. For this reason, the estimation value of ΔV was small (Table 3). Lastly, type E reproducing the structure of epoch (6) also consists of two layers, and the gradient of upper layer is much stronger. The effect of strong layer is dominant, and the distribution of ΔT grad is similar to that of type C because the depth of layer is shallower than middle.
The numerical calculation suggests that modeling error occurs if the true depth of gradient layer is shallower or deeper than that of the assumed, such as types C and E. It makes the error increase at the point far from the origin. If the buoy stays stable near the array centroid, the modeling error hardly occurs because the contribution of x s can be ignored. In this case, the estimated value as the gradient parameter includes the information of depth and thickness of layer. Therefore, the assumption for sound speed structure in this study is basically appropriate for observation near the array centroid. If the buoy stays stable at the point far from the centroid, the modeling error could be proven fatal. Yokota and Ishikawa (2019a) proposed estimating the contributions of x u and x s to ΔT grad separately. This method can reduce the modeling error irrespective of the ranging position. If we intended to introduce their idea for a point survey, we would be faced with another problem:  insufficient number of observation data because we have only three units. Accordingly, the depth of gradient layer is required to be determined by using the Monte Carlo method or some other calculation techniques in the case of point survey far from the array centroid.

Conclusions
We have been developing GNSS-A observation system on the moored buoy to monitor seafloor crustal deformation continuously. The sound speed structure in the near a warm current has heterogeneity that is main cause of seafloor positioning error. Some previous works proposed sound speed model assuming a sloping structure. We applied this model to our analysis using data acquired through vessel and obtained horizontal gradient of sound speed with seafloor unit position. To validate the estimated gradient, we measured sound speed at multiple points around the unit array and obtained the three-dimensional structure. A comparison with measured structure showed that the estimated gradient was appropriate.
Since the buoy we applied in this research is not equipped with the function of measuring sound speed, we needed to use the profile measured in the past. Then, we also compared the estimation result for sound speed variations using profile observed on other days. The maximum estimation error for averaged sound speed and magnitude of horizontal gradient were 0.11 m/s and 0.03 m/s/km, respectively. Therefore, the estimated ranging error for observation using buoy will be less than 9.0 cm. This suggests that the presented analytical method can detect seafloor displacement within a required accuracy.
We indicated that it is difficult to separate array displacement from horizontal gradient of sound speed when the buoy completely stops at the unit array centroid. We also suggested that the misinterpretation of the gradient layer depth causes modeling error, and it possibly affects array displacement determination when the moored buoy stays stable at the point far from the array centroid.