A new GNSS-acoustic positioning software implementing multiple positioning functions considering nadir total delays

Global Navigation Satellite System–acoustic (GNSS-A) positioning is an important geodetic observation technique for detecting sea�oor crustal deformation. After the 2011 Tohoku-oki earthquake, GNSS-A observational networks were extended along various subduction zones, and observational systems have been improved, especially for sea-surface platforms, such as the introduction of an unmanned vehicle, the Wave Glider. The aforementioned development of GNSS-A observations has provided a large amount of observational data. Furthermore, GNSS-A positioning methods were recently developed considering the lateral heterogeneity of the sound speed structure. Thus, it is important to develop a software that makes it easy for widespread use of the latest GNSS-A positioning methods. However, there is currently only one open-source GNSS-A positioning software, which may hinder the entry of various researchers into GNSS-A positioning analyses. Here, we developed a new GNSS-A positioning software, henceforth called “SeaGap” (Software of enhanced analyses for GNSS-acoustic positioning), that executes various positioning methods from the conventional kinematic positioning technique to the latest Markov Chain Monte Carlo (MCMC)-based static positioning technique. We introduce their methodology and demonstrate its application to actual observational data. The software newly added optional prior distributions to the unknown parameters expressing the heterogeneity of a sound speed structure into the MCMC-based static positioning method, and we also applied the new method to actual observational data. In addition to the positioning functions, the software contains various auxiliary functions, including drawing. The developed software is written using the “Julia” language and is distributed as an open-source software.


Graphical Abstract
Introduction Spiess (1985) contrived the global navigation satellite system-acoustic (GNSS-A) positioning technique, which combines GNSS positioning on a sea surface platform and acoustic ranging (round-trip travel time measurement) between the sea surface platform and seafloor transponders.One GNSS-A site is composed of multiple (from three to six) seafloor transponders that form a triangular or square-shaped array with the array diameter of nearly its water depth.The GNSS-A positioning analysis is composed of three processes: (1) Positioning a kinematic GNSS antenna attached to a sea surface platform and transforming the GNSS antenna position in a local ENU coordinate, (2) Transforming the GNSS antenna positions into the sea surface transducer positions considering attitudes of the sea surface platform and the sea surface transducer position relative to the GNSS antenna (hereafter, ATD (Antenna TransDucer) offset), (3) Positioning the seafloor transponders by minimizing travel time residuals between the observed travel times and the modeled travel times using the acoustic ranging data.Although we can obtain an absolute position of an individual transponder in the process (3), "array displacement" has been often estimated to obtain seafloor crustal deformation.The array displacement is generally estimated by following step; firstly, individual seafloor transponder positions are solved (i.e., individual transponder positioning) and are fixed as initial seafloor transponder positions.Then, a displacement of the seafloor transponder array relative to the initial seafloor transponder positions (i.e., translation of the seafloor transponders; it corresponds to the array displacement) is solved (i.e., array positioning, Fig. 1).From the positioning results of the array displacement, crustal deformations and plate motions are discussed.Spiess et al. (1998) first reported the geophysical results of GNSS-A observations of offshore inter-seismic displacements along the Cascadia subduction zone.After the success of Spiess et al. (1998), the GNSS-A observations revealed crustal deformation due to a seismic cycle: inter-seismic deformation [e.g., Peru-Chile trench (Gagnon et al. 2005), Nankai trough (Yokota et al. 2016;Yasuda et al. 2017), Japan trench (Fujita et al. 2006;Sato et al. 2013a)], coseismic deformation [e.g., the 2004 M w 7.5 Kii Peninsula earthquake (Kido et al. 2006;Tadokoro et al. 2006), the 2011 M w 9.0 Tohoku-Oki earthquake (Sato et al. 2011;Kido et al. 2011)], and post-seismic deformation [e.g., the 2011 Tohoku-Oki earthquake (Tomita et al. 2017;Watanabe et al. 2021)].Furthermore, GNSS-A observations have also contributed to the detection of slow-slip events (e.g., Honsho et al. 2019;Yokota and Ishikawa 2020) and oceanographical studies (e.g., Yokota et al. 2020).GNSS-A observational systems are significantly improved since the first introduction, and one of the most important developments is the evolution of sea surface platforms.A side-mounted pole with a transponder attached to a research vessel (e.g., Fujita et al. 2006;Ikuta et al. 2008) or towing buoys (e.g., Kido et al. 2006) was typically used in the 2000s.The introduction of a research vessel with a hull-mounted transducer enabled to effectively obtain acoustic ranging data with low noise and elaborate geometric ship tracks (e.g., Sato et al. 2013b).A Wave Glider, which is an unmanned surface vehicle, successfully obtained GNSS-A observational data (e.g., Iinuma et al. 2021), and it has started to collect great amounts of the observational data.Furthermore, the extension of the GNSS-A observational networks is an important recent advancement.The GNSS-A observational results of the 2011 Tohoku-Oki earthquake demonstrate the importance of seafloor geodesy.After the event, GNSS-A observational networks largely extended along the Japan Trench (Kido et al. 2015) and the Nankai Trough (Yokota et al. 2016).Such extensions of the GNSS-A observational networks were conducted not only around Japan, but also in other regions [e.g., Taiwan (Chen et al. 2021), Cascadia (DeSanto et al. 2022), and Alaska (Brooks et al. 2023)].
Compared to the enrichment of GNSS-A observational networks and the accumulation of GNSS-A observational data, only one GNSS-A positioning software is publicly distributed at this moment: the "GARPOS" software which has the original version (Watanabe et al. 2020) and the updated version, "GARPOS-MCMC" (Watanabe et al. 2023).The shortage of GNSS-A analysis tools may hinder the entry of many researchers into analyzing GNSS-A observational data and developing GNSS-A observational studies.Furthermore, GARPOS only implements only a static positioning method that estimates the three-dimensional transponders' positions (either of individual or translational positions) using a large amount of acoustic ping data.The development of a kinematic GNSS-A positioning method, which estimates the array displacement using a small amount of acoustic ranging data, is also important for handling limited amounts of the observational data due to ship-time or for conducting real-time GNSS-A positioning (e.g., Imano et al. 2019).
In this study, a new tool for GNSS-A positioning was developed called SeaGap (software of enhanced analyses for GNSS-acoustic positioning) using a programming language "Julia", which has high readability and computational speed (Bezanson et al. 2017).Although GARPOS provides a single computational function, we can flexibly assign unknown parameters and perform various types of the GNSS-A positioning analyses such as individual transducer positioning, array positioning.On the other hand, the system of SeaGap is divided according to the analysis conditions [e.g., calculation method, kinematic/ static positioning mode, individual transponder/array positioning mode, and assumption of a sound speed structure (SSS)], and we call each divided analysis system as a "function" of SeaGap in this study.SeaGap prepares many computational functions not only for performing the various types of GNSS-A positioning analyses but also conducting useful post-processing, such as visualization and statistical analysis.In this study, the methodology, the GNSS-A positioning functions in SeaGap, and their applications are introduced.Details on the use of SeaGap are provided in the online manual.

Variety of GNSS-A positioning schemes
Since the GNSS-acoustic observation was contrived (Spiess 1985), the following two-step procedure has been often adopted; we first perform the individual transponder positioning method to obtain initial values for individual transponder positions and then perform the array positioning method to obtain array displacements as a seafloor displacement.Although this two-step procedure is important for obtaining precise seafloor crustal deformation (e.g., Kido et al. 2006;Watanabe et al. 2020), several studies did not adopt the two-step procedure (e.g., Fujita et al. 2006;Ikuta et al. 2008).These studies performed the individual transponder positioning method and then successfully detected a temporal change of the centroid of the estimated transponder positions as a seafloor displacement.Note that most recent studies have adopted the two-step procedure, but the two-step procedure assumes geometry of the seafloor transponder array is invariant with respect to time (i.e., rigid motion of the array).If this assumption is inadequate due to large earthquake or strong ground motion, the two-step procedure cannot be applied (e.g., Sato et al. 2011;Kido et al. 2011).
Based on the above individual transponder positioning and the array positioning methods, various GNSS-A positioning schemes have been developed during these twenty years.The positioning scheme is strongly related to the characteristics of the collected observational data.GNSS-A observational data were collected using two types of surveys: a fixed-point survey and moving survey.A sea surface platform keeps its position above the center of the seafloor transponder array for the fixed-point survey, whereas it moved around for the moving survey (Additional file 1: Fig. S1).A fixed-point survey is a more traditional survey (e.g., Spiess et al. 1998;Kido et al. 2006) and is specialized to detect horizontal array displacements (Kido 2007).Spiess et al. (1998) indicated that, as long as the horizontally stratified SSS (HS-SSS) was valid, influence of the sound speed fluctuation appeared to be the same degree in all transponders with the same incident angle and did not affect estimation on the horizontal array displacement in the case of the fixed-point survey.Moreover, the fixed-point survey data provide theoretically high sensitivity to detection of the horizontal array displacements (Kido 2007) even if the initial transponder positions are not well-determined.Due to these advantages, fixed-point survey data are suitable for a kinematic array positioning method that estimates the horizontal array displacement for each acoustic ping.However, horizontal array displacements may contain systematic biases when the spatial heterogeneity of a SSS is significant.Contrastingly, moving survey data have sensitivity to estimating vertical positions (e.g., Sato et al. 2013b) and the spatial heterogeneity of SSS (e.g., Yokota et al. 2018;Honsho et al. 2019).To exploit this advantage, a static array positioning method was used to estimate the three-dimensional array displacements and sound speed gradients using moving survey data.Although fixedpoint survey data with diversity of the incident angles are theoretically sensitive in estimating the vertical array displacement and spatial heterogeneity of SSS (Tomita et al. 2019), positioning accuracy of this type of the data has not been well-investigated.In addition, moving survey data were utilized to determine individual seafloor transponder positions as initial values using a static positioning scheme.

Underwater acoustic ranging models using nadir total delay
As indicated in Sect."Introduction", the seafloor transponder positions are eventually estimated by minimizing residuals between the observed and the modeled travel times.Considering the travel time of the n th shot number for k th transponder in the acoustic ranging data, the simplest observation equation is written as where T obs n,k and T cal are the observed and calculated round-trip travel times, respectively.The k th seafloor transponder position is p k , and the sea surface trans- ducer position is u which is calculated from the GNSS antenna position with the ATD offset b 0 .v x, y, z, t is a four-dimensional SSS; x , y , and z are the east-west (EW), north-south (NS), and depth positions (i.e., Eastward, Northward, Downward)  are the observational times when an acoustic signal is transmitted and received by the sea surface transducer, respectively.GNSS-A observation typically measures round-trip travel times by the mirror transponder system for canceling out various external effects such as due to ocean currents and a drift of the clock.The details of this acoustic ranging system are documented by Spiess et al. (1998) andFujimoto (2014).For the round-trip travel time calculation for T cal , outward and return one-way travel times are separately calculated, and then they are summed up.
Although it is plausible to calculate travel times under the complex four-dimensional SSS following Eq. 1, it is difficult to detect the complex underwater SSS properly.Moreover, the high computational cost of the travel time calculation in a four-dimensional SSS based on the eikonal equation also makes it difficult to solve Eq. 1 directly.Therefore, most previous studies typically have calculated travel times assuming an HS-SSS constructed from a reference sound speed profile v 0 (z) for fast computation (e.g., Chadwell and Sweeney 2010).In practice, v 0 (z) can be obtained by sound speed measurements (e.g., a CTD profiler) or statistical models (e.g., World Ocean Atlas (1) y, z, t , 2018;Garcia et al. 2019).Assuming an HS-SSS, Eq. 1 is written as where T n,k is a travel time delay due to the assump- tion of the HS-SSS.As T n,k often provides systematic error in the positioning results, many studies developed various modeling approaches to reduce T n,k .The con- ventional studies have modeled T n,k under an HS-SSS with temporal variation (e.g., Fujita et al. 2006;Kido et al. 2006;Ikuta et al. 2008), while the recent studies have included the effects of horizontal gradients in the models (e.g., Yokota et al. 2018;Honsho et al. 2019;Watanabe et al. 2020;Kinugasa et al. 2020).Note that we can express the observation equations for both the individual transponder positioning and the array positioning by formulation of p k as where p 0 k is the initial transponder positions, and δp k and δp are the unknown parameter vector for the indi- vidual transponder positioning and the array positioning, respectively.δp k is an individual deviation for the k th transponder from p 0 k , while δp is the common translation from p 0 k .Here, we defined slowness to describe the travel time delay as with Honsho et al. (2019).The slowness is a reciprocal of the sound speed, and then we defined the slowness based on v x, y, z, t and v 0 (z) as s x, y, z, t and s 0 (z) , respectively.Considering the perturbation compo- nent of the slowness from the reference sound speed profile as s ( s = s 0 + s), the travel time delay is expressed as where L n,k is the round-trip acoustic path and dl is derivative along the acoustic path.Under an HS-SSS with temporal variation, the perturbation component of the slowness depends on depth and time, meaning �s x, y, z, t ≈ �s H (z, t) .For simplicity, approximating the acoustic path as a straight line, derivative along the acoustic path is expressed as dl = 1 cosξ dz using the inci- dent angle ξ .Following this approximation, the travel time delay under an HS-SSS with temporal variation T H is written based on Eq. 4 as (2) where C(t) is the travel time delay integrated along the nadir direction of a sea surface platform and corresponds to the nadir total delay (NTD) defined by Kido et al. (2008) and Honsho and Kido (2017).Note that NTD could vary depending on the water depth of each seafloor transponder Z k ; however, as the travel time delay typically appears in the shallow underwater portion, we approximate that the influence due to the water depths of the seafloor transponders was insignificant and introduced NTD C(t) which was independent from Z k .The NTD is an analogous concept to the Zenith Total Delay in the GNSS positioning (e.g., Marini 1972) with a mapping function of 1 cosξ .Although complex mapping functions were often used in GNSS positioning, this simple mapping function has been used for GNSS-A positioning (Tomita et al. 2019).Using the concept of NTD, Eq. 2 under an HS-SSS with temporal variation can be rewritten as with where t n,k corresponds to 2 by assuming that the SSS does not significantly change in short time (< ~1 min between the transmission and the reception of an acoustic signal).Equation 6indicates the observed travel times are transformed with the mapping function M n,k and that this transformation practically provides a weight for each observational value depending on its acoustic path length.For the accurate estimation of NTDs, it is important to collect travel times from multiple transponders, because estimation of NTD is an ill-posed problem when only using a single seafloor transponder and cannot improve the positioning accuracy.
For the introduction of horizontal heterogeneity in a SSS, we modeled the slowness structure as superposition of the HS-SSS with temporal variation and the perturbation following Honsho et al. (2019): �s x, y, z, t ≈ �s H (z, t) + �s G x, y, z, t , where s G indicates the horizontal heterogeneity.Assuming the horizontal heterogeneity was expressed as a gradient at each depth layer, �s G x, y, z, t ≈ δ�s EW (z, t)x + δ�s NS (z, t)y where δ�s EW and δ�s NS are spatial gradients (i.e., spa- tial derivative of the slowness) in EW and NS directions, (6) respectively.In this case, the travel time delay due to the horizontal gradients T G n,k is expressed as with where δ�s G = δ�s EW , δ�s NS .x is decomposed into two factors of horizontal positions of the sea surface transponder u hor , and u hor n,k corresponds to u hor t n,k .The azimuth angle of the acoustic path is φ .Although Eq. 9 also approximates the acoustic path as a straight line as well as Eq. 4, influence of this approximation on the travel time delay calculation was confirmed to be small (Honsho et al. 2019).According to Eq. 8, G s • u hor and G d • h indicate travel time delays integrated along the nadir and the horizontal directions of the acoustic path, respectively, which are caused by the horizontal gradient.Moreover, G s • u hor and G d • h are sensitive to the shallow and deep sound speed layers, respectively (e.g., Yokota et al. 2018Yokota et al. , 2022;;Honsho et al. 2019).Thus, for convenience, we called G s and G d as the shallow gradient and the deep gradient in this study, respectively.Using Eqs. 3, 5 and 8, the observation equation in a SSS with horizontal gradients are written as The shallow and deep gradients, which are defined in Eqs. 10 and 11, practically consider a multiple-layered gradient structure (MLGS) by integrating the travel time delays due to the individual gradient layers (e.g., Yokota et al. 2022).However, to perform stabilize the estimation, approximation of a single-layered gradient (8) structure (SLGS) has been often utilized (e.g., Honsho et al. 2019;Kinugasa et al. 2020;Yokota et al. 2022).Honsho et al. (2019) assumed that single gradient layer existed from the sea surface to a certain depth D (here- after, referred as gradient depth) and derived the relationship between the gradients and gradient depth by applying the depth of D to Eqs. 10 and 11: According to this assumption, the spatial pattern of the deep gradient is defined equally as that of a shallow gradient, and the intensity of the deep gradient is expressed by the gradient depth.Thus, if a shallow gradient is given, then the deep gradient can be expressed in terms of a single parameter, i.e., the gradient depth, which implies that the unknown parameters are reduced.

Travel-time calculation
To perform GNSS-A positioning techniques following the observation equations in Sect."Underwater acoustic ranging models using nadir total delay" (e.g., Eq. 2), it is important to calculate synthetic travel times T cal in an HS-SSS defined by a reference sound speed profile v 0 (z) accurately (plausibly, ~1 mm along the path line ≈ ~± 1 microsecond; e.g., Chadwell and Sweeney 2010).SeaGap executes two types of travel time calculation methods developed by Tomita and Kido (2022), which are named the "exact" and "approximate" travel time calculations.The details of these travel time calculations were already shown in Tomita and Kido (2022); therefore, an overview of these methods is given below in this paper.
For the "exact" travel time calculation, a synthetic travel time in an HS-SSS based on v 0 (z) is calculated by the shooting method considering the Snell's law in a spherical frame.To express the spherical frame, the Earth's radius is provided as a Gaussian mean radius depending on the site latitude, which considers ellipsoidal curvature of Earth.Although the dependency of the azimuth of an acoustic path should be considered in the spherical frame approximation as stated by Chadwell and Sweeney (2010), the "exact" travel time calculation does not consider this effect.The range error due to this effect is less than ~ 1 mm when the horizontal distance of the acoustic path is within a few kilometers (~2 mm for even the case within the horizontal distance of ~10 km), and mis-modeling in the array displacement due to this effect is less than sub-millimeters for a GNSS-A site with the water depth of 5 km (Chadwell and Sweeney 2010).Thus, influence of the azimuth of an acoustic path does not critical in typical cases of GNSS-A observations.Moreover, since consideration of the azimuth dependency makes it complex to optimize the "approximate" travel time to the (13) "exact" travel time (this optimization is addressed later), the "exact" travel time calculation does not consider the azimuth dependency.
To calculate the "approximate" travel time, the travel time in a spherical Earth without considering the bend of an acoustic path is first calculated (we call this type of travel time as "the spherical travel time").As the difference between the spherical travel time and the "exact" travel time mostly depends on the incident angle and the transponder height, the correction terms, which are expressed by the 8th order polynomial functions depending on them, are introduced.The correction terms can be optimized to the set of an input sound speed profile and the individual seafloor transponder depth in advance.Then, the approximate travel time, which is the summation of the spherical travel time and the correction terms, coincides with the "exact" travel time within ~7 × 10 -9 s (~0.01 mm).Because of the pre-calculation of the correction terms, the approximate travel times can be quickly calculated.The calculation for the GNSS-A positioning implemented in SeaGap adopt the "approximate" travel time calculation to perform fast computation even for a Markov Chain Monte Carlo (MCMC) technique.
Besides, SeaGap prepares functions to calculate both "exact" and "approximate" travel times for a given sound speed profile v 0 , a sea surface point u , and a seafloor point p .Thus, synthetic travel-time data sets can be eas- ily produced and utilized for numerical simulations.

GNSS-A positioning functions in SeaGap
Here, we introduced GNSS-A positioning functions equipped in SeaGap, which employ the observation equations in Sect."Underwater acoustic ranging models using nadir total delay" and the travel time calculation techniques in Sect."Travel-time calculation".These functions are summarized in Table 1.Each function has its own observation equation under different analysis conditions.The analysis conditions include following 4 major factors: the calculation method [a nonlinear least squares (NLLS) method (Gauss-Newton method) or an MCMC method], types of the unknown parameters (a kinematic array positioning method or static positioning methods), the assumed underwater SSS (HS-SSS, an HS-SSS with a MLGS, or an HS-SSS with an SLGS).

A kinematic array positioning function by NLLS
For the acoustic units developed in the Scripps Institute of Oceanography (Spiess et al. 1998) and Tohoku University (Fujimoto 2014), all seafloor transponders replied to an acoustic signal from the sea surface, so that travel times for all transponders were simultaneously collected by a single acoustic shot from the sea surface transponder.Contrastingly, the acoustic units of the Japan Coast Guard and Nagoya University collect travel time data by calling an individual seafloor transponder (Ishikawa et al. 2020;Kinugasa et al. 2020).The former system is quite useful for conducting kinematic array positioning, because we can estimate the horizontal array displacement using a shot group data set collected by a single acoustic shot; therefore, kinematic array positioning was developed by the Scrips Institute of Oceanography and Tohoku University.Tohoku University documented various important observation results using a kinematic array positioning method that considers the NTD (e.g., Kido et al. 2006;Tomita et al. 2015).
SeaGap provides a kinematic array positioning function that estimates horizontal array displacement for each acoustic shot group in an HS-SSS (function name: "kinematic_array").The acoustic shot groups can be set arbitrarily by providing flags in an input data file; thus, even if the travel times are not collected simultaneously, we can execute the kinematic positioning method.According to Eqs. 3 and 6, the observation equation for each acoustic shot group is as follows: where δp hor i and C i are the unknown parameters of this equation, which are horizontal array displacement from the initial transponder positions p 0 k and the NTD for the acoustic shot group, respectively.i indicates the shot number belonging to each acoustic shot group.Although p 0 k has three components (EW, NS, and UD), we fixed the vertical component of the array displacement vector to be zero: To solve this equation, the number of observations must be at least three.Therefore, ( 14) the kinematic array positioning method was performed when the number of travel-time data points was three or more for each acoustic shot group by default.Optionally, the minimum integer of the travel time data points can be chosen from integers of three or more when executing the kinematic array positioning method.Although the unknown parameters are solved by NLLS, the convergence criteria can be changed by users.Tomita et al. (2019) developed kinematic array positioning methods to estimate the three-dimensional array displacement for a multi-angled transponder site or to estimate using an extended Kalman filter; however, these functions are not included in the current version of SeaGap.

Static array positioning functions by NLLS
Static array positioning techniques have been mainly developed by the Japan Coast Guard (e.g., Watanabe et al. 2020), Nagoya University (e.g., Tadokoro et al. 2018), and Tohoku University (e.g., Honsho and Kido 2017;Tomita and Kido 2022).The conventional static array positioning methods (e.g., Watanabe et al. 2014) have estimated temporal evolution of the sound speed fluctuation in the HS-SSS as well as the array displacement.Based on the sensitivity of moving survey data to the spatial heterogeneity of a SSS, it has become popular to estimate sound speed gradients using a static array positioning method (e.g., Honsho et al. 2019;Watanabe et al. 2020;Kinugasa et al. 2020;Tomita and Kido 2022).SeaGap performs various static array positioning methods with and without considering sound speed gradients.
The static array positioning methods of SeaGap express the temporal fluctuation of the NTD using cubic B-spline functions as

Table 1 GNSS-A positioning functions implemented in SeaGap
"Mode" includes "Kinematic" (kinematic positioning), "Static" (static positioning), "array" (array positioning), and "individual" (individual transponder positioning)."Calculation Method" is NLLS (non-linear least-squares) or MCMC (Markov Chain Monte Carlo) methods."Structure" indicates the assumed sound speed structure, which is "HS-SSS" (horizontally stratified sound speed structure), "MLGS" (HS-SSS with multiple layered gradient structure), or "SLGS" (HS-SSS with single layered gradient structure) where j and c j are the j th cubic B-spline basis and its coefficient, respectively.The cubic B-spline bases were distributed at equal intervals, and the total number of cubic B-spline bases J controlled the roughness of the NTD fluctuation.Using Eqs. 3, 6 and 15, the observation equation for the static array positioning in an HS-SSS is written as The static array positioning function based on Eq. 16 is implemented in SeaGap (function name: "static_array").To execute this function, the users need to assign the total number of cubic B-spline bases J .A simple way to optimize the total number of cubic B-spline functions is to utilize an information criterion.SeaGap calculates the Akaike's information criterion (AIC) (Akaike 1974) and Bayesian information criterion (BIC) (Schwartz 1978) values, and easily optimizes the total number of cubic B-spline bases.The static array positioning function solves the three-dimensional array displacement δp and the coefficients of the cubic B-spline functions c by NLLS.
SeaGap also provides a static array positioning function simultaneously optimizing the position of the sea surface transducer relative to the GNSS antenna b 0 (function name: "static_array_TR").Its observation equation is the same as that in Eq. 16; however, b 0 is also estimated as an unknown parameter.Ideally, it is useful to solve this parameter simultaneously using all campaign observational data as long as the same sea surface platform was employed, regardless of the observation sites or the observation periods (Honsho and Kido 2017).However, such a type of the estimation complicates the program code.It is a simple but powerful method to estimate b 0 for the observational data for individual campaigns and then average the solutions of b 0 .
The static array positioning function under an HS-SSS considering MLGS (function name: "static_array_grad") is also provided.The horizontal gradients are modeled by two factors: the shallow and deep gradients as shown in Eq. 12.However, almost full contributions of the shallow gradient delay G s • u hor can be even explained by the temporal fluctuation of NTD C(t) which is the flex- ible cubic B-spline function, when the observational data (15 are obtained by a single sea surface platform (Honsho et al. 2019).This is because the shallow gradient delay expresses the travel time delay integrated along the nadir direction in the gradient structure (Eq.10).As the NTD is the travel time delay integrated along the nadir direction in an HS-SSS, there is a kind of trade-off relationship between the shallow gradient delay and the NTD; thereby, Honsho et al. (2019) called the shallow gradient as "NTD gradient" considering the similarity of the shallow gradient delay and the NTD.Then, it was found that the shallow gradients did not strongly affect the estimation of the array displacement when employing a single sea surface platform (Tomita and Kido 2022).Hence, "static_array_grad" eliminates the shallow gradients from the unknown parameters.Using this assumption and Eq. 15, the observation equation (Eq.12) is rewritten as The temporally constant gradient structure is assumed in Eq. 17 for the stable estimation ( G d (t) ≈ g d ).This observation equation was also introduced in Tomita and Kido (2022), and they examined its performance to the actual observational data.By contrast, GARPOS solves temporal fluctuations of shallow and deep gradients by providing constraints on their temporal fluctuations (Watanabe et al. 2020).Tomita and Kido (2022) demonstrated that the positioning method based on Eq. 17 was comparable to GARPOS in terms of the accuracy of the array positioning results when the observational data included the moving survey data obtained from geometrically well-distributed sea surface points.However, these methods are unsuitable for processing observational data including fixed-point survey data (e.g., Honsho et al. 2019).

An individual transponder positioning function by NLLS
Individual transponder positioning is similar to seismological hypocenter determination (e.g., Hirata and Masu'ura 1987) and is performed as a type of static positioning.To eliminate the temporal fluctuation of the average sound speed by precisely estimating the NTD, it is important to simultaneously estimate multiple seafloor transponder positions not to estimate single seafloor transponder position.Assuming an HS-SSS, the observation equation of the individual transponder positioning was written based on Eqs. 6 and 15: (17) SeaGap implements the individual transponder positioning function following Eq.18 as "static_individual".
For the accurate determination of δp k , it is effective to simultaneously employ multiple campaign data obtained at a certain site and then solve the individual transponder positions δp k and the array displacements among indi- vidual campaigns (e.g., Honsho and Kido 2017).However, SeaGap does not have such a sophisticated function to avoid complicating programs.Hence, it is effective to estimate δp k for each campaign datum and to average the solutions of δp k for the robust determination.

Static array positioning functions by MCMC
Using Eq. 12 and the condition of SLGS (Eq.13), Honsho et al. (2019) provided the observation equation under an HS-SSS with temporally constant SLGS: The gradients were assumed to be constant with time during the observation period, which indicated G d (t) ≈ g d and G s (t) ≈ g s , for stable estimation (e.g., Honsho et al. 2019;Tomita and Kido 2022).As mentioned above, the shallow gradient g s tends to have a tradeoff relationship with the temporal fluctuation of NTD.To avoid the trade-off relationship, Honsho et al. (2019) separated the NTD fluctuation C(t) into the long- term component (L-NTD) with the time scale of several hours or more and the short-term component (S-NTD) with the time scale of several minutes.They applied the polynomial functions and cubic B-Spline functions to L-NTD and S-NTD, respectively, then a smoothing constraint was imposed on temporal variation of S-NTD.Moreover, the gradient depth was fixed at a constant value a priori.The introduction of L-NTD enabled us to avoid the trade-off relationship and to estimate the shallow gradient.
Tomita and Kido (2022) also introduced L-NTD modeled by a quartic polynomial function to estimate the shallow gradient and also estimated the optimal gradient depth through an MCMC method.To avoid imposing the (18) smoothing constraint on the temporal variation of S-NTD, Tomita and Kido (2022) divided Eq. 19 into two observation equations: where γ indicates the coefficients of the quartic polyno- mial function.The unknown parameters in Eq. 20 are the array displacement δp , the cubic B-spline functions c , and the gradient depth D .Besides, those in Eq.21 are the shallow gradient g s and the coefficients for the quartic polynomial function γ expressing L-NTD.Although we call NTD expressed by the cubic B-spline functions c in Eq. 20 as S-NTD to distinguish it from L-NTD in Eq.21, S-NTD in Eq. 20 is different from S-NTD defined by Honsho et al. (2019) and is the combination of the NTD and the contributions due to the shallow gradients as same with Eq. 17. Optimizing Eqs.20 and 21 iteratively, we can estimate the shallow gradient g s avoiding the trade-off relationship with NTD in Eq.21 and then estimate the array displacement δp under the optimal deep gradient field D 2 g s which is composed of the shallow gra- dient (spatial pattern of the gradient in SLGS) g s and the gradient depth (intensity of the gradient) D.
Tomita and Kido ( 2022) developed an MCMC method to iteratively optimize Eqs.20 and 21.In the MCMC method, information regarding unknowns is expressed using a probability density function (PDF).According to Bayes' theorem (Bayes 1763), a posterior probability density for an unknown parameter vector x with a given data vector d (the travel time data) is expressed as where p(d|x) is the likelihood for a given x , and p(x) is the prior probability density of x .The likelihood func- tion for each observation equation (Eqs.20 and 21) was assumed to be a normal distribution of the travel time residuals with the standard deviation of 10 and 2 are the scaling factors for Eqs.20 and 21, respectively, which are included in the unknown vector x .The scaling factors were introduced to effectively optimize the observational error for each equation (e.g., Kubo et al. 2016;Tomita et al. 2021) and to balance the two observation equations (Eqs.20 and 21).Although the scale factors in Tomita and Kido (2022) were defined as the variance of the travel time residuals, this study defined the scale factors as the standard deviation to intuitively understand their values.Assuming the likelihood function for each observation equation follows a normal distribution, the likelihood is expressed as with where W i is the inverse of the variance-covariance matrix and is a diagonal matrix, and r i is a vector of the travel time residuals for Eq.20 ( i = 1 ) and Eq.21 ( i = 2 ).A uniform distribution with a sufficiently wide range was assumed for the prior PDFs of all parameters except the gradient depth D .For the gradient depth, a uniform distribution from the sea surface to the sea bottom (i.e., 0 ≤ D ≤ Water depth ) was assumed.
To sample x from the posterior PDF during the MCMC iterations, the random walk Metropolis-Hastings (M-H) algorithm (Metropolis et al. 1953;Hastings 1970) was employed.This algorithm adds perturbation to the current unknown parameters and then evaluate the perturbated unknown parameters to be accepted or rejected using an acceptance probability.For the m th MCMC iter- ation, the current state x (m−1) (the unknown parameter vector generated at the previous iteration) is perturbated, and the perturbated state x ′ (a candidate of the unknown parameter vector) are evaluated.If x ′ is accepted, the state is updated as x ′ ( x (m) = x ′ ); if x ′ is rejected, the cur- rent state is remained ( x (m) = x (m−1) ).The initial state x (0) indicates the initial values for the unknown param- eters.The l th component of the perturbed parameter is given by the random walk process as follows: where u l is a random value obtained from a uniform dis- tribution with range of [−0.5,0.5] and σ l is a step width for the l th unknown parameter.For the effective sam- pling, the unknown parameters in Eq. 20 are updated at ( 23) each odd MCMC iteration, whereas those in Eq.21 are updated at each even MCMC iteration.Thus, assuming that the unknown parameter vectors in Eqs.20 and 21 are expressed as x 1 and x 2 , respectively, x ′ is expressed as follows: Because the average of the S-NTDs in Eq. 20 generally shows a correlation with the vertical array displacement, it is occasionally useful to add a uniform perturbation to all the S-NTDs to promote the convergence of the solutions; then, assuming that the l S th component of the perturbed parameter corresponds to the parameter expressing S-NTD, the parameter update is given as follows: where u S−NTD and σ S−NTD are a random value obtained from a uniform distribution with range of [-0.5,0.5] and a step width, respectively, which are common among all S-NTD parameters for each MCMC iteration.This perturbation process for the S-NTD can be used optionally.
The acceptance probability at the m th MCMC iteration is written as follows: where q denotes the proposal probability density for gen- erating another state from the current state.As the proposal probability densities q x (m−1) |x ′ and q x ′ |x (m−1) are symmetric in this case, the acceptance probability is defined as the ratio of the posterior probability densities at the given states x (m−1) and x ′ .The acceptance prob- ability for the current state x (m−1) and the candidate state x ′ can be easily computed using Eqs. 22-24, 28.At each step, as clarified by introducing Eq. 26 into Eq.28, the acceptance probability in Eq. 28 corresponds to that for each block in a multiple-block M-H method (e.g., Chib 2001) when the unknown parameters are divided into two blocks as x 1 and x 2 .Although the original multiple- block M-H method counts one step when optimizing all blocks, this study counts one step when optimizing each block to monitor the convergence process for each block (i.e., each observation equation: Eqs.20 and 21) as shown ( 26) x ′ 1 , x 2 : at the odd step x 1 , x ′ 2 : at the even step . ( later in Sect."Static array positioning method using a MCMC technique".A random value u from a uniform distribution with a range of [0, 1] is produced for each MCMC iteration, and then accept the candidate when α x ′ |x > u .This MCMC method can be executed by the "static_array_mcmcgrad" function in SeaGap.Although we introduced the basic strategy of the MCMC method here, the detailed methodology and the flow-chart of the optimization process were shown in Tomita and Kido (2022).
SeaGap also provides an improved version of this function, called "static_array_mcmcgradc".As indicated by Tomita and Kido (2022), their method was not applicable to observational data without moving the survey data.The function "static_array_mcmcgradc" adds further constraints to the shallow gradient and the gradient depth by providing prior distributions.By introducing these prior distributions, the array displacement can be stably solved (see Sect. "Static array positioning method using a MCMC technique").
Regarding the shallow gradient parameters, a normal distribution with the mean of zero and the arbitrary standard deviation was given as a prior PDF.If the moving survey data are insufficient, it would be difficult to optimize the shallow gradient.Thus, when the data do not have sufficient sensitivity to constrain the shallow gradient, it should be zero.Tomita and Kido (2022) analyzed the 212 campaign data sets which were provided by Japan Coast Guard (6 sites: KAMN, KAMS, MYGI, MYGW, FUKU, CHOS) during 2011-2020 with elaborate geometric of sea surface tracks (Watanabe et al. 2021), and they estimated the EW and NS components of the temporally constant shallow gradient for each campaign data set.Additional file 1: Fig. S2 shows the histogram of the intensity of the shallow gradients combining both the EW and NS components for all campaign data sets.Considering the spread of the estimated shallow gradients in Additional file 1: Fig. S2, the default value of the standard deviation the "static_array_mcmcgradc" function was set to 0.15 ms/km, which could be changed optionally.
Regarding the gradient depth, a normal distribution is given as a prior PDF, because the posterior PDFs of the gradient depth are typically obtained as a normal distribution in Tomita and Kido (2022).Moreover, the range of the gradient depth is restricted from zero to the water depth by providing the prior probability density of zero when a candidate of the gradient depth is out of the range.The gradient depth was typically assumed (Yasuda et al. 2017;Honsho et al. 2019) or was estimated (Tomita and Kido 2022) to be ~0.5-1km.Additional file 1: Fig. S3 shows histograms of the gradient depths estimated by Tomita and Kido (2022) and their standard deviations using the 212 campaign data sets denoted above.These histograms show that the gradient depth has a large variation among the observational campaigns but a small variation during each campaign.As the average condition for the gradient depth, a normal distribution with a mean of 0.65 km and a standard deviation of 0.1 km is given as the default value for the prior distribution.These values were almost the mean values of the individual parameters, as shown in Additional file 1: Fig. S3.This prior distribution is useful for stabilizing solutions with insufficient moving survey data, although it also brings the solutions closer to those of the analysis method with a fixed gradient depth (e.g., Honsho et al. 2019) when the data have sufficient sensitivity to detect the gradient depth.Thus, the prior distribution parameters should be set carefully when sufficient moving survey data are available.Furthermore, the prior distribution parameters were obtained from relatively shallow GNSS-A sites (water depth of approximately 1-2 km) along the Japan Trench, and large outliers were sometimes detected.Thus, the validity of these values should be investigated in the future using deep-sea sites with sufficient moving survey data.

Data format
The GNSS-A positioning functions of SeaGap shown in Sect."GNSS-A positioning functions in SeaGap" were conducted using the following four contents of files: (1) offset between a GNSS antenna and sea surface transducer, (2) initial seafloor transponder positions, (3) sound speed profile, and (4) observation data.
Files (1), (2), and (3) contain information on b 0 , p 0 k , and v 0 in Eq. 2, respectively.File (4) contains travel times, observational period, GNSS antenna positions, and attitudes of the sea surface platform when the acoustic signals are transmitted and received at the sea surface platform.The initial seafloor transponder positions (File 2) and the GNSS antenna positions (File 4) were mapped into an ENU coordinate using the tangent-type transverse Mercator projection with its standard parallel intersecting the array center for horizontal components.The vertical component remains ellipsoidal height.The format of these data files is similar to that of the existing GNSS-A positioning software, GARPOS (Watanabe et al. 2020), but the above projection system is different from GARPOS.The information additionally required for Sea-Gap is the flag of the group number of acoustic shots in file (4) for kinematic array positioning (Sect."A kinematic array positioning function by NLLS").
In addition to the above files, a parameter setting file is required when MCMC-based array positioning methods are applied.This file lists the initial values, lower and upper limits, and step widths of the unknown parameters of the MCMC-based array positioning method.These values can be assigned by users.By setting the same values for the initial value and the lower and upper limits for a particular unknown parameter, the analysis can be performed with the unknown parameter fixed at the initial value.The parameter setting file can be easily created from the results of static array positioning using NLLS by a SeaGap function exclusively for this purpose.However, their values are plausible to be adjusted manually so that the acceptance ratio of the M-H algorithm (the ratio of the number that the candidate state x ′ is accepted to the total MCMC iterations) is approximately 23% (Gelman et al. 1996).

Applications
SeaGap functions were applied to actual observational data obtained by a charter ship "No. 3 Kaiyo-Maru" with a hull-mounted transducer at the G20 site (longitude: 142.0826°E, latitude: 36.1575°N,water depth: 2742.7 m) on November 13, 2015.As the observational data were obtained from the acoustic units of Tohoku University, travel times from all transponders were simultaneously collected in a single acoustic shot.We employed this data set, because it includes both the moving survey and fixed-point survey and has a clear underwater gradient structure as demonstrated by Tomita and Kido (2022).Evident outliers in the observational data were eliminated in advance.The track and time series of the sea surface platform positions are shown in Fig. 2a, b, respectively.

Kinematic array positioning
The kinematic array positioning function, named "kin-ematic_array" (Eq.14), was performed to each acoustic transmission.Figure 2c, d shows the estimated array displacements in the map view and time series, respectively.As the moving survey was conducted during the first 3 h, the estimated array displacements showed a larger dispersion than in later periods.This is because the horizontal positioning accuracy is degraded when the sea surface platform was apart from the array center (e.g., Imano et al. 2019) and a wide spatial-scale and temporally constant sound speed heterogeneity causes systematic biases depending on the sea surface platform position.The estimated array displacements show temporal variation also during the fixed-point survey, which is considered to be small-scale sound speed heterogeneity with short time periods, such as nonlinear cascades derived from internal waves (Matsui et al. 2019).In addition to the array displacements, the NTD ( C i in Eq. 14) was estimated as shown in the upper panel of Fig. 2d.Using the estimated array displacements during the fixed-point survey, the average array displacement of this campaign is obtained as 2.9 ± 3.7 cm (1σ) in the EW component and −1.6 ± 4.8 cm (1σ) in the NS component relative to the initial seafloor transponder positions.

Static array positioning by a non-linear least square method
The static array positioning method assuming an HS-SSS, named "static_array" (Eq.16), was performed.As explained in Sect."Static array positioning functions by NLLS", the temporal evolution of the NTD was modeled using cubic B-spline functions in the static array positioning methods of SeaGap.The total number of cubic B-spline bases expressed the NTD temporal smoothness.Thus, the AIC and BIC values were calculated to optimize the total number of cubic B-spline bases.The AIC and BIC are derived from different statistical assumptions, and which indicator should be employed depends on the data and the purpose.In general, the BIC is employed to select a correct model among the given models, while the AIC is employed to find the best model for predicting future observations (Chakrabarti and Ghosh 2011).We then adopted the BIC to search for the optimal total number of cubic B-spline bases by default.Practically, the AIC selected too large number of the cubic B-spline bases in terms of efficiently performing the inversion calculation for our observational data.Figure 3 shows the BIC values for various numbers of cubic B-spline bases, and the optimal total number of cubic B-spline bases was determined as 73 which corresponds to ~12 min interval).Although Fig. 3 shows multiple local minimums in BIC, these were caused depending on the consistency between the sampling interval of the travel time data and the interval of the cubic B-spline bases.The travel time data were obtained at the sampling interval of 30 s.If the interval of the cubic B-spline bases is a multiple of the sampling interval, the cubic B-spline curve well fits with the data and lowers the BIC value.Therefore, this tendency of the BIC values appeared in Fig. 3.
Each cross symbol in Fig. 4a indicates the projected travel time residuals T cal n,k in Eq. 16) in the nadir direction, and the black curve in Fig. 4a shows the NTD fluctuation modeled by the cubic B-spline functions.Figure 4b shows the projected travel time residual obtained by subtracting the modeled NTD fluctuations.As shown in Fig. 4a, the cubic B-spline functions expressed the temporal and spatial variations of the NTD during the moving survey (the first 3 h), and the final travel time residuals (the projected travel-time residuals that the NTD was subtracted) (Fig. 4b) did not show significant discrepancies between the moving and fixedpoint surveys.The RMS of the final travel-time residuals was ∼ 3.5 × 10 −5 sec.The estimated array displacement of this campaign relative to the initial seafloor transponder positions is obtained as 2.6 ± 0.2 cm (1σ), −1.6 ± 0.2 cm (1σ), and −3.4 ± 0.6 cm (1σ) in the EW, NS, and UD components, respectively.These standard deviations are obtained from the diagonal components of the covariance matrix.They represent small values within 1 cm, although the repeated positioning precision for each campaign typically ranges from a few to ~ 10 cm (e.g., Tomita et al. 2017;Yokota and Ishikawa 2020).This inconsistency appears, because the observation equation does not sufficiently consider systematic error sources such as the heterogeneity of an actual SSS, and the variance of the model parameters as a result of the inversion analysis is not necessarily an appropriate index to express actual measurement errors.However, the standard deviations from the covariance matrix are useful for assessing the resolution of the array displacement, depending on the geometric configuration of the sea surface platform track and seafloor transponder positions (e.g., Imano et al. 2019).Particularly, the sensitivity of the vertical array displacement strongly depends on the geometric configuration.Thus, the relative difference in the uncertainty of the vertical component among multiple campaigns provides valuable information for performing a time-series function fitting of the vertical array displacements among the campaigns.

Static array positioning method using an MCMC technique
The static array positioning method assuming an HS-SSS with temporally constant SLGS through the MCMC technique (Eqs.20 and 21) with the additional prior distributions ("static_array_mcmcgradc") was performed.The prior distribution for the gradient depth was a normal distribution with a mean of 0.65 km and a standard deviation of 0.15 km, and that for shallow gradients was a normal distribution with the mean of zero and the standard deviation of 0.1 ms/km.The total number of MCMC iterations and burn-in period were assigned as 1.2 × 10 6 and 2 × 10 5 , respectively, and the sampling interval was set to be 5. Figure 5 shows the variations in the log-likelihoods, log-prior probability densities, and the RMSs of the projected travel time residuals (the NTD and the gradient components were subtracted) for Eq.20 during the MCMC iterations.Additional file 1: Fig. S4 shows the variations in the primary unknown parameters (three components of the array displacement, shallow gradients, and gradient depth) during the MCMC iterations.These figures confirm the convergence of the unknown parameters during the MCMC iterations.The plots in these figures were randomly extracted from all MCMC samples for quick drawing and visibility.Note that the log-likelihood is divided into two components which correspond to i of Eq. ( 23) in Fig. 5 to check the convergence process for each observation equation (Eqs.20 and 21).
The sampled histograms of the unknown parameters of L-NTD occasionally show multiple steep peaks that may indicate inefficient update of these parameter during the MCMC method because of the strong covariance among the L-NTD parameters.However, this did not significantly affect the L-NTD modeling, as discussed in the next paragraph.Although it is plausible to introduce another sampling technique, such as Hamiltonian Monte Carlo, to cope with this issue, the current version of Sea-Gap does not implement it.
Figure 6 shows the fitness of the optimal models expressed by Eqs.20 and 21 and the spatial variation in the sound speed heterogeneity.To produce Fig. 6, An optimal value for each unknown parameter was calculated as mean of the MCMC samples for the corresponding parameter.The top panel of Fig. 6a shows the fitting of the model expressed by Eq. 20: The colored crosses show the projected travel-time residuals dem- onstrate contributions of the shallow gradients; their spatial pattern is shown in Fig. 6b.The positive and negative travel time residuals in Fig. 6b indicate lower and higher sound speeds, respectively.The bottom panel of Fig. 6a shows the final projected travel time residuals (the S-NTD and the equivalent deep gradients were subtracted) in Eq. 20.
Figure 7 shows the sampling results for the primary unknown parameters.As shown by the diagonal components in this figure, the parameters appear to follow a normal distribution.The horizontal array displacements and gradient depth were clearly correlated with each other, whereas shallow gradients had smaller correlations with the other parameters.The vertical array displacement did not show a correlation with the other parameters in Fig. 7; however, it did show a correlation with the parameters expressing the NTDs as shown in Additional corresponding to the solution obtained by NLLS.The solutions of the MCMC method in the horizontal components (Fig. 8a, b) are systematically different from those of the NLLS and exhibit larger variances than those of the NLLS.These differences are caused by the assumption of a sound-speed gradient in the MCMC method, and such differences due to the methods were discussed in detail in Tomita & Kido (2022).Contrastingly, the vertical array displacement of the MCMC method is comparable to that of the NLLS method, which implies that the assumption of the sound speed gradient does not significantly affect the vertical array displacement.The mean and standard deviation of the array displacements obtained as the MCMC samples are 6.1 ± 0.4 cm (1σ), −5.6 ± 0.5 cm (1σ), and −4.5 ± 0.5 cm (1σ) in the EW, NS, and UD components, respectively.The uncertainties in the horizontal components are slightly larger than those obtained by the NLLS in the HS-SSS by considering the sound speed gradient; however, they are still too small to express the actual measurement discrepancy among the campaigns for the same reasons explained in Sect."Static array positioning by a non-linear least square method".
Finally, the influence of the additional prior distributions for shallow gradients and gradient depth on the observational data without moving survey data was investigated.As discussed in Tomita and Kido (2022), the observational data without moving survey data had insufficient sensitivity to constrain shallow gradients and gradient depths, and these parameters were wrongly estimated.The observational data of G17 site (longitude: 142.7154°E, latitude: 36.8990°N,water depth: 4232.4 m) on October 16, 2017 were obtained only from the fixedpoint survey.The track of the sea surface platform in the observational data is shown in Additional file 1: Fig. S6.
Figure 9a shows the MCMC samples of the primary unknown parameters when the MCMC method was performed without the additional prior distributions ("static_array_mcmcgrad"). The shallow gradients had multiple peaks, and their intensity (~1-2 ms/km) were unnaturally high compared with the previous studies (< ~0.1 ms/km, Additional file 1: Fig. S2).Therefore, we considered that the sampled shallow gradients were not well-determined, because the observational data did not have enough sensitivity to solve them; consequently, the gradient depth and the horizontal array displacements were not also well-determined.In such a case that the survey data do not have sufficient sensitivity to estimate the shallow gradients, the positioning results assuming an HS-SSS would show more reasonable solutions as indicated by Tomita and Kido (2022).In terms of applying a uniform positioning method to all campaign observation data sets for obtaining comparable solutions among the data sets, it is not appropriate to vary the assumption of the positioning method (an HS-SSS or an HS-SSS with horizontal gradients) from one data set to another.We considered it was plausible to constraint the shallow gradients to be zero when the observational data are not sensitive to the shallow gradients, which would provide the solutions close to the those estimated by the positioning method assuming an HS-SSS.Moreover, as the shallow gradients were typically estimated to be close to zero in the previous studies (Additional file 1: Fig. S2), this constraint was considered to be plausible.
Figure 9b shows the solutions obtained by the MCMC method with the additional prior distributions of ("static_ array_mcmcgradc").The prior distribution for the gradient depth was a normal distribution with a mean of 0.65 km and a standard deviation of 0.15 km.The prior distribution for the shallow gradients was a normal distribution with the mean of zero and the standard deviation of 5 × 10 -2 ms/km.To avoid the unnatural deviation of the shallow gradients as shown in Fig. 9a, we provided the standard deviation much smaller than the statistic value of 0.15 ms/km.Note that we confirmed that this value did not strongly affect the estimation when applying it to the observational data including the moving survey data (Additional file 1: Fig. S7).Additional file 1: Fig. S7 shows the histograms of the MCMC samples using the same data set and analysis conditions with Fig. 7, except the standard deviation of prior distribution for the shallow gradients was 5 × 10 -2 ms/km.The mean values of the sampled horizontal array displacements were comparable to the solutions obtained by the NLLS assuming an HS-SSS (black curves), while the standard deviations of the sampled horizontal array displacements were much larger than those of the solutions assuming the HS-SSS (Fig. 9b).The larger standard deviations were attribute to the introduction of the horizontal gradients.The shallow gradients and the gradient depths were obtained following the prior distributions as the observational data did not have enough potential to constrain them.The histogram of the vertical array displacement shows multiple peaks, but this might be derived from inefficient sampling due to the trade-off relationship between the vertical array displacement and the NTD components.
When the observational data do not include the sufficient moving survey data, the additional prior distributions enable us to obtain the robust solutions, which are close to those assuming an HS-SSS, even by the framework considering horizontal gradients in the SSS.Although the horizontal gradient structure was not properly estimated in this case, the robustness for the positioning results due to the additional prior distributions is effective in universally applying a fixed analysis method to various observational data sets.In terms of this universal application to various data sets, it is important to explore appropriate values of the additional prior distributions in the future.

Other features of SeaGap
In addition to the aforementioned positioning functions, SeaGap has various other functions.One of these functions is outlier removal.Additional file 1: Fig. S6 shows an example of outlier removal from the observational data at G20 on November 13, 2015, as used in Sect."Applications".This function first calculated exact travel times, and then we obtained travel time residuals (blue and red crosses in the left panels).Next, we calculated a smoothed time series of the travel time residuals using a running average filter or running median filter.In Additional file 1: Fig. S7, a running median filter with a time window of 7 samples is employed to smoothen the traveltime residuals (magenta curve in the left panels).If outliers exist in the observed travel times, the original travel time residuals are apart from the smoothed travel time residuals curve.Subtracting the smoothed travel time residuals from the original travel time residuals, we calculated the standard deviation of the residuals (blue and red crosses in the right panels).The function removed the travel time data with their residuals beyond the arbitrary multiples of the standard deviation as outliers (red crosses in the left and right panels).In Additional file 1: Fig. S6, the residuals beyond the 4σ standard deviation were identified as outliers.The filter properties and criteria for identifying the outliers can be set arbitrarily.
SeaGap prepares a line-fitting function to simply calculate array displacement rates and their uncertainties.Figure 10 shows an example of the time series of the array displacements for multiple campaigns at G20 and the fitted lines with 95% credible intervals.The error bars in Fig. 10 demonstrated 10σ standard errors defined by the model covariance matrix.As indicated in Sect."Applications", the uncertainties obtained using static array positioning methods are generally small.Therefore, the 10σ standard deviation is empirically adopted as the error range (e.g., Honsho et al. 2019).Furthermore, SeaGap has various visualization functions; Figs. 2,3,4,5,6,7,8,9, 10 and Additional file 1: Figs.S4-S7 were generated using SeaGap's drawing functions.

Summary
We developed a new software for GNSS-A positioning "SeaGap" and applied it to actual observational data.The software has various GNSS-A positioning functions for both kinematic and static positioning, which cover the traditional kinematic array positioning method for individual acoustic shot groups (e.g., Kido et al. 2006) and the latest MCMC-based static array positioning method that considers a sound speed gradient (Tomita and Kido 2022).The format of the input data files was designed by referring to existing GNSS-A positioning software (Watanabe et al. 2020).All positioning methods in the software adopt the NTD to model the temporal fluctuation of the average sound speed, which leads to sophisticated observation equations.The MCMC-based static array positioning method (Tomita and Kido 2022) demonstrated unstable solutions when no moving survey data were available.We implemented a countermeasure for this issue into the positioning method of Tomita and Kido (2022), which constrained the unknown parameters expressing the sound speed gradient by providing additional prior distributions.If the observational data do not have sufficient sensitivity to solve the sound speed gradient, this new function provided robust solutions close to those assuming an HS-SSS.In addition to these positioning functions, SeaGap has various useful functions such as outlier removal, linear line-fitting, and drawing.
The existing positioning software, GARPOS, provides a single computational function, which can flexibly assign unknown parameters and enable us to perform various positioning analyses.The system of SeaGap is divided into multiple computational functions depending on the analysis conditions.GARPOS performs under an HS-SSS or an HS-SSS with temporally variational MLGS, while SeaGap performs under an HS-SSS or an HS-SSS with temporally constant SLGS.The former software is powerful to solve even a complex SSS, while SeaGap provides robust solutions using various observational data sets including a data set that contains fixed-point and moving surveys or a data set that only contain a fixed-point survey.
Regarding the limitations of SeaGap, the solutions may be unstable or biased when observational data deviating from the assumptions of the observation equations are employed.Because the SeaGap's positioning methods, which consider the sound-speed gradient, assume that the sound-speed gradient is constant over time, the temporal variation of the sound-speed gradient may provide inaccurate solutions.Thus, this issue should be considered when analyzing long-term observational data.For the employment of long-term observational data, the segmentation of observational data may be a good step to assess this issue.In addition, the SeaGap system assumes that the observational data are collected by single sea surface platform.Thus, it is a future task to handle observational data collected simultaneously by multiple sea surface platforms in SeaGap.The investigation of the prior distribution parameters that can be used universally for a variety of observational data is another issue.The developed software is available in an open-access repository, and is equipped with an online manual and tutorials.We hope that the increase in GNSS-A positioning tools and open access will promote the engagement of researchers in the seafloor geodesy field and help develop a methodology for GNSS-A positioning techniques in the future.

Fig. 1
Fig. 1 Schematic image of the GNSS-A observation.Incident angle ξ and azimuth φ of the acoustic path used in Methodology are shown.Red vectors show a common displacement among the seafloor transponders (i.e., array displacement)

Fig. 2 a
Fig. 2 a Track of a sea-surface platform of the observational data obtained at G20 on November 13, 2015 is shown by the cross symbols.The color represents the observational time.Gray triangles show the seafloor transponder positions.b Time series of the sea-surface platform positions shown in (a).c Estimated horizontal array displacements for individual acoustic shot groups using the observational data of a, b.The color represents the observational time.d Time series of the estimated horizontal array displacements and NTDs corresponding to the results of (c)

Fig. 3
Fig. 3 BIC values of various number of B-spline bases are shown in blue cross symbols, for the observational data obtained at G20 in Oct. 2015.Vertical red line shows the number of B-spline bases for the minimum BIC value

Fig. 5
Fig. 5 Convergence process of prior densities, likelihoods, summation of logarithms of prior densities and likelihoods, and RMSs during the MCMC iteration.a, b are logarithm of prior densities for the gradient depth and the shallow gradients, respectively.c, d are logarithm of likelihoods for Eqs.20 and 21, respectively.e is summation of (c) and (d).f is summation of (a), (b), and (e).g, h are RMSs of the projected travel time residuals for Eqs.20 and 21, respectively 1

Fig. 6 a
Fig. 6 a Travel-time residuals obtained by the MCMC-based static array positioning method are shown in time series.Top panel shows the projected travel-time residuals by the colored cross symbols and the black symbols are summation of the NTDs modeled as S-NTD and the equivalent deep gradients.Middle panel shows the projected travel-time residuals subtracting the equivalent deep gradients by the colored cross symbols.The black, red, and blue curves show S-NTD, L-NTD, and summation of L-NTD and the contributions of the shallow gradients.Bottom panel shows the final travel-time residuals subtracting NTDs.b Projected travel-time residuals subtracting L-NTD are shown on the sea-surface platform track

Fig. 7
Fig. 7 Histograms, heatmaps, and scatter maps for the primary unknown parameters sampled by the MCMC-based static array positioning method.The diagonal components show the histograms of the individual unknown parameters.The lower and upper triangle components show the heatmaps and the scatter maps among the unknown parameters, respectively

Fig. 8
Fig. 8 Cyan histograms in panels a-c show the posterior PDF for the EW, NS, and UD array displacements sampled by the MCMC-based static array positioning method, respectively.Each blue curve demonstrates a fitted normal distribution expressing the poster PDF.Each black curve demonstrates a normal distribution obtained by the NNLS-based static array positioning method assuming a horizontally stratified sound speed structure

Fig. 10
Fig. 10 Time series of the array displacements for multiple campaigns at G20 estimated by the MCMC-based static array positioning method.Estimation error for each campaign represents 10σ standard deviation.Each cyan line demonstrates a weighted linear fitting line.Each cyan area demonstrates 95% credible interval