The Retrieval of Wet Refractivity Index by Tomography using Spherical Cap Harmonics

High spatial and temporal variability of the tropospheric wet refractivity index, makes it difficult to present an accurate model for this variable. Up to now, Radiosonde stations data have been used for monitoring atmosphere parameters. Furthermore, because of the sparse distribution of radiosonde stations to monitor the lower levels of the atmosphere, the numerical weather models do not have enough accuracies for atmospheric parameters. Using the GPS tropospheric delay measurements and tomography approaches, the wet refractivity index can be estimated. In this research, three-dimensional and four-dimensional basis-function tomography is used to demonstrate the distribution of wet refractivity index of the troposphere. In this model, spherical cap harmonics are used for the horizontal distribution of the wet refractivity index, and empirical orthogonal functions are used for the vertical distribution of the index. In addition, temporal changes are considered by correlating the unknown coefficients using fourier series. The region of study is in the west California State, and the wet refractivity index is retrieved from the wet tropospheric delay measurements. To validate the results, radiosonde profiles were compared to the tomographically retrieved profiles. The results show that wet refractivity indices can be retrieved using functional models with RMSE about 2.4 ppm till 3.9 in the four-dimensional method. The comparisons show that the four-dimensional retrieved profiles show improvement up to 34 and 42 percentages in mid-day tomography epochs compared to the three-dimensional tomography results. Also it can be seen that in mid-night epochs, the three-dimensional tomography has higher accuracy compared to four-dimensional method because of low variation of wet refractivity indices.


Introduction
Due to the variability in the wet refractivity index in the atmosphere, a microwave signal experiences various velocities in its path through the troposphere. This variability in the wet refractivity index is because of the variation in humidity and water vapor content in the troposphere. Water vapor is an important constituent of the atmosphere, the knowledge of which has a significant role in weather forecasting [1]. Because of the lack of information about the water vapor distribution in the Vol. 15 / Special Issue 2022 / (51) atmosphere, numerical weather models suffer from low accuracies [2].
The Global Navigation Satellite Systems (GNSS) signals are delayed by the water vapor content in their path through the atmosphere. This delay is divided into wet and dry components. dry component can be approximated from the temperature and pressure at the station. The wet component is largely dependent on the water vapor content of the atmosphere. The tropospheric dry delay can be measured by weather information at a station and using models such as Saastamoinen [3] or Hopfield [4]. The wet tropospheric delay can be estimated from GNSS observations using precise positioning techniques [5]. The wet tropospheric delay measurements contain information about the structure of the wet refractivity index in the atmosphere; therefore, spatial and temporal structure of the wet refractivity index can be estimated from the wet tropospheric delay measurements in a dense network of GNSS stations using a tomographic approach [6]. Tregoning [7] evaluated the capability of Global Positioning System (GPS) for estimating the water vapor and other tropospheric components by comparing the estimates with radiosonde measurements. Bevis [8] and Emardson [3] compared water vapor estimated from GPS observations with radiometer data. These studies only determined the amount of water vapor in the atmosphere, but the structure and distribution of the water vapor was not estimated. In further studies, the distribution of the tropospheric parameters was estimated by tomographic approaches in regional networks. For instance, Flores [9] implemented a tomographic approach to estimate the water vapor in Hawaii region. Due to the inhomogeneous distribution of the stations used in a network, and lack of GNSS signals passing through the atmosphere, and also because of the observational errors, a tomography problem is usually an irregular and ill-posed problem. Several methods have been presented for regularization of the problem. For example, Hirahara [10] added some constraints in horizontal and vertical directions by means of non-iterative regularization methods. Nilsson [11] showed that the vertical estimation accuracy will be significantly improved if the altitude difference of the network stations is in the range of 800 meters. Bender [6] compared the convergence time of several iterative regularization approaches in a study aimed at investigating climate change in Germany. In a more recent research, Ding [12] exploited the combination of iterative and non-iterative approaches for the regularization. Xia [13] investigated the optimum number of iterations for Algebraic Reconstruction Techniques (ART) and IART (improved ART) methods. Because of the irregularity and ill-posedness of the tomographic approaches, using constraints or a priori values for the unknowns is inevitable in all the above studies. Knowledge of proper constraints and adequate priori information about the unknown parameters are limited to the iterative and non-iterative regularization approaches. To avoid this limitation, unknown parameters can be expanded to orthogonal basis functions. This way, the number of unknowns which are only the coefficients of the orthogonal basis functions, are far less than the actual unknowns (i.e. the wet refractivity indices). An example of these basis functions is the spherical harmonics representing a global field. This basis function is appropriate for representing a field in a global level scale. To accurately represent a field in a regional scale, the degree of spherical harmonics, i.e. the number of unknown coefficients, must be increased. Increasing the number of unknown parameters not only increases the computational time required but also makes the problem more irregular. Haines [14] applied some changes to spherical harmonics and utilized them to model a field in a region of sphere. He named them Spherical Cap Harmonics (SCP). [15] [16] used the SCHA (Spherical Cap Harmonics Analysis) to estimate the horizontal distribution of the electron density in the ionosphere and used empirical orthogonal functions for vertical distribution of the electron density. They presented a three-dimensional model for the electron density parameter. Farzaneh [17] utilized spline basis functions and empirical orthogonal functions to recover the threedimensional structure of the ionosphere. Schmidt [18] used B-splines to correlate the unknown coefficients of the spherical harmonics in the time, and estimated timevariable coefficients. In this study the observations from a set of CRTN (California Real Time Network) stations located in the west of the United States are used for implementation of 4D and 3D functional tomography method to determine the distribution of wet tropospheric refractivity index in the network of the study. It is demonstrated that the four-dimensional tomography approach yields more accurate estimates than the threedimensional approach in the mid-day epochs. To validate the results, radiosonde observations in the summer and winter are utilized in order to identify the seasonal effects on the functional model parameters.

Methodology
At first, the tomography method for reconstruction of atmospheric wet refractivity indices will be described. Then, implementation of SCH for derivind 3D and 4D models of wet refractivity profiles will be illustrated and results will be compared to radiosonde data.

Tomography
The equation of GPS tropospheric delay measurements is [2] : 10 (1) Where STD is the tropospheric slant total delay, SHD the tropospheric slant hydrostatic (dry) delay, SWD the tropospheric wet delay, L is the GPS signal length between the satellite and the receiver, the dry refractivity index, and is the wet refractivity index. Using weather information, refractivity indices are estimated [19]:

77.6890
(2) 71.2952 375463 Where and are partial pressure of dry weather and partial pressure of water vapor in terms of hPa, respectively. T is the temperature in Kelvin. considering a consistent value for wet refractivity indices in each voxel, the slant wet delay can be rewritten as follows [6]: Where , , and , , are the corresponding wet tropospheric refractivity and partial length of signal through each voxel. Providing all the wet tropospheric delay measurements, the observation equation is: Where Y contains the observation of wet tropospheric delay, A Consists of the partial derivative of the slant wet delay with respect to the wet refractivity indices of each voxel in order of , m is the number of observations, and n is the number of unknown parameters (i.e. wet refractivity indices). N is the vector of wet refractivity indices (i.e. the unknown parameters) of order 1. The columns of the design matrix are lengths of every signal passing through the specified voxel. Due to the large number of unknowns and the sparse distribution of observations which prevents the signals from passing through some of the voxels, Equation (5) is an ill-posed equation system. The illposedness can be depicted by the condition number as [20]: (6) Where λ and λ are the maximum and the minimum of the singular values of the design matrix. The Picard condition is not satisfied, if the condition number is large; and the least squares solution of unknowns cannot be obtained confidently [20]. Thus, to estimate a regular and unique solution of unknowns, a proper regularization approach must be used.

Spherical Cap Harmonics
The spherical harmonics expansion of a f λ, θ function on a sphere is expressed as: f λ, θ a cos mλ b sin mλ .
P cosθ Where a and b are the coefficients of spherical harmonics expansion, and P cosθ represents the Legendre function. If the region of study is part of a sphere, the spherical harmonic functions are not suitable bases for expansion of a function in the specified region, because the Legendre functions are orthogonal only on the surface of a sphere [14]. Haines [14] suggested using Spherical Cap Harmonics (SCH) to overcome this problem. In his suggested approach, an modified Legendre function with a non-integer degree n are used as bases to expand the functions of cap shaped region in part of a sphere. These Legendre functions and their corresponding gradients tends to the zero at edge of the cap-shaped area with half-angle [14] d θ dθ 0 , for k m even The value of the n can be determined using the above equation. Substracting modified Legendre functions, f(λ,θ) can be written as follows: Where k is the maximum degree of the SCH that represents the best reconstruction of the considered field.

Coordinate transformation
Translating the coordinates to a new coordinate system (in which the origin and pole are the center of the earth and the center of the spherical cap, respectively) is the first step of expanding a function to the spherical cap harmonics. To transform a geographic coordinate system into the spherical cap harmonics system, Equation (11) and (12) are used [16]: cos θ cos θ cos θ sin θ sin θ cos λ λ (11) tan Π λc sin θ cos λ λ sin λ cos θ cos θ sin θ cos λ λ

Computation of the Associated Legendre Function
The modified Legendre function with non-integer degree n can be written as: Details for computing A m, n can be found in [17]. The derivative of Equation (13) can be written as: Equations (8) and (9) can be solved by Equations (13) and (14); then, the non-integer n can be computed for different m and k. To estimate the horizontal distribution of the refractivity index, the unknown parameters of Equation (5) can be expanded into spherical cap harmonics as follows: λ, θ a cos mλ b sin mλ . cosθ (15) Equations (5) can be rearranged as: (16) Where A is a n N matrix that consists of the bese functions of spherical cap harmonics in which N is the number of spherical cap harmonics expansion coefficients and equals to k 1 . The matrix x contains the expansion coefficients.

Empirical Orthogonal Functions (EOF)
Spherical cap harmonic functions show the horizontal distribution of the wet refractivity index, but cannot indicate the vertical distribution of the refractivity index. To overcome this problem, empirical orthogonal functions forming a vertical basis of the space are used to represent the vertical distribution of the wet refractivity index [16]. The empirical orthogonal functions are obtained from the observations related to the desired quantity; therefore, to generate these functions, the refractivity index vertical profiles are required at different times . For this purpose, ERA-Interim numerical weather prediction model data have been used to calculate the wet refractivity index in the network of study. ERA-Interim is an global atmospheric reanalysis model that provides the atmospheric parameters since 1979. this model is grided data with spatial resolution of about 79 km in 69 different pressure levels. In this numerical model, data set of the temperature and humadity data at different pressure levels can be found. Therefore, with this data, wet refractivity indices can be calculated for different times and locations. Assuming the computed refractivity indices at different times t i 1,2, … , M and heights h j 1,2, … , N are indicated by N t , h , the matrix of wet refractivity indices can be written as [15]: In the above matrix, the columns are the time series of the refractivity index, and the rows are the refractivity indices at different heights. In matrix N t, h , by subtracting the mean value of each column from its elements, the mean value of the corresponding column will be equal to zero. Expressing the obtained matrix as N t , h , the covariance matrix can be determined as follows: In Equation (18), the eigenvectors of the matrix S are the empirical orthogonal functions. The covariance matrix indicates the behavior of random variables. The larger the amount of an eigenvalue, a more significant role the corresponding empirical orthogonal functions play in the system behavior. The other empirical orthogonal functions corresponding to small eigen values are due to the random errors in the data [21]. The vertical distribution function of the refractivity index can be written as below using the empirical orthogonal functions [21]: (19) Where α are the expansion coefficients of the vertical function. By combining Equations (15) and (19)

Fourier series
To correlate the wet refractivity indices in the time, the Fourier series can be used. The h(x) function can be expanded as below using the fourier series [1]: Where are the unknown coefficients correlated in the time, and are the basis functions. For the regional signal modeling, one can write , where are the fourier series base functions. Thus: Where each coefficient represents two individual unknowns. The fourier series base functions are obtained by the following equation [1]: To estimate the unknowns correlated in the time using fourier series, Equation (16) can be rewritten as: where q 1 , 2 , … , Q are the observation epochs, i 1 , 2, … , I are the number of observations at each epoch, a is a matrix including the spherical cap harmonics and empirical orthogonal functions, x is the unknown coefficient of spherical cap harmonics expansion and the coefficient of the empirical orthogonal functions at time q. Therefore, for all the observations at different epochs, Equation (25) is rewritten as [18]: Where R is the observation error, and A is the coefficients matrix as A a , a , … , a . By showing the unknown parameters of the matrix X as c , t , these coefficints can be expanded using the fourier series basis functions [18]: Equation (27) is expressed as below for all the unknown coefficients: (28) Letting N Q K 1 , the matrix C will have the dimensions , and is expressed: The vector u with dimensions k 1 is expressed as: ɸ , ɸ , . . . , ɸ Substituting Equation (30) into Equation (26) yields: Where Y y , y , … , y is the observation matrix with dimensions I Q, E e , e , … , e is the error matrix, and U u , u , … , u has dimensions k Q. To compute the unknowns, Equation (34) should be rewritten using a tensor product approach [18]: Recognizing X U ⊗ A and defining β vecC yields: Where P is the weight matrix of the observations. It is assumed that the observations are independent and have equal weight. The time-correlated unknown coefficients can be estimated utilizing Equation (33). Since Equation (33) is an ill-posed system, regularization methods need to be used. To regularize the problem and to solve the Equation (33), Algebraic Reconstruction

Vol. 15 / Special Issue 2022 / (51)
Technique (ART) can be used [1]. This method uses the confluence of the assigned hyper-planes to each observation in order to estimate the unknowns in an iterative procedure. The iteration number of the regularization methods were discussed by [13]. They showed that the number of optimum iterations in ART approach is 71 and increasing the iteration number makes the solution semi-convergent and biased.  (4) shows The altitude distribution of the GPS stations in the network vertical resolution. The difference between the maximum and the minimum station height is about 2400 meters. The lowest altitude of stations in the network is about 27 meters below sea level; thus, the minimum height of the network is assumed 50 meters below sea level. in altitudes higher than 9 km, the value of the wet refractivity indices tends to zero. So, the maximum altitude is considered to be about 10 km [1]. Therefore, the minimum and the maximum heights studied are -50 m and 9950 m, and the vertical spacing between the voxels is chosen to be 500 meters. The area of the study region 210 km × 210 km and vertical resolution of voxels is considered as 40 kilometers. Therefore, based on the network gridding choices, the number of unknown parameters of the tomography solution (i.e., the number of voxels in the network) is 500. As shown in Figure (3), there are 45 GPS stations in the network; the horizontal network has 17 voxels in which at least one GPS station exists. The only radiosonde station in the network is NKX which is located 21 km from the nearest GPS station and can be seen in figure (3). The KNKX weather station is also located near the radiosonde station, and can be used to measure weather parameters at the NKX station. Gamit and GLOBK programs have been used to estimate SWD (Slant Wet Delay) that are input observation for tomography approach. The values of the SWD at 30 second intervals are estimated from the GPS precise positioning solution. ZTD (Zenit Total Delay) has been estimated at two-hour intervals. The Global Pressure and Temperature (GPT) empirical model is used for the priori information to estimate ZHD (Zenit Hydrostatic Delay) and extract ZWD (Zenit Wet Delay) for ZTD and the Global Mapping Functions (GMF) are used for the mapping functions of the hydrostatic and wet tropospheric delay components.

Tomography region and GPS data
Tropospheric horizontal gradients are also applied in the least-squares solution. The International GNSS Service (IGS) final orbit products with sampling of about 15 minutes are used for positioning and interpolation of the GPS satellite positions at 30-second time intervals that is necessary for computing the passing GPS signal length at different voxels. The observation of about onehour are considered for estimation of the tomography problem in a manner that the unknowns will be placed at the last epoch of each one-hour interval. The EOFs are constructed from the ERA-Interim data. The obtained wet refractivity indices will be compared with the radiosonde data using RMSE (Root Mean Squared Error) as follows [13]: Where N is the number of unknowns in the retrieved profile, N is the wet refractivity index estimated from the tomography solution in voxel i, and N is the corresponding radiosonde wet refractivity. To regularize the problem and to solve Equation (33), the ART approach is employed .

Results and Discussion
In this study, a functional model is used to estimate the wet refractivity indices in three and four dimensions. For this purpose, the mathematical model is developed at two observation epochs in summer and winter. The computed wet refractivity indices have been compared with radiosonde data for the estimation of the optimum number of the SCH coefficients and the number of the EOFs in Equation (20). In addition, the order of fourier series in Equation (27) is assessed in different status and different seasons. The radiosonde data are only observed twice a day at 12 and 24 o'clock in UTC. The observation epochs are therefore considered at 12 and 24 o'clock on 2 nd July 2011 and 8 th February 2012 to study the impact of the different seasons on the parameters that must be computed in an optimal sense. From the radiosonde data, precipitable water values in millimeter for tomography epochs are listed in table (1). Information of the wind speed on earth's surface can be collected from KNX weather station that is near the radiosonde station. This information has been shown in table (2).
The estimated wet refractivity indices are assessed in three and four dimensions. The number of EOFs in the Equation (20) was considered to be 1, because about 99% of the whole data was retrieved by the first EOF. One epoch of the radiosonde data was selected from each season to determine the optimum degree of the parameter in equation (20) and (27) for three and four-dimentian models. The considered epochs are 12 o'clock for 2 nd July 2011 and 8 th February 2012. To evaluate the optimal model parameter SCH degree in equation (20), the order of fourier series (kj=m) in equation (27) needs to be determined. First, to determine optimal fourier series degree in equation (27) different values are considered for SCH degree and order of fourier series, then, for different orders of fourier series and the considered parameter for SCH degree, the wet refractivity indices are estimated using functional model. Then, RMSE for estimated wet refractivity indices compared to radiosonde data has been calculated. The RMSE values have been computed for different considered parameters to indicate the best set of functional model parameters at each considered epoch. The results are shown in Figure (5) and the best model parameters are displayed in Table (3). As illustrated in Figure (5), the minimum value of RMSE corresponds to k 2 and for 2 nd July 2011 and k 2 and for 8 th February 2012 epochs. Thus, the optimum model parameters for the two epochs are k 2, , and EOF 1. Based on the estimated parameters, the solution of the tomography problem in 3D and also 4D methods can be obtained. By implementation of Equation (20), results for 3D functional model are shown in Figure (6) and Table (4 8 3.6 Table 3. the best functional model parameters   8 3.30 Also, the results for four-dimensional model can be shown in Figure ( (4) and (5), the accuracies are improved by the four-dimensional approach compared to the threedimensional approach for the 12 o'clock epochs, but these improvements are less for 24 o'clock epochs. The improvement in the accuracies at 12 o'clock can be explained by the fact that since the tomography time intervals are 1-hour windows and the wet refractivity changes related to the water vapor are higher in the midday than the midnight as the studied region is a coastal area and the higher temperatures at midday, correlating the unknowns in time improves the accuracy of the estimations. For the 24 o'clock epochs, because the water vapor (and correspondingly the wet refractivity index) is less variable between the 1-hour tomography time intervals at midnight, correlating the unknowns in time does not improve the accuracy compared to the midday epochs. From figure 7( (2), this epoch has the highest wind speed. Thus, the difference between the radiosonde profile and the tomography solution can be related to the distance of radiosonde balloon from NKX station. According to the weather information of the wind speed in the station KNKX in Table (2), the wind speed had the minimum value on 2 nd July 2011 at 12 o'clock. As shown in Figures 5(a) and 7(a), the difference between the computed wet refractivity index profile and the radiosonde data in the lowest altitude is small in both three and four dimensions for this epoch. The difference between the tomography solution and the radiosonde data is generally larger in the lower levels of the atmosphere than in the upper layers both in three and four dimensional cases. The reason could be that the variability of the wet refractivity index is higher in lower levels than the higher levels of the atmosphere, and thus, the estimated solution has larger errors.

Conclusion
In this research, the tomography functional models were developed to retrieve the atmospheric refractivity indices. The ill-posed tomography problem was solved using functional model that led to reduction of unknown parameters. The wet refractivity index profiles were eventually estimated. For this purpose, a fusion of the base functions from spherical cap harmonic and the empirical orthogonal was performed to indicate the unknowns in three dimensions and afterward, Fourier series have used to expand the unknowns in the time domain and represent the four-dimensional model. The spherical cap harmonics were used to represent the longitudinal and latitudinal distribution of the unknowns. The empirical orthogonal functions generated from weather information in the network were used to show the height variability, and the fourier series were used to represent the unknowns over time. The results in both three and four dimensions were compared with the radiosonde station data in the network to evaluate the accuracy of the proposed approach. Based on the results, the parameters of the functional model are different in various seasons. This difference can be caused by the water vapor density changes in the atmosphere; the higher the density, the higher numbers of parameters are needed to estimate the unknowns. Due to the higher variability of the wet refractivity index in the midday, the four-dimensional model, which correlates the unknown in time, performs more accurately than the three-dimensional model at midday. It was demonstrated in this study that the wet refractivity indices can be estimated by developing the functional model and decreasing the number of unknowns.