Influence of Sea Level Anomaly on Underwater Gravity Gradient Measurements

Considering the theoretical research needs of gravity gradient detection and navigation, this study uses the right rectangular prism method to calculate the disturbing gravity gradient from sea level anomalies in the range of 5° × 5° in the Kuroshio extension area of the western Pacific with large sea level anomalies. The disturbing gravity gradient is obtained in different directions within a depth of 50 m below the mean sea level based on the principle of the disturbing gravity gradient. The calculation results show that the sea level anomalies at local positions significantly impact the underwater gravity gradient measurements, with the maximum contribution exceeding 10 E and the maximum difference between different locations exceeding 20 E. The change of the sea level anomaly over time also significantly impacts the measurement of the underwater gravity gradient, with the maximum change value exceeding 20 E. The impact will have a corresponding change with the seasonal change of the sea level anomaly. Therefore, the underwater carrier needs to consider the disturbing gravity gradient caused by sea level anomalies when using the gravity gradient for underwater detection and navigation.


Introduction
Gravity field is an inherent physical property of the earth. Due to differences in seabed topography and crustal density at different locations, the gravity field parameters at different locations in the ocean vary. Passive navigation can be realized using this feature. Gravity gradient denotes the rate of change of gravity in the entire space and is the second derivative of the gravity potential. Gradient anomalies can reflect the details of the field source and exhibit higher resolution than gravity. During detection and navigation using the gravity gradient, a carrier does not actively transmit signals to the outside, thus providing natural advantage for submarines with high concealment.
First, a high-precision and high-resolution gravity gradient reference map needs to be prepared and assembled into the navigation system. The gravity gradiometer is used to measure the gravity gradient data at the carrier location in real time. The measured data are compared with the data of the previously determined gravity gradient reference map. Subsequently, the location of the carrier can be determined by matching and solving according to a certain algorithm, and then, the output results of the inertial navigation can be corrected. The spatial resolution and accuracy of a gravity gradient reference map directly determine the accuracy of navigation and detection. Since the gravity gradiometer has not yet been employed in large-scale engineering applications, the gravity gradient background field needs to be calculated using other models. One possible method is using the Digital Elevation Model (DEM) to forward the gravity gradient [1][2][3][4]. When using DEM data to calculate the gravity gradient, the mean sea level or the geoid, which is approximated by the mean sea level, is used as the datum [5][6][7]. However, certain differences are present between the actual and mean sea levels.
Yan, Z. et al. conducted a quantitative analysis and calculation of the high-frequency wave of the sea surface, and they determined that for the sea state above level 4 (the average wave height is greater than 1.25 m), the fluctuation of the sea waves will cause significant disturbance in the measurement value of the gravity gradient on the underwater carrier [8]. In addition to the high-frequency waves, there are also low-frequency changes in sea level. Sea level denotes the average height of the sea surface over a period of time. It is the level that the ocean can maintain at a certain moment, assuming that no sea surface fluctuation is caused by tides, waves, or other factors. Based on the observation time, the sea level can be divided into daily mean sea level, monthly mean sea level, annual mean sea level, long-term mean sea level, etc., which are affected by weather, climate, or celestial movement and change. Moreover, long-term global trend changes, periodic changes, and large spatial differences can be observed at different locations [9][10][11]. This study aims to determine how the mass anomaly caused by the difference of the sea level height affects the gravity gradient, which could provide a theoretical reference for the practical application of gravity gradient detection and navigation in the future.

Data and Methods
This study uses daily gridded sea level data from 1993 to the present from global ocean satellite observations that were provided by the Copernicus Marine Environment Monitoring Service (CMEMS) and the Copernicus Climate Change Service (C3S). The dataset provides gridded daily global estimates of sea level anomalies based on satellite altimetry with a spatial resolution of 0.25 • × 0.25 • , including absolute dynamic terrain, sea level anomaly, absolute geostrophic velocity, and geostrophic velocity anomaly. Moreover, sea level anomalies were calculated relative to the 20-year average reference period (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012). The dataset is based on a satellite constellation with a stable number of altimeters ( Figure 1) to ensure long-term stability of the ocean observing system. After the steps of satellite acquisition, preprocessing, input data quality control, multi-task cross-calibration, along-orbit product generation, gridded product generation, and final quality control, estimates of sea level anomalies, absolute dynamic topography and corresponding geostrophic velocities were obtained [12][13][14]. DEM data to calculate the gravity gradient, the mean sea level or the geoid, which is approximated by the mean sea level, is used as the datum [5][6][7]. However, certain differences are present between the actual and mean sea levels. Yan, Z. et al. conducted a quantitative analysis and calculation of the high-frequency wave of the sea surface, and they determined that for the sea state above level 4 (the average wave height is greater than 1.25 m), the fluctuation of the sea waves will cause significant disturbance in the measurement value of the gravity gradient on the underwater carrier [8]. In addition to the high-frequency waves, there are also low-frequency changes in sea level. Sea level denotes the average height of the sea surface over a period of time. It is the level that the ocean can maintain at a certain moment, assuming that no sea surface fluctuation is caused by tides, waves, or other factors. Based on the observation time, the sea level can be divided into daily mean sea level, monthly mean sea level, annual mean sea level, long-term mean sea level, etc., which are affected by weather, climate, or celestial movement and change. Moreover, long-term global trend changes, periodic changes, and large spatial differences can be observed at different locations [9][10][11]. This study aims to determine how the mass anomaly caused by the difference of the sea level height affects the gravity gradient, which could provide a theoretical reference for the practical application of gravity gradient detection and navigation in the future.

Data and Methods
This study uses daily gridded sea level data from 1993 to the present from global ocean satellite observations that were provided by the Copernicus Marine Environment Monitoring Service (CMEMS) and the Copernicus Climate Change Service (C3S). The dataset provides gridded daily global estimates of sea level anomalies based on satellite altimetry with a spatial resolution of 0.25° × 0.25°, including absolute dynamic terrain, sea level anomaly, absolute geostrophic velocity, and geostrophic velocity anomaly. Moreover, sea level anomalies were calculated relative to the 20-year average reference period (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012). The dataset is based on a satellite constellation with a stable number of altimeters ( Figure 1) to ensure long-term stability of the ocean observing system. After the steps of satellite acquisition, preprocessing, input data quality control, multi-task crosscalibration, along-orbit product generation, gridded product generation, and final quality control, estimates of sea level anomalies, absolute dynamic topography and corresponding geostrophic velocities were obtained [12][13][14]. Gravity gradient forward modeling is based on Newton's law of gravity and utilizes known mass anomalies to calculate the gravity gradient. The methods for the gravity gradient forward modeling from the terrain mainly include right rectangular prism, polyhedron, direct numerical integration, and fast Fourier transform method. Comparison of the Gravity gradient forward modeling is based on Newton's law of gravity and utilizes known mass anomalies to calculate the gravity gradient. The methods for the gravity gradient forward modeling from the terrain mainly include right rectangular prism, polyhedron, direct numerical integration, and fast Fourier transform method. Comparison of the various methods showed that the right rectangular prism method affords results with the highest accuracy [2,3,16]; thus, this study uses the right rectangular prism method for calculations.

Spatial and Temporal Distribution of Sea Level Anomalies
Since the sea level has a significant one-year cycle and the carrier generally does not operate underwater for more than one year, this study chooses 2020 as the representative year. The daily mean sea level anomaly data of 366 days from 1 January 2020 to 31 December 2020 are used to calculate and analyze the impact of sea level height anomaly on the underwater gravity gradient measurements. The spatial distribution of the annual mean sea level anomaly in 2020 ( Figure 2) was obtained by time-averaging the 366-day data. The results show that most of the world's oceans were above the long-term mean sea level in 2020, which is consistent with the current trend of rising sea levels due to global warming. Significant regional differences in the sea level are also observed.  various methods showed that the right rectangular prism method affords results with the highest accuracy [2,3,16]; thus, this study uses the right rectangular prism method for calculations.

Spatial and Temporal Distribution of Sea Level Anomalies
Since the sea level has a significant one-year cycle and the carrier generally does not operate underwater for more than one year, this study chooses 2020 as the representative year. The daily mean sea level anomaly data of 366 days from 1 January 2020 to 31 December 2020 are used to calculate and analyze the impact of sea level height anomaly on the underwater gravity gradient measurements. The spatial distribution of the annual mean sea level anomaly in 2020 ( Figure 2) was obtained by time-averaging the 366-day data. The results show that most of the world's oceans were above the long-term mean sea level in 2020, which is consistent with the current trend of rising sea levels due to global warming. Significant regional differences in the sea level are also observed.  The difference between the maximum and minimum values in the 366 days at each grid point is used to obtain the annual range of global sea level anomalies in 2020 ( Figure  3), representing the variation range of the sea level height in 2020. Figure 3 shows that the high-value areas of the annual range are mainly concentrated in the western boundary flow area of the ocean and the Southern Ocean, with amplitudes exceeding 1 m, which indicates that the sea level significantly changes in these areas during the year. The sea area near the southern Cape of Good Hope in Africa (39.4° S, 19.6° E) has an amplitude of 2.11 m, which exhibits the maximum annual range of global sea level anomalies in 2020. The difference between the maximum and minimum values in the 366 days at each grid point is used to obtain the annual range of global sea level anomalies in 2020 ( Figure 3), representing the variation range of the sea level height in 2020. Figure 3 shows that the high-value areas of the annual range are mainly concentrated in the western boundary flow area of the ocean and the Southern Ocean, with amplitudes exceeding 1 m, which indicates that the sea level significantly changes in these areas during the year. The sea area near the southern Cape of Good Hope in Africa (39.4 • S, 19.6 • E) has an amplitude of 2.11 m, which exhibits the maximum annual range of global sea level anomalies in 2020.
The above analysis shows that the sea level anomaly varies across different regions of the world with a maximum difference of more than 1 m, and the amplitude of changes within a year also differs, with a maximum change of 2 m. Therefore, this study aims to determine whether the existence of sea level anomalies and their variation with the spatiotemporal distribution affect the underwater gravity gradient measurements.  The above analysis shows that the sea level anomaly varies across different regions of the world with a maximum difference of more than 1 m, and the amplitude of changes within a year also differs, with a maximum change of 2 m. Therefore, this study aims to determine whether the existence of sea level anomalies and their variation with the spatiotemporal distribution affect the underwater gravity gradient measurements.

Simulation Calculation of Gravity Gradient
Considering the large amplitudes of the annual mean and annual range in the Kuroshio extension area in Figures 2 and 3, the right rectangular prism method is used to calculate the disturbing gravity gradient from sea level anomalies in the range of 5° × 5° (34° N-39° N, 152° E-157° E) in the Kuroshio extension area with large sea level anomalies. Moreover, the data of the day with the highest sea level anomaly in the region (8 October 2020) is selected from the 366-day data for subsequent calculation ( Figure 4 and Table 1).

Simulation Calculation of Gravity Gradient
Considering the large amplitudes of the annual mean and annual range in the Kuroshio extension area in Figures 2 and 3, the right rectangular prism method is used to calculate the disturbing gravity gradient from sea level anomalies in the range of 5 in the Kuroshio extension area with large sea level anomalies. Moreover, the data of the day with the highest sea level anomaly in the region (8 October 2020) is selected from the 366-day data for subsequent calculation ( Figure 4 and Table 1). The above analysis shows that the sea level anomaly varies across different regions of the world with a maximum difference of more than 1 m, and the amplitude of changes within a year also differs, with a maximum change of 2 m. Therefore, this study aims to determine whether the existence of sea level anomalies and their variation with the spatiotemporal distribution affect the underwater gravity gradient measurements.

Simulation Calculation of Gravity Gradient
Considering the large amplitudes of the annual mean and annual range in the Kuroshio extension area in Figures 2 and 3, the right rectangular prism method is used to calculate the disturbing gravity gradient from sea level anomalies in the range of 5° × 5° (34° N-39° N, 152° E-157° E) in the Kuroshio extension area with large sea level anomalies. Moreover, the data of the day with the highest sea level anomaly in the region (8 October 2020) is selected from the 366-day data for subsequent calculation ( Figure 4 and Table 1).   Due to the small longitudinal and latitudinal span of the study area, the calculation area is regarded as a plane, and a local Cartesian coordinate system is established on the mean sea level with the x, y, and z axes pointing toward the east, north, and upward  Figure 5). The projection of the computation point on the mean sea level is set as the origin, which will change the horizontal coordinates of the computation point P to 0, subsequently simplifying the calculation. Ignoring the influence of the centrifugal force generated by the earth's rotation, the sea level anomaly data are divided into rectangular prisms based on the longitude and latitude grid, whose length and width are determined as the average length and width of all grids and the upper and lower surfaces of the prism are determined as the sea level anomaly surface and mean sea level, respectively. Then, the horizontal position of the computation point is at the center of each grid ( Figure 6).

Lat. and Lon.
Resolution Grids Max (m) Min (m) STD (m) Mean (m) Due to the small longitudinal and latitudinal span of the study area, the calculation area is regarded as a plane, and a local Cartesian coordinate system is established on the mean sea level with the x, y, and z axes pointing toward the east, north, and upward directions ( Figure 5). The projection of the computation point on the mean sea level is set as the origin, which will change the horizontal coordinates of the computation point P to 0, subsequently simplifying the calculation. Ignoring the influence of the centrifugal force generated by the earth's rotation, the sea level anomaly data are divided into rectangular prisms based on the longitude and latitude grid, whose length and width are determined as the average length and width of all grids and the upper and lower surfaces of the prism are determined as the sea level anomaly surface and mean sea level, respectively. Then, the horizontal position of the computation point is at the center of each grid ( Figure 6).  To simplify the formula, using , , and to represent x, y, and z, respectively, we obtain the gravity gradient generated by the ith prism at the computation point P: Here, denotes Newton's gravitation constant, which is taken as 6.67 × 10 −11 m 3 ·kg −1 ·s −2 ; is the seawater density of the ith prism, which is assumed to be a constant value of 1025 kg·m −3 ; ℎ is the sea level anomaly value of the ith prism; and (a , b ) is the horizontal coordinate of the center point of the ith prism. Additionally, Δ and Δ are the average length and width, respectively, of each grid in the study area, which are cal-  Due to the small longitudinal and latitudinal span of the study area, the calculation area is regarded as a plane, and a local Cartesian coordinate system is established on the mean sea level with the x, y, and z axes pointing toward the east, north, and upward directions ( Figure 5). The projection of the computation point on the mean sea level is set as the origin, which will change the horizontal coordinates of the computation point P to 0, subsequently simplifying the calculation. Ignoring the influence of the centrifugal force generated by the earth's rotation, the sea level anomaly data are divided into rectangular prisms based on the longitude and latitude grid, whose length and width are determined as the average length and width of all grids and the upper and lower surfaces of the prism are determined as the sea level anomaly surface and mean sea level, respectively. Then, the horizontal position of the computation point is at the center of each grid ( Figure 6).  To simplify the formula, using , , and to represent x, y, and z, respectively, we obtain the gravity gradient generated by the ith prism at the computation point P: Here, denotes Newton's gravitation constant, which is taken as 6.67 × 10 −11 m 3 ·kg −1 ·s −2 ; is the seawater density of the ith prism, which is assumed to be a constant value of 1025 kg·m −3 ; ℎ is the sea level anomaly value of the ith prism; and (a , b ) is the horizontal coordinate of the center point of the ith prism. Additionally, Δ and Δ are the average length and width, respectively, of each grid in the study area, which are cal-   To simplify the formula, using x 1 , x 2 , and x 3 to represent x, y, and z, respectively, we obtain the gravity gradient generated by the ith prism at the computation point P: Here, G denotes Newton's gravitation constant, which is taken as 6.67 × 10 −11 m 3 ·kg −1 ·s −2 ; ρ i is the seawater density of the ith prism, which is assumed to be a constant value of 1025 kg·m −3 ; h i is the sea level anomaly value of the ith prism; and (a i , b i ) is the horizontal coordinate of the center point of the ith prism. Additionally, ∆x 1 and ∆x 2 are the average length and width, respectively, of each grid in the study area, which are calculated to be 22,371.490 and 27,742.137 m. Furthermore, r = is the distance between the computation point P (x 1 , x 2 , x 3 ) and element point P' (x 1 , x 2 , x 3 ), j and k are the subscripts of x, which can be taken as 1, 2, and 3, representing the x, y, and z directions, respectively. Substitute x 1 = 0, x 2 = 0, and x 3 = Z into Equation (1), place ρ before the integral sign, and replace x 1 , x 2 with u, v to represent the horizontal coordinates of the element point. Next, integrate Equation (1) in the vertical direction to obtain six components: Z denotes the depth from the mean sea level, which is a negative value. Selecting a certain depth and substituting it into the formula can yield the components of the gravity gradient. When the selected depth Z = 0, Equations (2)-(7) can be simplified: Each gravity gradient component of prism I at point P can be calculated using Equations (2)-(7), and its value at point P can be obtained by adding the disturbing gravity gradient values of all prisms in the calculation area at point P. Then, the disturbing gravity gradient caused by the sea level anomaly in the study area can be obtained through point by point calculation of the computation point of all grid centers [17].
The horizontal coordinate of the center point of the ith prism (a i , b i ) in Equation (1) should cover all mass elements since all mass elements will affect the gravity gradient at the computation point; however, this is difficult to achieve in practical calculations. Equation (1) shows that the gravity gradient rapidly decreases with increasing distance, denoting that the contribution of distant mass elements to the gradient value is very small. Therefore, the range of the calculation area needs to be determined to consider the balance between the calculation efficiency and error. Jekeli, C. developed a systematic algorithm method to determine the range of terrain data [18]. Wu, L. took the relative error of 1% as the truncation distance to calculate the gravity gradient effect [7]. We adopt Wu's method and select the locations corresponding to the highest value, the lowest value, and any other value in the study area as test points, and then, gradually increase the side length of the integration area to calculate the vertical disturbing gravity gradient at different depths below the mean sea level caused by the sea level anomalies (Figure 7 and Tables 2-4). Figure 7 shows that with increasing integration area, the vertical disturbing gravity gradient tends to be stable and the change decreases. When the gradient tends to stabilize, the side length of the integral region increases as the observation depth increases. To achieve a relative error of less than 1%, the side length of the three test points at 10 m below the mean sea level in the integration area needs to reach 24 m, which is 2.4 times the observation point depth. The side length of the three test points at 30 m below the mean sea level needs to reach 68 m, which is 2.3 times the observation point depth. The side length of the three test points at 50 m below the mean sea level needs to reach 108 m, which is 2.2 times the observation point depth. Furthermore, the gravity gradient forward modeling result is related to the integration range, and the selection of the integration range is related to the calculation height (depth). Therefore, to obtain relatively accurate forward modeling results, different integration ranges need to be considered when calculating at different heights (depths).
by point calculation of the computation point of all grid centers [17].
The horizontal coordinate of the center point of the ith prism (a , b ) in Equation (1) should cover all mass elements since all mass elements will affect the gravity gradient at the computation point; however, this is difficult to achieve in practical calculations. Equation (1) shows that the gravity gradient rapidly decreases with increasing distance, denoting that the contribution of distant mass elements to the gradient value is very small. Therefore, the range of the calculation area needs to be determined to consider the balance between the calculation efficiency and error. Jekeli, C. developed a systematic algorithm method to determine the range of terrain data [18]. Wu, L. took the relative error of 1% as the truncation distance to calculate the gravity gradient effect [7]. We adopt Wu's method and select the locations corresponding to the highest value, the lowest value, and any other value in the study area as test points, and then, gradually increase the side length of the integration area to calculate the vertical disturbing gravity gradient at different depths below the mean sea level caused by the sea level anomalies (Figure 7 and Tables 2-4). Figure 7 shows that with increasing integration area, the vertical disturbing gravity gradient tends to be stable and the change decreases. When the gradient tends to stabilize, the side length of the integral region increases as the observation depth increases. To achieve a relative error of less than 1%, the side length of the three test points at 10 m below the mean sea level in the integration area needs to reach 24 m, which is 2.4 times the observation point depth. The side length of the three test points at 30 m below the mean sea level needs to reach 68 m, which is 2.3 times the observation point depth. The side length of the three test points at 50 m below the mean sea level needs to reach 108 m, which is 2.2 times the observation point depth. Furthermore, the gravity gradient forward modeling result is related to the integration range, and the selection of the integration range is related to the calculation height (depth). Therefore, to obtain relatively accurate forward modeling results, different integration ranges need to be considered when calculating at different heights (depths).

Experimental Results and Analysis
When the relative error is less than 1%, 24, 68, and 108 m can be taken as the side lengths of the integration area to calculate the disturbing gravity gradient at 10, 30, and 50 m below the mean sea level, respectively. Since the average grid side length of the sea level anomaly data is about 20 km, which is significantly larger than the side length of the integration area, the gravity gradient value of each grid is only calculated from the mass anomaly in the central area of the grid, i.e., only one prism is calculated at each grid. Only the three components of Γ xx , Γ yy , and Γ zz are calculated and analyzed below since the Γ xy , Γ xz , and Γ yz of a single prism on the vertical line where the midpoint is located are zero [19]. Figure 8 shows that the sum of Γ xx , Γ yy , and Γ zz is close to zero, satisfying the Laplace equation, and thus, the calculation process can be considered to be correct. The chart shows that the sea level anomaly considerably impacts the disturbing gravity gradient value at 10 m below the mean sea level, up to 16 E, and the difference between the disturbing gravity gradient at different locations can reach 22 E. At 30 and 50 m below the mean sea level, the spatial distribution of Γ xx , Γ yy , and Γ zz is similar to the distribution at 10 m below the mean sea level, but the magnitude is lower; thus, only the statistics of the calculation results are listed. Table 5 shows that the disturbing gravity gradient caused by the sea level anomaly significantly decreases with increasing depth. The vertical disturbing gravity gradient can still reach 3 E at 50 m below the mean sea level, and the maximum difference of the disturbing gravity gradient at different locations can reach 4 E. Therefore, based on the measurement accuracy of the gravity gradiometer 1 E [20][21][22][23][24][25], the local sea level anomaly will considerably affect the underwater gravity gradient measurements at the depth of 50 m below the mean sea level.
The above results show that in addition to the spatial differences, the sea level anomalies will change with time. The annual range of the sea level anomalies in this area (Figure 9) is used to forward calculate the disturbing gravity gradient. The maximum and minimum values of the sea level anomaly are used as the upper and lower surfaces of the prism, respectively. The calculation process is the same as above. We obtain the change value of the disturbing gravity gradient caused by sea level anomaly within one year ( Figure 10). Figure 10 and Table 6 show that if the carrier measured the gravity gradient for one year in 2020 at 10 m below the mean sea level in the study area, the variation of the vertical disturbing gravity gradient caused by the sea level anomaly variation with time could reach a maximum of 24 E and minimum of 3 E, and the variation of the horizontal disturbing gravity gradient could reach a maximum of 10 E and minimum of 1 E. At 50 m below the mean sea level, the maximum variation of the vertical disturbing gravity gradient is more than 4 E, while the minimum variation is less than 1 E. We select the location of the maximum annual range in the study area (36.6 • N, 154.4 • E) and calculate the disturbing gravity gradient for one year in 2020 ( Figure 11). We find that the disturbing gravity gradient will be larger in summer and autumn compared to the other seasons, which is closely related to the seasonal changes in sea level anomalies. Therefore, within a depth of 50 m below the mean sea level, the change of the local sea level anomaly with time will have different effects on the underwater gravity gradient measurements, and this effect will change according to the seasonal changes of the sea level anomalies.  The above results show that in addition to the spatial differences, the sea level anomalies will change with time. The annual range of the sea level anomalies in this area ( Figure  9) is used to forward calculate the disturbing gravity gradient. The maximum and minimum values of the sea level anomaly are used as the upper and lower surfaces of the prism, respectively. The calculation process is the same as above. We obtain the change value of the disturbing gravity gradient caused by sea level anomaly within one year (Figure 10). Figure 10 and Table 6 show that if the carrier measured the gravity gradient for one year in 2020 at 10 m below the mean sea level in the study area, the variation of the vertical disturbing gravity gradient caused by the sea level anomaly variation with time could reach a maximum of 24 E and minimum of 3 E, and the variation of the horizontal disturbing gravity gradient could reach a maximum of 10 E and minimum of 1 E. At 50 m below the mean sea level, the maximum variation of the vertical disturbing gravity gradient is more than 4 E, while the minimum variation is less than 1 E. We select the location of the maximum annual range in the study area (36.6° N, 154.4° E) and calculate the disturbing gravity gradient for one year in 2020 ( Figure 11). We find that the disturbing gravity gradient will be larger in summer and autumn compared to the other seasons, which is closely related to the seasonal changes in sea level anomalies. Therefore, within a depth of 50 m below the mean sea level, the change of the local sea level anomaly with time will have different effects on the underwater gravity gradient measurements, and this effect will change according to the seasonal changes of the sea level anomalies.

Conclusions
Herein, we processed the daily grid data of sea level observed by global ocean satellites and quantitatively calculated the disturbing gravity gradient caused by sea level anomalies using the right rectangular prism method. The results show the following: (1) The gravity gradient forward modeling results are related to the integration range, and the selection of the integration range is related to the calculation height (depth). Therefore, to obtain relatively accurate forward modeling results, different integration ranges need to be considered when calculating at different heights (depths). (2) Based on the measurement accuracy of the gravity gradiometer 1 E, within 50 m below the mean sea level, sea level anomalies at local positions will significantly affect the underwater gravity gradient measurements, with a maximum contribution exceeding 10 E and the maximum difference between different locations exceeding 20 E. Moreover, the change of the sea level anomalies with time will significantly impact the underwater gravity gradient measurements, with the maximum change value exceeding 20 E, and the impact will accordingly change with the seasonal change of the sea level anomalies. Therefore, underwater carriers need to consider the disturbing gravity gradient caused by sea level anomalies when using gravity gradient for underwater detection and navigation.