Machine Learning-Based Estimation of Hourly GNSS Precipitable Water Vapour

: Water vapour plays a key role in long-term climate studies and short-term weather forecasting. Therefore, to understand atmospheric variations, it is crucial to observe water vapour and its spatial distribution. In the current era, Global Navigation Satellite Systems (GNSS) are widely used to monitor this critical atmospheric component because GNSS signals pass through the atmosphere, allowing us to estimate water vapour at various locations and times. The amount of precipi-table water vapour ( 𝑃𝑃𝑃𝑃𝑉𝑉 ) is one of the most fascinating quantities, which provides meteorologists and climate scientists with valuable information. However, calculating


Introduction
Water vapour is a crucial component of the earth's atmosphere, and it plays a significant role in climate and weather systems.As a result, monitoring this variable atmospheric greenhouse gas has a direct impact on short-term and long-term studies of the weather, including climate change and weather forecast.Precipitable Water Vapour () is a valuable water vapour product defined as the vertically integrated atmospheric water vapour in a column of a unit area [1,2] and can be estimated with conventional techniques, such as radiosonde measurements, water vapour radiometer, or data derived from numerical weather models.There are still some drawbacks to these techniques, such as their high cost and poor spatial-temporal resolution, even though they have many advantages.These limitations can be mitigated by GNSS due to its continuous scanning of the troposphere at a low cost with a higher spatial-temporal resolution in all-weather conditions [1,[3][4][5].Thereby, GNSS observations can be used to estimate  with long-term stability, high accuracy, and high spatiotemporal resolution due to the development of GNSS technology and the increasing number of permanently tracking GNSS receivers [1,[6][7][8][9][10].
Calculating GNSS  requires estimation of Zenith Wet Delay ().To do so, GNSS measurements are first processed with the GNSS software (here, Bernese version 5.2 [11]) to estimate Zenith Tropospheric Delay ().Then, the Zenith Hydrostatic Delay (  ) can be calculated by utilizing well-known hydrostatic models such as Saastamoinen by introducing the required meteorological parameters [12].In the next step,  is obtained by subtracting  from  , and then  is converted to  using a conversion factor (Π) [1,13].Since this factor strongly depends on weighted mean temperature (  ),  estimation is affected directly by   , whose value differs depending on the geographical location (, , ℎ), season, and weather condition [7,14,15].
One of the earliest and most influential works on GPS meteorology was by Bevis et al. [16] which paved the way for further research on GNSS-based  estimation and modelling.The authors found that GPS data provided accurate and continuous measurements of  that could potentially be used to enhance the Numerical Weather Model accuracy, predicting severe rainstorms and studying the climate [8,[17][18][19][20][21][22][23][24].Based on this exemplary study, some other researchers have done valuable work in the past three decades to model  with different methods at various spatial scales and temporal resolutions.In 2019, Zhang et al. [25] developed a real-time GNSS  monitoring system utilizing the Chinese national GNSS network.Real-time GNSS  were generated with high accuracy and performance using the precise point positioning technique.A validation of the GPS  using radiosonde data found the system to be highly accurate, with a mean bias of −0.1 mm and root mean square of 1.7 mm.A method for analysing  using a low-cost GNSS receiver installed onboard the ISABU vessel, operated by the Korea Institute of Ocean Science and Technology, was developed by Sohn et al. 2020 [26].Radiosonde data and GNSS  were found to be in close agreement in their validation.In spite of their results, shipborne GNSS offers the potential for accurate  derivation via kinematic precise point positioning.In 2022, Li et al. [27] used multiple linear regression to establish multi-factor  conversion models, focusing on the correlations between , , surface temperature, and atmospheric pressure.Their results indicated that the multi-factor models significantly enhanced the models' robustness and applicability across the China region.
Some researchers have also explored the use of machine learning methods to model .For example, a generalized regression neural network was used by Zhang and Yao in 2021 [28] to fuse  from GNSS, Moderate-Resolution Imaging Spectroradiometer (MODIS), and European Centre for Medium-Range Weather Forecasts Reanalysis 5 (ERA5) in North America.In the proposed method, sparse high-quality data could be used to enrich low-quality data to enhance their quality and utility.Consequently, a homogeneous  product can be derived with reduced temporal and spatial variations in accuracy, thus raising data usability for weather and climate monitoring.A stacked machine learning model was proposed in 2022 by Zheng et al. [29] to map  into  without meteorological parameters at the global scale.To validate the  values, ERA5 and radiosondes measurements were utilized, which showed a consistency within 2.5 mm.According to the study results, time-critical meteorological applications can benefit from this model, and other methods that are capable of sensing  can take advantage of it for determining water vapour in real -time.
As already stated,   is a crucial parameter in  retrieval.In general,   modelling methods can be categorized into surface meteorological factor models and non-meteorological factor models based on the use of surface meteorological data [30,31].In the first method, a linear or nonlinear model of   is developed based on the relationships between this parameter and surface meteorological variables, such as surface temperature (  ), pressure (  ), and water vapour pressure (  ) [31].The seminal study by Bevis et al. (1992) [1] pioneered establishing the linear model between   and   , which has been the basis for water vapour modelling studies since then, see for instance [7,[32][33][34][35][36][37].In addition, this technique has also been used in numerous studies to obtain a more accurate global or local model of   [19,[37][38][39][40][41].For the second method, instead of using in situ meteorological data, other variables such as geographical location and time are taken into account to model   [30,31,42].There have been a number of empirical models developed in recent years, including different models of Global Pressure and Temperature (GPT) [43][44][45], Global Weighted Mean Temperature [46,47], Global Tropospheric Model [48], Hourly Global Pressure and Temperature model [49], GTm_R [34] model, and the Global Weighted Mean Temperature model [50].
In the following, we summarize some recent studies on modelling   on local and global scales.In 2017, Manandhar et al. [35] proposed a simplified model for the conversion factor for retrieval of  from GPS signals based on latitude and day of year.The proposed model was compared with temperature-dependent models using data from 174 stations. values derived from the proposed model were validated with temperature-dependent models using three databases.Based on the obtained results, the simplified model can be applied universally, since it is computationally efficient, and has high accuracy for a wide range of geodetic applications.In 2020, Wan et al. [51] utilized radiosonde profiles from 12 Antarctic and 58 Arctic stations between 2008 and 2015 to establish a model relating   and   in polar regions.Two regional   models were developed, one based on linear regression and the other based on the quadratic function.These models were found to be more accurate than the global GPT2w (GPT 2 wet) [44] model.Even though the quadratic function   model had slightly a higher accuracy, both could be used for retrieving s.In 2021, Long, Hu et al. [42] established   models over China using the neural network technique with three different schemes, namely, a non-meteorological-factor   model, a single-meteorological-factor   model and a multi-meteorological-factor   model.In this study, the three-layer feedforward neural network method was used with the help of an ensemble learning in order to combine multiple models and consequently improve the accuracy and robustness of the predictions.The new models were found to be very capable of capturing regional spatiotemporal variations in   and simulating the interactions between   and a variety of surface meteorological parameters.Zhao et al. [52] developed a global conversion factor model of nonmeteorological parameters using gridded data taken from the 2006 to 2014 Global Geodetic Observing System.In addition to the geographical location, the changes in atmospheric water vapour on an annual basis, semiannually, and day by day are taken into account in this model.Based on European Centre for Medium-Range Weather Forecasts Reanalysis (ERA5) reanalysis data from 1990 to 2018, Zhang et al. [31] proposed a principal component analysis method for modelling   in Greenland where in situ meteorological data are limited.Their model validated using 11 radiosonde stations (2015-2019) showed a bias of −0.110K and root mean square error of 4.447 K. Additionally, the principal component analysis model performed better than GPT3 and Global Tropospheric Model, required fewer parameters, and improved computation efficiency.
Since Artificial Neural Networks (ANNs) and Random Forest (RF) are well-established machine learning algorithms that have the potential to learn complex relationships between input and output variables and can handle high-dimensional data well, this study is motivated to explore their effectiveness for Tm and  models.Therefore, based on three radiosonde stations in this study, the   model for Austria is developed using RF, ANN, and linear regression models.We then analyse the model performance against the radiosonde control station along with the Bevis and GPT3 models.Next, the hourly  for 39 GNSS stations located in the EPOSA (Echtzeit Positionierung Austria) network is modelled using RF and ANN algorithms.In the next step, the retrieved  is validated by the test period and four radiosonde stations in different period (August 2022 to April 2023).
This paper is organized as follows.Section 2 describes the data sources that were used for this study, including the GNSS  , radiosonde observations, and the ERA5 model.Section 3 discusses the techniques that are used for modelling   and .We investigate both the   and  data series in Section 4. Finally, we present our outlook and conclusions gained from this study in Section 5.

Data
In this section, we first describe the GNSS  time series employed in this study.Then, radiosonde observations are defined, which are used for both   modelling and  validation.In the end, the pre-processing step is described, which removes outliers in the data.Detailed information on all the datasets used in this study is provided in Tables 1 and 2.

GNSS ZTDs
A total of 39 permanent GNSS stations of the EPOSA network are considered in this study for the period between June 2018 and July 2022.GNSS sites in this network range in height from 172 m to 2215 m, extending from eastern to western Austria.The hourly GNSS  data are used in this study over two different time periods: training (January 2018-June 2021) and testing (July 2021-July 2022).To estimate  from GNSS measurements, we have used the Bernese software version 5.2 employing double-difference mode (please see [11,53] for more details).In Figure 1, the distribution of the various multi-GNSS (GPS + GLONASS) stations is visualized along with the location of the four control radiosonde stations (RS11035, RS11120, RS11010, and RS11240).

Radiosonde Observations
Radiosondes are highly accurate standard sensors in meteorology that measure a variety of meteorological parameters, including temperature (  = ±0.5°),pressure (  = [±1, ±2] ℎ), relative humidity (  = ±5 %), and geopotential height [54][55][56].Figure 2 depicts the distribution of radiosonde stations on the Austrian territory used in this study.To train and test the regional   model, we take into account the period of January 2010-June 2022.In addition, it should be highlighted that radiosonde profiles associated with the radiosonde stations have been downloaded free of charge via the website of the University of Wyoming.It is important to note that the radiosonde stations shown here correspond to the radiosonde stations in Figure 1.However, three of them are used for training while one is used to test the   model.Table 2 shows the location of the radiosonde stations in the case study, as well as the number of launches at different times (00, 06, 12, and 18) over the experimental period.

Data Pre-Processing
The presence of outliers in the measurements is inevitable and can negatively affect the performance and accuracy of a model if they are not addressed appropriately [30,[57][58][59].In data pre-processing, the Interquartile Range () method is widely used to detect outliers.It is based on the idea of quartiles, namely the first quartile ( 1 ), the second quartile ( 2 ), and the third quartile ( 3 ).Based on this, the  is defined as the difference between  1 and  3 , as follows: In the next step, using Equations ( 2) and ( 3), the lower and upper bounds for the data,  and , are determined as follows: where  is usually equal to 1.5.All data outside these bands are considered outliers and removed from the dataset.As a result, model performance will be enhanced.

Methodology
First, we describe Tm modelling using different techniques and data sources.Next, we discuss modelling with GNSS and radiosonde data sources.The main characteristics and features we employ in the RF and ANN models for   and  parameters are also described in those sections since we will apply these methods to model the   and  parameters later in this study.The main strategy of this research is demonstrated in Figure 3.

Tm Modelling
Here, Bevis and GPT3, two of the most commonly used models, are described with their main properties and formulas.Next, we present four different regional models based on regression (linear and polynomial) and machine learning for Austria.
In general, by considering the meteorological data of a radiosonde profile in different layers,   (K) can be calculated as follows [42,60]: here,  (hpa) and  (K) are water vapor pressure and temperature derived from the radiosonde profile, respectively; ∆ℎ (m) and  refer to the height difference between the adjacent data points, and the total number of layers, respectively.

Empirical Models
As part of this study, we examine and compare the accuracy of two widely-used   models, namely the Bevis and GPT3 models.The primary features of the Bevis and GPT3 models are outlined in Table 3.

• Bevis model
The Bevis model is an empirical model that estimates   using the linear relationship between   and   as follows [1]: To determine this model, Bevis et al. [1] utilized data from 8718 radiosonde stations in the United States over the years 1990-1991.In spite of the fact that this model has been widely used for many years in GNSS meteorology, it has some limitations, such as its simplicity as well as time and location independence.

•
GPT3 model GPT3 is the latest version of the empirical model GPT provided on a global grid of 5°× 5° and 1°× 1° [61,62].Various meteorological parameters can be obtained through this model, including pressure, temperature, and weighted mean temperature.Due to its high accuracy and simplicity [61], GPT3 is widely used in geodetic and meteorological fields.In Equation ( 6), the formula for estimating   for a given  (day of the year) is provided [45,61]: 0 represents the mean   value, while ( 1 ,  1 ) and ( 2 ,  2 ) represent the annual and semi-annual variations.These coefficients were estimated by means of an least-square adjustment using ERA-Interim data from ECMWF covering the years 2001-2010 [62].The GPT3 model can then be used to interpolate all different meteorological parameters, including   , for a specific location and time.In spite of this, it is still an empirical model, and thus it may not accurately reflect the actual state of the atmosphere in different regions or weather conditions [34,41,56].

Developed 𝑍𝑍 𝑚𝑚 Models
To develop a more accurate and consistent model for Austria, we established four models.Two of these models use traditional regression techniques, linear and polynomial, and the other two models use machine learning techniques, RF and ANN.The main characteristics of these models are presented in Table 4.
For the calculation of the  and  coefficients, we used data obtained from three radiosonde stations (RS11035, RS11240, RS1120) located in different parts of Austria between January 2010 and June 2021 (see Figure 2).By deriving linear regression coefficients, Equation ( 7) can be fitted for the Austrian region as shown below: •

Polynomial Model
As follows, this model is defined by the geographic location of the desired point (, ℎ) plus   : Using the same data (RS11035, RS11240, RS1120) from January 2010 to June 2021, the unknown coefficients ( 0 ,  1 ,  2 ,  3 ) were derived using the least-square method.Introducing these coefficients, the regional model for Austria states as follows: • Machine Learning Models (RF and ANN) In the 1950s, artificial intelligence pioneer Arthur Samuel defined machine learning as a subset of artificial intelligence [63][64][65].According to Samuel's definition, it is the field of study that enables computers to acquire knowledge without explicit programming [63].Plenty of complex problems can be modelled with this method, such as image classification, object detection, crop monitoring, wind speed prediction, groundwater modelling, and climate modelling [66][67][68][69].A number of scholars have used machine learning methods to model   on global and regional scales (e.g., [56,70,71]).To enhance the quality of   for Austria, we also apply two well-known machine learning techniques, namely RF and ANN.
Accordingly, the following model is considered for both RF and ANN algorithms: here, () refers to the machine-learning based estimation model for   .In order to obtain the machine learning models, again we used measurements at three different radiosonde stations (RS11035, RS11240, RS1120) located across Austria between January 2010 and June 2021 (see Figure 2).
-RF Method Beriman proposed the RF method as a supervised machine learning algorithm in 2001 [72], which can be used in classification and regression problems (for e.g., [71,[73][74][75][76]).To perform this method, a set of decision trees is created using random subsets of data and features [72,77].By creating such an uncorrelated set of trees, the risk of bias and overfitting is reduced [72].Once multiple decision trees have been calculated, the outputs are then summed for regression and voted for classification to come up with a single result [72,78].A schematic diagram of how this method is used to estimate the final output is shown in Figure 4. Equation ( 10) is implemented with the Scikit-learn machine learning library [79], module randomforest.regressor.The main parameters of this method have been tuned using the grid search method (see [80,81] for more details) by considering five cross-folds.Table A1 in Appendix A summarizes the hyper-parameters of the RF technique model.

-ANN Method
The idea behind neural networks was that they would simulate the human brain in some way [82].An ANN model consists of interconnected artificial neurons that can determine a general relationship between inputs and outputs without having a prior knowledge of the variables [82,83].As shown in Figure 5, the ANN architecture consists of an input layer, one or more hidden layers, and an output layer, with several neurons connected between them and running simultaneously [84].A variety of tasks can be accomplished using this method, such as environmental remote sensing [85], prediction of tropospheric wet delay and rainfall [84,86], and classification of crop yield [87].There are several hyperparameters that are critical for ANN, including the number of hidden layers, the number of neurons in each layer, the activation function, and the dropout rate.One and two hidden layers have been examined in the present study in order to determine how they affect the model's accuracy.Moreover, the Bayesian optimization method from the Python keras-tuner toolkit [88] has been used to determine the optimal set of hyperparameters.The method consists of predicting the optimum hyperparameter sets for a complex computational function from a probabilistic model [89,90].By using surrogate models, like Gaussian processes, the Bayesian optimization technique approximates a true objective function by using a distribution over the objective function [89][90][91].With Bayesian optimization, the model's performance can be evaluated computationally inexpensively even when the search space is large, as it smartly selects which hyperparameters to examine next, thus reducing the number of tests [92,93].Additionally, the validation mean square error (val-MSE) loss has been used to assess the performance of the model on unseen data in order to prevent overfitting.To do this, 20 per cent of the training dataset was used to monitor val-MSE's behaviour during training.
In Appendix A Table A2, we detail the best parameter settings for an ANN with one hidden layer to model   .Furthermore, another hidden layer was added in order to compare the performance of the double hidden layer neural network with the single hidden layer neural network.Table A3, shows the best hyperparameter sets for the double hidden layer neural network.
In the numerical results section, we evaluate the results obtained with the different methods described here in order to find the best model of   for the Austrian region.This model will then determine  in this area in the next step.

𝑃𝑃𝑃𝑃𝑉𝑉 Modelling
We first define the  based on GNSS data and radiosonde measurements.Following that, four different regional models of  are detailed adapted for the Austrian region based on RF regressions and ANNs techniques.

𝑃𝑃𝑃𝑃𝑉𝑉 Derived from GNSS Data
In order to derive  from GNSS data, we have first estimated  using Bernese GNSS software in baseline mode [11].Tropospheric zenith delay parameters were estimated on top of the apriori tropospheric delay model as piecewise linear functions.Loose constraints of 5 m for the initial parameter and 1 m for subsequent parameters were applied.A summary of the primary inputs and configurations used to estimate  in this study can be found in Table 5.In the next step,  has been computed using the Saastamoinen model as follows [12]: where,  (km),  (rad), and   (hpa) are the ellipsoid height, latitude, and surface pressure of a GNSS station, respectively.Using  as the result of subtracting  from , the  can be derived by multiplication of  with the conversion factor (Π) as noted below: With whereby,   (1000 kg/m 3 ) and   (461.51J/K×kg) denote the density of liquid water and the specific gas constant of water vapour, respectively.Moreover, the atmospheric refraction constants  2 ′ and  3 have respective values of 22.97 K/hpa and 375,463 K 2 /hpa [94].The   used in Equation ( 14) has been calculated using the best model in Section 3.1, which is investigated in the numerical results.

PWV Derived from Radiosonde Data
As already stated in Section 2.2, radiosonde provides accurate and reliable atmospheric parameter measurements.This allows us to calculate  by using Equation (15) [25]: In this equation, ∆ (hpa) and  (m/s 2 ) represent the pressure difference between adjacent layers and gravity acceleration, respectively.  (kg) refers to specific humidity and can be obtained as below [54]: here   (hpa) is the water vapour pressure and can be computed as follows [95]: whereby  (K) is the temperature at each pressure level.It is important to note that, for the  machine learning method developed from s of the GNSS stations, Equation ( 15) serves as an external validation reference.

Developed 𝑃𝑃𝑃𝑃𝑉𝑉 Models
As for   , we have used RF and ANN techniques in order to determine an hourly  in conjunction with GNSS  for the Austrian region.This is mainly aimed at improving weather forecasting and climate monitoring by delivering timely and accurate .Equations (18,19) show the general structure of the machine learning method, which describes the relationship between input parameters and output : where  1 () and  2 () refer to the machine-learning based estimation model for timely  in presence of four and five input variables, respectively.From now on, we will refer to the model with five parameters as Scheme#1 and the model with six parameters as Scheme#2.Further, to develop the machine learning model, we used data from 39 GNSS stations situated across Austria between January 2018 and June 2021 (see Table 1 and Figure 1).

-RF method
To develop the RF models in Scheme#1 and Scheme#2, we relied on Scikit-learn machine learning library [79], module randomforest.regressor.The grid search method has also been used for finding the best set of hyperparameters.Using grid search with five cross-folds, the best hyperparameters for Scheme#1 have been determined.Table A4 in Appendix B shows the obtained result.Additionally, Table A5 details the best set of hyperparameters for Scheme#2 obtained through the grid search with five cross-folds.In the numerical result section, the results gained using the RF method for Scheme#1 and Scheme#2 are discussed.
-ANN method ANN parameters have been tuned using Bayesian optimization techniques; for more information, please refer to Section 3.1.2.Accordingly, we have implemented ANN with one hidden layer for Scheme#1 and Scheme#2.The hyperparameters for the ANN method in Scheme #1 are listed in Table A6 in Appendix B. Further, a set of best hyperparameters for estimating hourly  using the ANN method of Scheme#2 is presented in Table A7.
The double hidden layer model has also been implemented for Scheme#1 to assess the impact of the extra hidden layer on ANN performance.Due to the fact that Scheme#2 did not differ significantly from Scheme#1 in terms of numerical results (please see Section 4.2.1),we did not implement a double hidden layer for Scheme#2.Table A8 details the best hyperparameter sets for the two-layer hidden neural network in Scheme#1.
Using the various approaches discussed here, we assess the numerical outcomes in order to determine the strengths and limitations of each, leading to the selection of the  model that performs best for Austria.

Statistical Metrics
To measure model fitness, the -squared indicator is used, which indicates how well the model fits the data [51,96].Generally, the closer the -squared value is to 1, the better the model fitted and, consequently, the more credible the fit.A formula for determining this statistical metric is as follows [96]: In this formula,  stands for Sum of Squared Residuals, and  stands for Sum of Squares Total.Moreover,    ,    , and   � represent the observed value ( or   ), the predicted value, and the average of the observed value.According to the Equation ( 20), -squared reflects the proportion of the variance explained by the model to its total variance.
Additionally, we have employed a set of statistical indicators to measure the model's performance quantitatively.The following indicators are included in the study: Mean Absolute Error (), Root Mean Square Error (), and Pearson correlation coefficient () [51,56,97,98]:

Numerical Results
The performance of both,   and  models is examined in this section.The developed model for   is evaluated internally using the radiosonde test data and externally using data from a check station (RS11010).A  model is also validated internally using GNSS test data and externally with radiosonde stations (RS11035, RS11120, RS11120, and RS11240).

Accuracy Evaluation of 𝑍𝑍 𝑚𝑚
Different statistical metrics like  and  are used here to evaluate the performance of the   models described in Section 3.1.2.We start by assessing the model's reliability and precision on the test dataset.In the next step, we will validate the efficiency of different models using RS11010 data.Based on the results, we will select the best   model to be applied to  calculation.

Internal Model Testing
Internal validation is a crucial step in the machine learning procedure in order to validate the performance and robustness of the trained model.The goal of this step is to demonstrate how well the developed model can generalize to unseen data and meet users' expectations regarding accuracy and reliability.Here, we have used the test dataset, which was separated from the whole dataset (see Table 1), covering the period July 2021-July 2022.A summary of the statistical indicators for the different   models is shown in Table 6.According to the reported results, the GPT3 model ( of 3.58 K and -squared of 50%) delivers slightly less accurate results in comparison with other models.Furthermore, the RF model provides better results with a  of 2.38 K, which is up to 30% better than for other methods but still comparable to ANN.A further analysis shows no significant difference between ANNs in both single-hidden ( of 2.45 K, -squared of 77%) and double-hidden layers ( of 2.46 K, -squared of 77%).This demonstrates that the   problem may not require more than one hidden layer because there is a nearly linear relationship between the defined features and the input data have relatively simple patterns.According to the training speed of the RF model, it would be more practical to use the RF model since it can achieve almost the same results as ANN within a shorter period of time.The purpose of external model testing is to verify that the model can adapt to changes in conditions and maintain its quality and reliability over time.In order to accomplish this, data from RS11010 have been used during August 2022-April 2023 to evaluate the performance of the discussed   models.
In Figure 6, the results of the empirical (Bevis and GPT3) and regression models (Polynomial and Linear) are shown.According to this figure, polynomial regression outperformed other models with an  of 2.70 K and a  of 3.25 K.Moreover, GPT3 has the weakest performance among the other models, degraded by about 17% in terms of .Figure 7 illustrates the developed models using RF and ANN methods.It can be seen from this figure that the RF model shows a slight improvement in  compared to the two other models, namely ANN with a single hidden layer (1 HL) and ANN with double hidden layers (2 HL).Moreover, the ANN model with one and two hidden layers produces almost identical results.Therefore, the results confirm once again that using two hidden layers for modelling   in this study has no significant advantage over one hidden layer.Based on the performance of different strategies to model Tm, RF appears to be the superior model in this case study with  improvements of approximately 2-19%.As a result, we employ the RF model in the next step to model  in Austria.

Accuracy Evaluation of 𝑃𝑃𝑃𝑃𝑉𝑉
We first assess the reliability of the model using the test dataset in this section.Afterwards, four radiosonde stations in Austria were analyzed to determine the accuracy of the model for predicting hourly .

Internal Model Testing
The GNSS  dataset covering July 2021-2022 was used to ensure the robustness of the  estimation using RF and ANN methods.Based on , , , and -squared (see Section 3.3 for more details), we have calculated the differences between  derived directly from GNSS and  estimated from RF and ANN models.As shown in Table 7, there is no significant difference between Scheme#1 and Scheme#2, which means adding longitude as an extra feature could only affect the model accuracy (, and ) by less than 2%.This indicates that the significant variations in both water vapour and pressure, and accordingly troposphere, tend to be more dependent on latitude and height than longitude.Furthermore, the double hidden layer was relatively ineffective in improving ANN model performance, and even one hidden layer was better than the double hidden layer.Accordingly, Scheme#1, using the ANN model with a single hidden layer, outperformed other models with  improvements by 1-15%.For the purpose of comparing the effectiveness of the RF and ANN models when applied to new, unseen data, we used data from four radiosonde stations (see Figure 1) from August 2022 to April 2023.This allows us to determine the accuracy and reliability of both models in real-world applications.Figure 8 illustrates the relationship between the predicted  using the RF method and radiosonde for RS11035 for Scheme#1 and Scheme#2.It can be seen from this figure that both methods have nearly the same overall pattern of points, although Scheme#2 has a slightly lower  and  than Scheme#1.However, in the case of RS11120 (see Appendix C, Table A9), there were some improvements in Scheme#1 over Scheme#2, which may be due to the fact that it is located in a more mountainous area than other radiosonde stations which might play a significant role in influencing  values.Nevertheless, similar to RS11035, two other radiosonde stations (RS11010 and RS11240) confirmed that Scheme#2 performs generally slightly better than Scheme#1.To acquire a better understanding of the models' performance, the , , and  of all radiosonde stations included in the external validation process were averaged.Based on the presented results in Table 8, Scheme#1 using the ANN with a single hidden layer outperformed the other models by almost 5% in  and 1-5% in .This result suggests that the developed ANN model (Scheme#1) can predict the  amounts in Austria with an average  of less than 2.5 mm.Moreover, as there was no significant difference between the RF and ANN models, the RF model may also be a good replacement for the ANN model in this case study.

Conclusions
In this study, a non-meteorological model was developed for estimating  in the Austrian region using a neural network technique.In order to accomplish this, first, we developed a regional model based on the location of the intended station, the surface temperature, the time of interest, and the day of the year.For the best   model, we compared the empirical (GPT3 and Bevis) models to the regression (linear and polynomial) as well as RF and ANN models.Based on the data from three radiosonde stations in different parts of Austria (RS11035, RS11120, and RS11240) from January 2010 to June 2021, regression and RF/ANN models were developed.A comparison of the RF model and the other models at the location of RS10101 revealed that the RF model provided superior results with  improvements of approximately 2-19% in the prediction phase from August 2022 to April 2023.Following the consideration of RF   , three models were constructed, namely RF, ANN with a single hidden layer, and ANN with a double hidden layer, in two separate schemes: Scheme#1 (, , ℎ, , ), and Scheme#2 (, , , ℎ, , ).These models were developed using  data from 39 GNSS stations covering the entire Austria territory between January 2018 and June 2021.A two-stage validation process was then performed on the developed models: an internal validation and an external validation.As a result, the performance of ANNs with one or two hidden layers was not significantly different.Additionally, longitude did not appear to be an important factor in improving , as tropospheric variation is influenced by latitude and height.According to the statistical results, it is possible to predict the hourly  in the Austria region without the need for sensed meteorological parameters with an average  of less than 2.5 mm.feedback and constructive suggestions.Moreover, we appreciate the guidance and support provided by the editors throughout the review process.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The best set of hyperparameters for RF and ANN techniques for modelling

Figure 1 .
Figure 1.Distribution of GNSS stations and four radiosonde control stations in Austria.

Figure 2 .
Figure 2. Distribution of three radiosonde stations to train the   model and one radiosonde control station to test the   model in Austria.

Figure 3 .
Figure 3. Main procedure for estimating   and  in this research.

Figure 4 .
Figure 4.The flowchart of the RF method in regression mode.

Figure 5 .
Figure 5. Schematic diagram of ANN architecture for regression.

Figure 7 .
Figure 7. Scatter plots of observed   against predicted   using different models: (a) RF, (b) ANN with one hidden layer, and (c) ANN with two hidden layer.

Figure 9
Figure9depicts the distribution of predicted  using ANN models compared to RS11035-derived .This figure shows that the estimated  using Scheme#1 with one hidden layer overcomes two other models with a  and  of 2.15 mm and 2.74 mm, respectively.Other radiosonde stations also support the same inference, as shown in Appendix C.

Figure 9 .
Figure 9. Scatter plots of RS11035  against predicted  using ANN model with: (a) single hidden layer for Scheme#1, (b) double hidden layer for Scheme#1, and (c) single hidden layer for Scheme#2.

Table 1 .
Information of the used data to model   and  parameters.

Table 2 .
Information on the used radiosonde stations in this study.

Table 3 .
Characteristics of the two selected models, Bevis and GPT3.

Table 4 .
Characteristics of the developed models based on traditional regression and machine learning techniques.

Table 5 .
Bernese GNSS software processing settings in this study.

Table 6 .
Statistical results of , , , and -squared of different   models for test dataset over the period of July 2021-July 2022.

Table 7 .
Statistical results of , , , and -squared of different  models for test dataset over the period of July 2021-July 2022.

Table A1 .
Hyperparameters tuning process of the final RF regressor for   modelling.

Table A2 .
Final tuned parameters for the single hidden layer ANN for   modelling.

Table A3 .
Final tuned parameters for the double hidden layer ANN for   modelling.The best set of hyperparameters for RF and ANN techniques for modelling

Table A4 .
Hyperparameters tuning process of the final RF regressor for  modelling in Scheme#1.areused when building trees.If False, the whole dataset is used to build each tree[False, True]

Table A5 .
Hyperparameters tuning process of the final RF regressor for  modelling in Scheme#2.areused when building trees.If False, the whole dataset is used to build each tree[False, True]

Table A6 .
Final tuned parameters for the single hidden layer ANN for  modelling in Scheme#1.

Table A7 .
Final tuned parameters for the single hidden layer ANN for  modelling in Scheme#2.

Table A8 .
Final tuned parameters for the double hidden layer ANN for  modelling in Scheme#1.Statistical results for different  models at RS11120, RS11010, and RS11240

Table A9 .
Statistical results of , , and  of different  models for RS11120 over the period of August 2022-April 2023

Table A10 .
Statistical results of , , and  of different  models for RS11010 over the period of August 2022-April 2023

Table A11 .
Statistical results of , , and  of different  models for RS11240 over the period of August 2022-April 2023