On the performance of equiangular mascon solution in GRACE-like missions

Mass concentration (mascon) solutions for GRACE (Gravity Recovery and Climate Experiment) data are widely used in various regional-to-global mass change studies. The current advances in the mascon solution have mainly concentrated on improving the spatial resolution of the solution, enhancing the applied least-squares regularization, and the characterization of the solution errors. Most of the mascon solutions are obtained on the equal-area grid, inducing complexities in creating the grid and its presentation. In this regard, estimation of the mascon solutions on equiangular grids can be appealing. Furthermore, in the equal-area methods, there is no global criterion to determine the size of the mascon areas. The mascon size is usually chosen in a subjective manner which hampers the objective application of different mascon solutions. In view of these challenges, two main questions are addressed in this study: i) what kind of modifications should be made in computation scheme of the mascon solution if equiangular grids are used to account for different areas of the grid patches, and ii) in case of non-equiangular solutions, how to define an objective criterion for the patch sizes based on the resolution of both the observation and the signal of interest. We investigate the performance of the high-resolution mascon-based approach, proposed by Abedini et al. [2021], which uses GRACE-like observations similar to level-1 data for a period of one month over the Greenland region. Two main practical issues are studied on the estimation of the surface density changes as follows. First, we show that for equiangular grids, the area of the patches should be accounted for in the regularization by introducing area-affected weights for the unknown parameters. We investigate the effect of three different area-affected weighting strategies on the derived solution. Secondly in order to obtain proper size for the patches, a novel approach is presented to investigate the performance of the mascon solution using the analysis of the resolution matrix entries. The proposed resolution analysis is used to obtain the optimal patch size for the discretization of the area of interest. Based on the results, it is demonstrated that the minimum legible patch size in the Greenland area for the current settings of the GRACE observations is 0.5 degree in the NS direction and a latitude-adaptive grid-size rather than equiangular grids at high latitude regions in the EW direction.


Introduction
In recent years, many studies have been performed in the hydrological applications using data extracted from satellite's missions like Gravity Recovery and Climate Experiment (GRACE) [Awange et al., 2011;Bettadpur, 2007;Loomis et al., 2019;Luthcke et al., 2013;Sabaka et al., 2010;Save et al., 2012;Scanlon et al., 2016;Watkins et al., 2015;Wiese et al., 2016].The main mission of the GRACE satellites was to monitor the time-variable gravity field, which indicates mass changes due to water movement in the Earth's water cycle [Scanlon et al., 2016].The mass changes could be due to the changes in the (under) ground water storage, ice-sheet melting and ocean mass distribution [Luthcke et al., 2013;Rodell et al., 2009;Velicogna et al., 2014].The Earth's surface density can generally be estimated by analysis of GRACE data with two different methods.
The first group uses spherical harmonics (SH) to present the mascon solutions [Bettadpur, 2007;Wahr et al., 1998], and the second group of methodologies uses the Stokes coefficients to estimate the mass concentration from the satellite-to-satellite-tracking (SST) data [Luthcke et al., 2013;Save et al., 2016;Watkins et al., 2015].A comprehensive review of different methods and their pros and cons have been provided by Save et al. [2016], Scanlon et al. [2016], Jing et al. [2019] and Abedini et al. [2021].
In our previous study [Abedini et al., 2021], we suggested a mascon approach that uses the piecewise constant surface density functions, as local basis functions, to estimate the surface density.While other traditional mascon methods (e.g., in Luthcke et al. [2013], Rowlands et al. [2010], Save et al. [2016] and Watkins et al. [2015] have applied global basis functions, our main motivation was to avoid the series truncation that is usually applied in these mascon methods.The truncations of high frequency terms result in aliasing and leakage errors affecting the signals of interest at lower frequencies.The advantage of our proposed mascon method is that it avoids the truncation and henceforth aim at improving the leakage problem at the cost of extensive computation time [Abedini et al., 2021]. Disregarding what methodology is used, most of mascon solutions are obtained on approximately equal-area grids.For example, the JPL mascon solution has been obtained on series of equal-area 3° spherical cap mascon to cover the surface of the Earth, Goddard Space Flight Center (GSFC) group applied 2° × 2° arc-degrees equal-area mascon regions (~49 000  2 ), and the CSR group used a geodesic equal-area grid with collection of hexagonal cells, as size of each cell is chosen to be equivalent to an equiangular patch at equator (see Save et al. [2016], Scanlon et al. [2016] and Jing et al. [2019] for more information on these methods).
Nevertheless, there are two fundamental problems in equal-area methodologies.The first problem is related to the implementation complexity compared to simple equiangular latitude-longitude-based grids [Save et al., 2016].For example, for the spectral analysis of the mascon solutions (i.e., when Fast Fourier transform (FFT) is needed to be applied to a solution), an equiangular sampling is required [Klees et al., 2008].Furthermore, some of the equal-area mascon solutions [such as Save et al., 2016] do not cover the entire area of interest, resulting in gaps between spherical cap mascon.Also creating the equiangular grids and its presentation are much easier compared to the equal-area methods.Secondly, in the equal-area methods, the area/size of patches/cells are usually selected arbitrarily, and there is no objective criterion capable of accounting for both the resolution of the observation system and the required resolution of the solution.This limitation induces some degree of subjectivity in the selection of the patch sizes in different methods.
Following our previous study, this paper is rather an experimental investigation on the two issues explained above.Here, we first show that, in case of equiangular grids, the size of the mascon should be accounted for in the applied regularization by introducing area-affected weights for the unknown parameters.This is motivated by our previous study [Abedini et al., 2021], where we showed that the patch area is directly related to the variance of the estimated surface densities: the smaller the area is, the higher the estimated variance will be.Therefore, in the estimation procedure, we suggest introducing weights for the patches proportional to their area/size.Specially, when using equiangular patch sizes, in the high latitude regions, the area-affected regularization can be significant.
We investigate four different area-affected weighting strategies on the derived solution.Secondly, given the characteristics of the GRACE mission, we investigate a challenging problem on determining the size of the patches of which the finest resolution can be achieved for the mascon solution over the Greenland region.To obtain the optimal mascon size for non-equiangular grids, a novel approach is presented based on the resolution and spectral analysis of the obtained mascon solution.In this approach, we analyze the patterns of the rows of the resolution matrix to define the legible range for the mascon solutions.

Abbas Abedini et al.
The rest of the paper is organized as follows.In Section 2, a brief review of the mascon methodology is given followed by the explanation of the simulation setting.Sections 3 elaborates on the regularized solution and the proposed resolution analysis using the resolution matrix and its approximation by a sinc function.The results and discussions are given in Section 4. Section 5 provides a summary and makes some conclusions.

Theoretical background and further development
This section consists of four subsections.The first two subsections briefly review the mascon method proposed by Abedini et al. [2021].The last two sections describe further development on the regularization of the mascon solution along with its precision description.

Data set description
The area of interest (AOI), in this research, includes the Greenland and its surrounding ocean waters, which are considered to be a region within latitude (58°, 86°) and longitude (-72°, -14°) and a 16-degree extended area around the region to analyze its effect on the estimated surface density changes (SDCs) when the number of observations increases.The data set used is the simulated level-1 GRACE data for a period of one month.The observations are range rates between the two GRACE satellites at different time instants on the passing arcs over/around the Greenland region.That is, the arcs of the satellites' path are extracted from those that are crossing over the AOI to analyze its effect on the estimated surface density changes (SDCs).The mean SDC over the region is assumed to be  ≈ 250 / 2 in the coastal areas [Harig and Simons, 2012].In the inland, the mean SDC is less than this value, and outside Greenland it is assumed to be zero.Here, the unknowns are the surface density changes, which are estimated on equiangular surface patches.Two different spatial resolutions (2°×2° and 1°×1° patch-sizes) will be investigated to perform the analyses.Different analyses will be performed on these data sets.

Review of mascon method
This study is based on the mascon approach presented by Abedini et al. [2021].The method avoids the series truncation, usually used in other mascon methods, and hence aims at mitigating the usual leakage problem at the cost of extensive numerical burden.In the following, the concept of this method is briefly explained.
Let us assume that the vector  contains  observations and it is written as the summation of the actual range rates  of the two GRACE satellites and an additional noise component , i.e.  =  +.By taking the partial derivatives of the range rates with respect to the surface densities , the design matrix is obtained as  = , which introduces a vector of  unknowns of surface density changes  = .The overdetermined system of observation equations, with  number of observations and  number of unknowns ( > ), is written as where  is the  ×  covariance matrix of  and  is the dispersion operator.The element of the th row and th column of the design matrix  is defined as: [] =  /.The least squares estimate of the unknown  parameters is then obtained as where  = ⁻ ¹ is the weight matrix of the observation.There are some complications about the above least squares problem, which will be elaborated in the present contribution.Among them, the linear system in Eq. ( 1) is usually ill-posed.In other words, the design matrix  is nearly rank deficient and therefore the normal matrix  =  cannot be regularly inverted.To find a solution for the ill-posed problem, two regularization methods, namely the ordinary Tikhonov method and the so-called area-affected approach were suggested by Abedini et al. [2021].

𝜕𝜌 𝜕𝜎
The Tikhonov regularized solution is provided by [Hansen, 1996;Tihonov, 1963]: where  is the unknown regularization parameter.Eq. (3) provides a regularized solution of Eq. ( 1) for which an appropriate value for  is to be determined.In this regularization approach, all patches have imposed the same weight in the regularization via the  ×  identity matrix  in the right-hand side of Eq. ( 3).Alternatively, it is possible to give different weights to patches based on their spatial coverage (or area).To tackle this problem, it is perhaps easier to describe the regularization as a least square solution of an extended model with pseudo observations, which will be explained in details in the next subsection.

Prior information and regularization methods
One way to obtain a regular solution for the ill-posed problem is to take possible prior information on the unknown parameters into account [Koch and Kusche, 2002].The regularization can be performed by introducing the pseudo observations into the linear model in Eq. ( 1).This provides an extended linear model of observation equations as (4) with the joint covariance matrix: (5) where  ( × ) and  ( × ) are the cofactor matrices of the original and pseudo observations, respectively, and ₀² is the variance component of the original observations whereas ² is the regularization variance of the pseudo observations.The parameters ₀² and ² are a-priori unknown and they should be estimated/determined.The covariance matrix in Eq. ( 5) can be inverted to (6) where  = ₀² ⁻¹ is the weight matrix of the observations  and  = ⁻¹ is the weight matrix of ₀, the vector of pseudo observations;  = ⁻² is the regularization parameter.The regularized least squares solution in Eq. (3) can then be modified to Further, the least squares estimate of the residuals are obtained as (8) where  =  −  is the residuals of the original observations and  ₀ = ₀ - is the residuals of the pseudo observations.The covariance matrix of the estimated surface densities can be calculated after regularization based on the linear error propagation law as: where  = ₀⁻ ² ⁻¹ =  is the normal matrix and  =  + is the regularized normal matrix.The square root of the diagonal elements of   gives the theoretical standard deviation of the estimated densities, i.e. .
The ordinary Tikhonov regularization is a special case of Eq. ( 7) with  =  and ₀ = 0.This equation can also be used to define other regularization methods with zero-value pseudo observations ₀ = 0, or generally nonzero values ₀ if such prior information is available (from now on we assume ₀ = 0).Therefore, the matrix  in Eq. ( 5) or  in Eq. ( 6) can be used to define different regularization schemes.One of such methods is motivated as follows.
The patch's area is directly related to the variance of the estimated surface densities: the smaller the area is, the higher the estimated variance will be [Abedini et al., 2021].This makes sense because a smaller patch is based on a smaller number of observations and consequently the variance increases.A regularization method may take into consideration the areas of the patches.Therefore, as an alternative method to the ordinary Tikhonov regularization, the area-affected regularization has been proposed.This approach considers weights for different patches to accounts for the variable area among patches.The weight factor for each patch is defined as the inverse of patch area squared.The patch area for the surface element spanning from  to  +  and  to  +  on a spherical surface at radius  is computed as where  and  are the colatitude and the longitude, respectively.We denote the area of the surface patch  as, for  = 1, ... , .
Here, the area-affected regularization method is applied with three different scenarios of weightings.where  runs from 1 to ,  is the area of patch , and   is the average area over the entire 2° × 2° or 1° × 1° patches.The matrix ₀ introduces the Tikhonov regularization method.The matrix 1 has already been used by Abedini et al. [2021].Matrices 2 and 3 are proposed to somehow downweigh the areas' effect; 2 by applying a square root of 1, and 3 by its combination with the Tikhonov method.The effect of these regularization matrices on the results is discussed in Section 4.1.

Precision description of regularized solutions
Different weighting strategies of Eqs. ( 11)-( 14) are employed to obtain different regularized solutions.For quality description of the results, the precision of these regularized solutions is also of importance and hence the subject of discussion here.The variance ₀² of the original observations and the unknown regularization variance factor ² of the pseudo observations are a priori unknown and have to be determined to obtain a regularized solution.In principle, variance component estimation (VCE) methods should be utilized to obtain the two unknowns.Our results show that it is not possible to simultaneously estimate the two unknowns due to the Ill-posedness of the original problem.
We will therefore take a two-step procedure to reach this goal.In the first step, given an initial value for ₀², the regularization parameter  = ⁻² can be determined using the L-curve method or other methods numerically [Hansen, 2005], whereas in the second step, the variance ₀² is estimated using the least squares variance component estimation (LS-VCE).It is noted that, in practice, we may have more unknowns than only one variance component ₀² in .For example, combination of white and flicker noise in GRACE data has been reported by Guo et al. [2018].
Ignoring colored noise can result in an optimistic precision description of the estimated surface densities.This could be further investigated in a future study.

Equiangular mascon solution in GRACE-like missions
To implement the LS-VCE for estimation of ₀², after determining the regularization parameter , the covariance matrix  can be rewritten as ( 15) where ₀ is the known part of the covariance matrix and 1 is the only cofactor matrix of which its variance component should be estimated using LS-VCE [Amiri-Simkooei, 2007;Teunissen and Amiri-Simkooei, 2008].The least squares estimate of the variance component ₀² is (see Appendix) ( 16) where  =  - is the residuals of the original observations and is the trace of a matrix (sum of diagonal elements).
Because  and both  and  are functions of the unknown variance component ₀², the final  ₀² should be obtained in an iterative procedure.Experimental results show that only two iterations are sufficient to obtain a converged value.The standard deviation of the above variance component is given as here  is given from Appendix.This will then give ( 17) which provides a measure for the precision of the estimated variance component in Eq. ( 16).

Analysis of resolution matrix
It is well known that the regularized least-squares method provides biased estimates.The biasedness can be explained by a so-called 'resolution matrix' [Hansen, 2005].If we take the mathematical expectation (E) of the linear model  =  + , with E() = 0, we obtain E() = .Then the expectation of the least square solution in Eq. ( 7), with ₀ = 0, becomes which is not equal to , indicating a biased estimate of .The product  = ⁻¹ is called the resolution matrix.
When the design matrix A is regular, no regularization is required, and therefore it follows that  = , resulting  =  and E( ) = .In fact () =  is the Kronecker delta function ( = 1 if  =  and  = 0 if  ≠ ), defined on a discrete domain and a discrete analog of the Dirac delta function ().In such an ideal case, the transformation from  to E ( ) is accomplished without distortion/smoothing.This is what we know as an unbiased estimation.
Eq. ( 18), however, does not imply this.In this section we will further elaborate on this issue.

Singular value decomposition of resolution matrix
When the design matrix is ill-posed, a regularized solution is used.The singular value decomposition (SVD) of the normal matrix is [Teunissen, 1985] and [Hansen, 2005] where = [1, ... , ] is an  ×  matrix of eigenvectors of  and Λ is a diagonal matrix containing its eigenvalues ).When the regularization weight matrix is an identity matrix, This, with Eq. ( 20),  = ⁻ ¹ , and  =  = , gives where Ω = (Λ + )⁻¹ is also a diagonal matrix with entries ( 22) The fact that 0 ≤  ≤ 1 deviates from one makes distortions to the transformation from  to E( ).This distortion depends first of all on the magnitudes of the eigenvalues  and secondly on the regularization parameter .The smaller the eigenvalues are, or the heavier the regularization is (i.e. the larger the  is), the more significant distortion will be.One can show that the diagonal elements of the resolution matrix can be written as: (23) Because ₌₁ ² = 1`and 0 ≤  ≤ 1, it follows that 0 ≤  ≤ 1.If no or very soft regularization is required it follows that  ≅ 0 and hence  ≅ 1, which gives  ≅ 1.But if the value  is comparable to some of the eigenvalues , = 1, ... , ,  becomes (much) smaller than one.For the off-diagonal entries of  we have Since ₌₁  = 0 for  ≠ `and 0 ≤  ≤ 1, it follows that  ≠ 0. Again, it follows that when the value of  is comparable to some of the eigenvalues ,  = 1, ... ,, the values  can deviate much from zeros, indicating heavier regularization.
Having the above formulation, the i th entry of E( ) in Eq. ( 18) can be written as which with 0 ≤  ≤ 1 and  ≠ 0 verifies that parts of information in  will be lost, and simultaneously information from other ,  = 1, ... ,,  ≠  will affect , both indicating that E( ) is a biased estimate of .As a matter of fact indices  closer to  will affect  more than those that are far away.This will then have a smoothing effect on the estimate  .This smoothing effect is approximated using a sinc function in the next subsection.
For a geometrical interpretation of the regularization, further simplifications are assumed.If we partition  = [₁ ⫶ 2] and Ω blkdiag(Ω₁, Ω2), with blkdiag an operator that makes a block diagonal matrix from its input Ω₁, Ω2, it follows that (26) where ₁ and 2 are  ×  and  × ( -₁) matrices containing the first ₁ and the last  -₁ eigenvectors of .

Equiangular mascon solution in GRACE-like missions
where  ₁ = ₁₁ and  ₂ = ₂₂ are two orthogonal projectors [Teunissen, 1985].The actual  under  is projected to  as follows which consists of two vectors: ₁ =  ₁  is the recovered part, whereas ₂ =  ₂  is the missed or attenuated part because / tends to zero. Figure 1 shows the geometric interpretation of the above regularized solution.As can be seen, part of that belongs to the range space of U 2 , R(U 2 ), cannot be determined.Algebraically, this is referred to as under parameterization, which makes the estimate to be biased.To obtain more robust estimates, one may use a smaller regularization parameter indicating that the eigenvalues become larger.This is related to the size of the patches, addressed in the following subsection.

Patch-size criterion using sinc function
We propose a method to specify an optimal patch size (optimal discretization) for the mascon solution.By having the optimal patch size, our goal is to estimate densities on patches that reflect the gain (or resolving power) of the data and to avoid creating artificial densities that are not constrained by the available data and they are purely results of the interpolation imposed by the smoothing effect of the regularization.
When translated to the surface density changes (SDCs) estimated by the GRACE-like satellites, the above resolution matrix explains how information from one patch can affect/bias SDC estimates of the neighboring patches.Therefore, patch sizes play an important role in forming the eigenvalues  and hence the regularization parameter .The resolution matrix provides useful information to obtain an optimal patch size when estimating SDCs.It presents the combined effects of the mathematical model (through the design matrix ) and the a priori assumptions about the signal (through s 2 C) on the estimates of SDCs.
In one hand, the patches should be small enough to allow the correct reconstruction of the desired frequencies of the signal of interest.On the other hand, they should be big enough to contain sufficient and significant information distilled from observations at each estimated patch.
Abbas Abedini et al. theorem states the patches should not be larger than half the spatial wavelength of the highest-frequency component of the signal [Nyquist, 1928;Shannon, 1948].We therefore consider an upper limit for the size of the patches by computing the highest frequency of the estimated signal.Contrarily, there should also be a lower limit for the patch size.If the patches are too small, the estimated densities will be purely results of the regularizationimposed smoothing.The neighboring patches contain completely redundant information with almost equal (or highly correlated) estimates.It makes sense then to choose the patches sufficiently big to avoid the high correlation between adjacent patches.
To define these lower and upper limits, we analyze the entries of the resolution matrix of the applied regularization for both 1-degree and 2-degree solutions.We later show that the entries of resolution matrix follow  The sinc function provides information on the regularization because a sinc filter is equivalent to a low-pass filter with a cut-off frequency.In fact, the regularization removes the high frequency components of the signal and hence has a smoothing effect.The 2D sinc filter is written as where  and  show the easting and northing coordinates, and  and  are the cut-off frequencies in  and  directions respectively.This filter is the spatial equivalence of a 2D rectangular pulse in the frequency domain, which is obtained from the inverse Fourier transform of the following function: (31) Figure 3 shows a typical example of a 2D sinc function with B x = 1 and B y = 1.The cut-off frequencies can be interpreted as the highest signal frequencies that retained in the regularized solution.Based on the Nyquist criterion, the patches should be equal to half the spatial wavelength of the highest-frequency component of the signal.Such a patch size in the x and y directions is defined as  x = 1/2B x and  y = 1/2B y , respectively.The higher the cut-off frequencies and are, the smaller the spatial wavelengths  x and  y will be.In this study, by fitting a proper sinc function to the rows and columns of the resolution matrix, we can estimate the corresponding cut-off frequencies for different patches, and consequently we define the optimal patch sizes.

Results and discussions
The surface density changes are estimated for both the 2° × 2° and 1° × 1° equiangular grids over Greenland.The total numbers of patches are n = 406 and n = 1624 for the 2-degree and 1-degree grids, respectively.The simulated data contains 342 arcs from GRACE satellites that cross over Greenland region during one-month (see Figure4).
The total number of observations in this study is m = 77627.Therefore, we have two case studies namely Case 1 for 2-degree and Case 2 for 1-degree equiangular grids.The effect of different weighting strategies on the solution has been analyzed and presented in the following.

Different weighting strategies and number of arcs
The results from four different weighting strategies provided in section 2.3 have been compared based on a Monte-Carlo methodology.We have simulated independent Gaussian white noise vectors, with the standard deviation of  0 = 10 -8 m/s, and added them to the constructed exact observations (range rate vector).Each noise vector will then result in a new observation vector and therefore a new solution, which is referred to as an independent run.The independent runs will be repeated 20 times, leading to 20 independent estimates of the surface density changes.Different statistics can be derived using the estimates of the independent runs.
Using the weight matrices of the regularization expressed in Eqs. ( 11) to ( 14), the unknown density changes have been estimated for the two cases, followed by empirical investigations on the sample solutions.
To investigate the precision of each of the four weight matrices, the standard deviations of 20 independent runs are used for comparison.Here we use both the empirical standard deviation and also the analytical one derived from Eq. ( 17). Figure 5 shows the empirical standard deviations of the estimated surface densities resulted from different weight matrices for Case 1.The standard deviations have been ordered based on the patch id-number, starting from the south-west of AOI increasing to the north-east of the AOI.
Therefore, the small numbers are associated with the patches in the south, whereas the large numbers are associated with the patches in the north.The empirical and analytical standard deviations obtain the same results from different weights.To get a more illustrative comparison, we plot the empirical and analytical standard deviations of the 2-degree solution of the four weighting strategies in Figure 6 and Figure 7, respectively.
Abbas Abedini et al.Equiangular mascon solution in GRACE-like missions For some weight matrices, a clear increase in the standard deviations is observed towards the north pole for both figures (Figs. 6 and 7).This indicates that the northern cells have larger standard deviations as they have a smaller number of observations per cell.The results show that the weight matrix P 1 provides better results than P 2 in the southern part of Greenland.

Weight matrix
Closer to the north pole, P 0 and P 2 provide superior results compared to P 1 .The statistics for the 2-degree case have been presented in Table 1.As can be seen, the standard deviation resulted from the the weight matrix generally provides the worst results.The best performance is obtained for P 3 , indicating that a down-weighted area-dependent weight matrix is superior over other weight matrices.This superiority is evident in terms of maximum standard deviations, indicating homogenous precision description of the SDCs for P 3 (see Max STD and Std Dev STD results in Table 1).Therefore, the optimum weight is indeed P 3 , proposed to somehow down weight the areas' affected by its combination with the Tikhonov method.

Equiangular mascon solution in GRACE-like missions
The investigation of the standard deviation of Case 2 (i.e., 342 arcs and 1-degree grid) for four different weight matrices have also been done.Again, the empirical standard deviations have been presented in Table 2, and Figs. 8 and 9 shows the spatial distribution of empirical and analytical standard deviations, respectively.The results indicate that the solution of Case 2 is more precise than Case 1 solution.Here also, the results from P 3 weighting strategy shows better performance and higher precision.
Figure 9 shows the standard deviations obtained from the regularized ordinary Tikhonov solution (P 0 ) to emphasize the effect of the optimal area-affected weighting.It highlights that considering an appropriate areaaffected weighting has a significant combined effect in case of 1-degree compared to the 2-degree solution.
In the previous surveys, we investigated the precision of the above two cases under four weighting methods.An important aspect of an estimate is its closeness to its actual value, which is referred to as the accuracy of that estimate and indicates the bias of the estimate.The measure for that is the root mean square error (RMSE), computed from the difference between the estimated values and their true values.The results are presented in Table 3 for Cases 1 and 2. Two observations are highlighted.1) The smallest RMSE have been consistently obtained for the weight matrix P 3 , which is also consistent with the precision results presented in the previous section.
2) The RMSE values of the 1-degree cells are larger than those of the 2-degree cells.This is because heavier regularization is used for the 1-degree cells than the 2-degree cells, which was explained in this section.

Distortion analysis of regularized solution using SVD
This subsection explains how the regularized solution becomes biased using the SVD theory explained in Section 3.1.Figure 10 shows the results of the estimation in 2-degree and 1-degree grids for the Cases 1 and 2. In the right frame, the estimation results of SDCs are presented for the 2-degree patches.The frames on the left show a zoomin view for the 2-degree and 1-degree patches along with the positions of the satellites and the true values of the surface densities.The P 3 weighting strategy has been used to obtain the results presented in this section.Two observations can be highlighted as follows.
The first observation is that the estimated surface density changes are subject to distortions due to the regularization effect.These distortions are clearly visible due to some negative (underestimated) and overestimated values, which systematically differ from the true values.The highest negative values are in the cells close to the border areas where a sharp transition happens from land to ocean on the surface density changes.The negative sign indicates that the distortion has been negatively affected by the positive values of their neighbouring cells, to be expressed by a sinc function (see next subsection).The second observation is that the distortion becomes more severe for the 1-degree cells than the 2-degree cells.This is related to applying a heavier regularization, which is in conjunction with SVD analysis of section 3.1.
To further study the mechanism of the regularization, we assume two schemes in which the true values are assumed to be known.The two schemes are implemented on the 2-degree cells in the entire region of the Greenland.
The first scheme is presented in Figure 11-a in which the true values of Figure 10 has been used.
The second scheme assumes some known values on regularly gridded cells and zeros elsewhere (Figure 12-a).
The goal is to retrieve the true values using the resolution matrix decomposed by the SVD.In principle,  should be equal to  indicating no distortion.This is however not the case because the regularized solution consists of two vectors ₁ = ₁₁, which is a recovered part, and ₂ = ₂₂, which is a missed/attenuated part of ₂ = ₂₂.Total results for 2-degree patches (right), and zoom-in view of a typical area for 2-degree and 1-degree patches along with true values of surface densities and positions of satellite tracks (left).

Equiangular mascon solution in GRACE-like missions
The results of the first scheme are presented in Figure 11.As the recovered part of , Figure 11-b shows that the main body of the true values has been retrieved, which has a similar pattern to the true values in Figure 11-a.The results are also very similar to those estimated in Figs. 9 and 10 from the real measurements.There is also a second term ₂ = ₂₂, which ranges from -40 to 40 kg/m 2 (Figure 11-c).This term cannot be recovered indeed (Figure 11-d).This unrecovered part is mainly in the border areas and inside Greenland where a sharp transition in the surface densities has occurred.The positive and negative signs of neighbouring cells are also visible in this figure.
To further clarify this issue, the results for the second scheme are presented in Figure 12.A few observations are highlighted.Compared with the true values, the recovered part experiences less distortions at the borders than in the middle (Figure 12-b).The distortion in the east-west (EW) direction is worse than in the north-south (NS) direction.This is concerned with the spatial resolution in these two directions; the NS resolution is ~220 km, whereas the EW resolution is ~120 km in the southern part of Greenland.The EW distortion becomes worse when approaching the north pole, having a resolution of ~150 km in the northern part of Greenland.This is also clear in the Figure 12-c, which shows that the unrecovered part has an inverse pattern.Further, looking at a central cell, it is observed that its zero-value neighbouring cells get positive values, become negative, and finally tend to vanish, whereas their true values were assumed to be zeros.This behavior can be expressed as a sinc function, which is elaborated in the next subsection.

Patch size analysis using resolution matrix
As already observed, the maximum errors occur in the border areas where a sharp transition happens from land to ocean on the surface density changes.The signal oscillation at the border points is known as the Gibbs phenomenon [Gibbs, 1898[Gibbs, , 1899]].This needs further investigation in the present contribution using the sinc function introduced in section 3.2.This is in conjunction with the sinc function in which signal at a particular cell can be transferred to its neighboring cells introducing distortion.
The main question arises here as to what is the highest achievable resolution when using the mascon solution over Greenland given the characteristics of the GRACE mission.Accordingly, an appropriate patch size for the Greenland region can be proposed.Later on, we will refer to two minimum and maximum legible patch sizes in this respect.It is therefore important to know how large the patches should be to reliably estimate the surface density changes in the mascon method.To investigate the optimal patch size, the resolution matrix of the applied regularization provides a practical tool (see section 3.2).In the following, we analyze the entries of the resolution matrix to come-up with a criterion for the optimal patch size of the mascon solution.
We select a patch with the coordinates 72N, 44W at the center of Greenland in the 2-degree solution.Figure 13 shows the row entries of the resolution matrix that correspond to the selected patch.These values can be interpreted as the contribution of other patches to the estimate of the selected patch.On the other hand, information in the selected cell will also be distributed among the other cells.As can be seen, the contribution of the neighboring patches is more significant than those located far away.Further, a vanishing periodic behavior can be seen around the selected patch.The behavior is similar to the sinc function we presented in section 3.2.To further clarify this issue, we make a cross section along the longitudes and latitudes to make the relevant curves in the 2D space.The plots in the below and at the right of the main plot show respectively the EW and NE profiles, as indicated by the white dashed lines in the main plot.The resolution pattern of the main plot looks like a 2D sinc function, similar to Figure 3.The red curves in the plots show two 1D profiles of the fitted 2D sinc function to the entries of the resolution matrix.These 1D profiles clearly indicate the sinc pattern.The observed sinc pattern confirms that the applied regularization acts as a sinc filter that is equivalent to a low-pass filter with a cut-off frequency.In fact, the regularization removes the high frequency components of the signal and therefore has a smoothing effect.This is related to the small eigenvalues of the normal matrix that can significantly be distorted and therefore parts of information will be lost (see Figure 1).
As discussed earlier, based on the Nyquist criterion, the patches should be smaller than half the spatial wavelength of the highest-frequency component of the signal.An appropriate patch size (we call it maximum legible patch size) in the x and y directions is defined as  x = 1/2B x and  y = 1/2B y , respectively.If we fit the sinc functions h(x) = 2B x sinc (2pB x x) and h(y) = 2B y sinc (2pB y y) to the rows of the resolution matrix in the longitude and latitude profiles, we can estimate the cut-off frequencies of each patch.The curve fitting has been performed using a nonlinear least squares problem by a trust region reflective optimization method [Coleman andLi, 1994, 1996].For the selected patch, these estimates are equal to  x = 5.83 and  y = 1.86 degrees (red line in Figure 12).That is, to meet the Nyquist criterion at the selected cell, the patch size in the x and y directions should be smaller than 5.83 and 1.86 degree.These maximum legible patch sizes (intervals) have been indicated by solid black lines in Figure 12.
To determine the so-called minimum legible patch sizes, the rationale is to choose the patches sufficiently big to avoid high correlation (i.e., > 0.9) between adjacent patches.For this purpose, we use the auto-correlation sinc function in a given direction t (t can be x or y): corr(t) = . It can be verified that for t = 1/8B t , the autocorrelation is equal to 0.9.That is, to avoid the correlation larger that 0.9 between patches in x and y directions, the patches should not be smaller than  x = 1/8B x and  y = 1/8B y , respectively.These values indicate the minimum legible patch sizes and are equal to quarter of the maximum legible patch size.In the profiles of Figure 12, black dashed lines have indicated the minimum legible patch sizes in both x and y directions.We therefore obtained the minimum and maximum legible patch sizes by fitting a 2D sinc function to the row entries of the resolution matrix in the longitude and latitude profiles.
We applied the proposed resolution analysis to the entire patches in the Greenland.The cut-off frequencies in both directions were estimated using the nonlinear sinc function fitting.Using the above strategy, for each cell, the minimum and maximum legible patch size were determined.for the minimum legible patch sizes have been scaled by a factor of four to obtain those for the maximum legible patch sizes.As can be seen, the 1-degree and 2-degree solutions provide comparable results for the patch sizes.
The results show that the minimum and maximum legible patch sizes in the NS direction (latitudinal resolution) are almost constant and are equal to 0.5 and 2 degrees over the AOI.However, the resolution analysis in the EW sin(2) 2 Abbas Abedini et al.Equiangular mascon solution in GRACE-like missions direction (longitudinal resolution) shows that the minimum patch size in this direction increases at higher latitudes, which is due to the meridian convergence.
To better illustrate the latter effect, we present the fitted sinc pattern at six different latitudes ordered from south to north in Figure 15.As can be seen, the width of the sinc function increases at higher latitudes.These results suggest that latitude-adaptive grid sizes rather than equiangular grids should be used for mascon solutions in highlatitude regions.One way to handle this problem is to use cells that are more or less regular.There are several exact solutions for the division of the sphere surface into regular spherical polygons.To approximate the sphere into a regular polyhedron, one possibility is to use the tetrahedron, which divides the area of interest into an imposed number of equal-area cells.Under such a condition, the cells can have very comparable forms, e.g.close to the square.There are a few general rules to reach this goal.For further information, we may refer to the study by Beckers and Beckers [2012].
The results suggest that patch sizes should not be smaller than 0.5 and 1.5 degrees, respectively in the NS and EW directions in the Greenland area.In the east-west direction, the cell size should increase to 4 degrees.We however note that because of the smoothing effect induced by the regularization, it is, in principle, possible to derive a solution even for patches smaller than the above values.But the obtained solution would seriously be affected/biased due to the regularization.Under such a condition, the majority of the SVD eigenvalues becomes smaller than the regularization parameter and therefore the distortion illustrated in Figure 1 becomes significant.An important question here is whether or not the use of the area-weighting strategy is still relevant if we use the proposed resolution-based method to define the patch sizes.We note that the proposed resolution method does not necessarily result in an equal-area grid.In fact, this method does not suggest a fix size for the patches, but a range of legible patch sizes in the latitude/longitude directions.One can somehow design an equal-area grid whose patches have a legible size.In this case, the proposed area-weighting method does not have any added value.In general, to answer the above question, we may first design a grid using the proposed resolution analysis and then calculate the area of its patches.If they have almost equal areas, the area-weighting strategy is not necessary anymore.Otherwise the area-weighting method can be effectively used.

Summary and conclusions
In this paper, we addressed two practical issues regarding the GRACE-based mascon estimation using local basis functions.First, we investigated four different area-affected weighting strategies for the mascon solution.The results confirmed that using the area-affected weights can improve and homogenize the quality of the solution.Also, based on the comparative analysis between the weighting strategies, we conclude that the P 3 weighting method, based on the square of each patch area plus mean area, is the optimal weighting strategy.
Secondly, we introduced a methodology to analyse the resolution matrix of the mascon solutions.The proposed methodology is twofold.First, it helps to find out how regularization can bias the mascon solution using the SVD decomposition of the resolution matrix, and second what would be the achievable resolution in the mascon solution over the Greenland region given the characteristics of GRACE mission.Note that, although we applied this method over Greenland and for mascon estimation, the proposed framework is rather general and can be used for any regularization problem.Based on the results of the resolution study we conclude that the minimum legible patch size in the Greenland area for the current setting of GRACE observations is 0.5 degree in NS direction.We therefore propose to use latitude-adaptive grid-size rather than equiangular grids for mascon solutions in high latitude regions.This means that in higher latitudes the bigger patches should be used.
There are likely trade-offs between the latitude-adaptive patch size using the grid-size selection criterion and the equiangular grids adjusted by the weight matrices.This requires further research in a future study in which the performance of these two methods on the precision and efficiency of the estimated surface density changes can be investigated.

Figure 1 .
Figure 1.Geometric interpretation of regularized solution consisting of two vectors 1 = 11, as a recovered part, and 2 = 22, as a missed/attenuated part.

Figure 2 .
Figure 2. The concepts of minimum and maximum legible patch size.

Figure 5 .
Figure 5.The Empirical standard deviations (kg/m 2 ) of estimated surface density changes regularized by four weight matrices for Case 1 (342 arcs and 2-degree grid).

Figure 4 .Figure 7 .
Figure 4.The arcs of the satellites' path are crossing over the AOI (Greenland region)

Figure 8 .
Figure 8.The Empirical standard deviations (kg/m 2 ) of estimated surface density changes regularized by four weight matrices for Case 2 (342 arcs and 1-degree grid).

Table 2 .
Different statistics in terms of mean, minimum, maximum and standard deviation (over n = 1624 patches) of the empirical standard deviation (STD) of estimated surface density changes over Greenland derived from 20 independent runs for Case 2. Units are in kg/m 2 .
Figure 13 visualizes the minimum legible patch sizes resulting from the resolution analysis applied to both 1-degree and 2-degree solutions in the NS and EW directions.The results for the maximum legible patch sizes are not presented because they have the same pattern as the minimum legible values; their sizes are a scaled (by a factor of four) version of the minimum legible patch sizes.Two profiles in the NS and EW directions of Figure 12 are presented in Figure 14a and 14b.Again, the results

Figure 12 .
Figure 12.Resolution analysis over a patch at the centre of the Greenland to obtain the maximum and minimum legible patch size.

Figure 13 .
Figure 13.Critical patch sizes obtained from the resolution analysis applied on both 1-degree solution (top a,b) and 2degree solution (bottom c,d) in both directions.

Figure 14 .
Figure 14.Minimum and Maximum legible patch sizes for a) north-south direction, b) east-west direction.

Figure 15 .
Figure 15.Fitted sinc pattern at six different latitudes (ordered from south to north) for 2-degree solution.

Table 1 .
Different statistics in terms of mean, minimum, maximum and standard deviation (over n = 406 patches) of the empirical standard deviation (STD) of estimated surface density changes over Greenland derived from 20 independent runs for Case 1. Units are in kg/m 2 .

Table 3 .
Root mean square error (RMSE) of estimated surface density changes over Greenland derived from 20 independent runs for two Cases and different weight matrices.Units are in kg/m 2 .