Regression Models for Soil Water Storage Estimation Using the ESA CCI Satellite Soil Moisture Product: A Case Study in Northeast Portugal

: The European Space Agency Climate Change Initiative Soil Moisture (ESA CCI SM) product provides soil moisture estimates from radar satellite data with a daily temporal resolution. Despite validation exercises with ground data that have been performed since the product’s launch, SM has not yet been consistently related to soil water storage, which is a key step for its application for prediction purposes. This study aimed to analyse the relationship between soil water storage ( S ), which was obtained from soil water balance computations with ground meteorological data, and soil moisture, which was obtained from radar data, as affected by soil water storage capacity (Smax). As a case study, a 14-year monthly series of soil water storage, produced via soil water balance computations using ground meteorological data from northeast Portugal and Smax from 25 mm to 150 mm, were matched with the corresponding monthly averaged SM product. Linear (I) and logistic (II) regression models relating S with SM were compared. Model performance ( r 2 in the 0.8–0.9 range) varied non-monotonically with Smax, with it being the highest at an Smax of 50 mm. The logistic model (II) performed better than the linear model (I) in the lower range of Smax. Improvements in model performance obtained with segregation of the data series in two subsets, representing soil water recharge and depletion phases throughout the year, outlined the hysteresis in the relationship between S and SM.


Introduction
Hydrological processes occurring on the Earth's surface and interacting with the lower atmosphere are key to understanding atmospheric dynamics and climate change and their effects on the global water cycle. In this context, soil moisture is one of the essential climate variables (ECVs) [1,2] playing an important role in the climate dynamics [3] and change [4,5], the water cycle and the soil-atmosphere energy exchange [6]. In fact, it influences hydrological and agricultural processes, runoff generation [7], drought development and many other atmospheric processes, such as actual evapotranspiration [8] and their spatial and temporal variability [9][10][11][12].
Climate and vegetation feedbacks are linked to soil moisture, meaning that soil moisture changes have a direct impact on natural ecosystems, forests and farmland productivity [2]. Advances in soil water dynamics research related to atmospheric drivers contribute towards improved water management in agriculture, increased irrigation efficiency and enhanced food security, especially in arid and semiarid areas [13][14][15], as well as providing Table 1. Summary of relevant literature using satellite-derived soil moisture products.

Studies Using: Data Assimilation in Land Surface Models Comparison with Ground Soil Moisture Data Soil Water Storage as Output Other Outputs
European Space Agency Climate Change Initiative Soil Moisture (ESA CCI SM) data - [26,29,[45][46][47][48][49][50]62] [ 25,36,37,[40][41][42][43][44]47,51] Other soil moisture remote sensing data [60,61] [ 24,31,33,34,63] [21, [30][31][32]40,41,66] In fact, water balance computations require the definition of a maximum soil water storage, or soil water storage capacity, whose effect on the relationship with SM is hardly reported in the literature [62,67]. This is an important and not yet fully explored issue that affects the use of satellite products measuring SM at smaller spatial scales, as the soil water storage capacity may present a wide range of spatial variation, depending on soil and crop characteristics [68]. On the other hand, a linear shape is commonly assumed in the relationships between the datasets issued from the different sources (satellite or ground), and regarding the different variables under consideration (soil moisture or soil water storage). Indeed, linearity is reported for the relationship between SM and the soil dielectric constant in most of the soil moisture range that is usually observed [69], as well as for that between total water storage and GRACE data [65]. No explicit support, other than statistical, is reported to justify the linearity assumption. Apart from the conventional approach of deriving statistical (data-driven) or functional (model-driven) relationships between datasets (matching remote data from different sources or matching remote data with ground sources), machine learning offers an alternative approach to generating largescale soil moisture estimates by also using radar products, with promising results [66]. Considering such a research background, this paper presents a case study in NE Portugal that was specifically aimed at analysing the statistical relationship between soil water storage, which was obtained from soil water balance computations with ground meteorological data, and soil moisture, as given by the ESA CCI SM product. As part of the analysis, the study also aimed at assessing the effect of soil water storage capacity in this relationship. In its broader purpose, the research supporting this study is intended to contribute towards enlarging the set of application and ground validation exercises of the ESA CCI SM product with ground data from Portugal. It is worth stressing that this exercise is performed in an especially important region, namely, the Mediterranean, where soil water shortage, and its prospected aggravation, is a critical water management challenge. Since it is based on long-term data series, the research also intended to respond to the growing demand for reliable tools for predicting soil water storage from satellite-borne data while accounting for actual soil crop variability.

Study Area
The study was carried out in the Bragança area, NE Portugal ( Figure 1). According to the Köppen-Geiger classification, the climate is temperate with dry and warm summers (Csb, Mediterranean type) [70,71], with the long-term annual averages of temperature and precipitation at Bragança being 12.3 • C and 758 mm, respectively [72]. Bragança is set in a mountain region with a roiling topography and strong regional climatic contrasts, ranging from humid and wet sub-humid at higher elevations (aridity index, i.e., the annual average precipitation to potential evapotranspiration ratio, AI > 0. 65) to semiarid in the topographic depressions (AI < 0.5) [73]. The regional surface area facing severe susceptibility to desertification and drought has significantly increased in the last few decades [74,75]. Threats to natural resources, such as soil and water, prevail all over the region, namely, soil degradation processes associated with high potential erosion risk, which is boosted by frequent forest fires [76][77][78]. In fact, almost 73% of the Bragança District area presents a semiarid (19.7%, 0.2 < AI < 0.5) or dry sub-humid (53%, 0.5 < AI < 0.65) climate [74], with high and moderate susceptibility to desertification, respectively, and 38% of the soils are degraded [71,75], which is dominantly concentrated in desertification susceptibility areas. The soil parent material is predominantly schist, although granite, mafic and ultramafic rocks are also found in the region [71]. The soils are mostly Leptosols (75.3%) [79], shallow and gravelly, dominantly medium-texture, acidic, with moderate to low organic  The soil parent material is predominantly schist, although granite, mafic and ultramafic rocks are also found in the region [71]. The soils are mostly Leptosols (75.3%) [79], shallow and gravelly, dominantly medium-texture, acidic, with moderate to low organic matter content [80]. The land suitability for agriculture is dominantly marginal and limited by a low effective soil depth, high soil stoniness, steep slopes and severe soil water shortage, with the deficit season lasting for up to 8 months in the warmer, drier and lower elevation areas [71,75].

Soil Water Balance Estimation with Ground Data
Soil water balance calculations were performed with ground data recorded at Quinta de Santa Apolónia weather station (QSA), which is managed by and located in the campus of the Polytechnic Institute of Bragança (41 • 47 48 N, 6 • 45 57 W, 681 m elevation) ( Figure 1). The Thornthwaite method [81] was adopted to estimate the reference evapotranspiration (ETo) using the QSA weather station monthly mean air temperature. Although the United Nations Food and Agriculture Organization (FAO) recommends the Penman-Monteith method [82], the Thornthwaite method requires fewer input data (air temperature and location latitude) and applies a more straightforward calculation routine; therefore, it requires less robust datasets, which is beneficial when dealing with missing values is a common challenge. In fact, [54] used the Surface Energy Balance Algorithm for Land (SE-BAL) to estimate ETo with the Penman-Monteith method from MODIS and TRMM satellite data, where they pointed out missing input data as a drawback for model validation. Moreover, in view of replication of the present study's methodological approach and for results upscaling purposes, Thornthwaite ETo estimation requires the most commonly available large-scale datasets, which are too restrictive for the Penman-Monteith application.
Soil water balance computations followed the procedures of the Thornthwaite-Mather method [83,84]. Requiring monthly ETo and P (precipitation) values, the method outputs estimates of monthly soil water storage (S) and its changes (∆S) throughout the year, together with real evapotranspiration (ETr), soil water deficit (WD) and superavit (WS). A maximum soil water storage must be defined a priori (Smax, millimetre equivalent depth), according to soil and crop characteristics, namely, available soil water capacity (AWC) and rooting depth (z) [83]. The reference Smax for climatic classification purposes is 100 mm [81].
The Thornthwaite-Mather method was applied using the QSA weather station monthly data for each year of the study period (January to December). The year totals of each input (P, ETo) and output variable (ETr, WD, WS) were also calculated. Deficit, superavit and recharge periods were identified, with the latter occurring when P − ETo > 0 and S < Smax, and the first two corresponding to the dry and wet seasons, respectively. For comparison purposes, the method was equally applied with the mean monthly P and ETo values of the study period and the climate normal . Procedures were repeated for Smax 150 mm (representing high soil water storage conditions), Smax 100 mm (representing reference Smax), and Smax 75, 50 and 25 mm (representing moderate-low, low and very low soil water storage, respectively), thus covering a wide range of soils (with different available water capacities) and/or crops (with different rooting depths). It should be added that soils of the study area are dominantly shallow, poor in organic matter and gravelly, and therefore better represented by the Smax lower range, while the upper range represents the scarce area assigned to deeper soils, such as Cambisols and Fluvisols [71,80].

ESA CCI Satellite Product
SM is an ESA CCI product (V 04.2) that utilises four active radar products (AMI WS, ERS-1/2 SCAT, MetOp-A+B ASCAT) and seven passive radar products (SMMR, SSM/I, TRMM, Windsat, AMSR-E, AMSR2, SMOS) microwave sensors; then, data is processed by the Microwave Remote Sensing Group at TU Wien [22,85]. The product is freely available in the Copernicus database (https://cophub.copernicus.eu/dhus/#/home), with a spatial resolution of 0.25 • , a time resolution of 1 day and the data series started in 1978. The GIS SNAP (Sentinel Application Platform), version 6.0, a freely available tool found at http: //step.esa.int/main/download/, was used for the satellite imagery processing.
Space and time attributes of the SM data were considered prior to comparisons with the ground data, which meant (i) appraising the representation in the SM pixel area of the ground climatic conditions, (ii) setting the study period for which the ground and satellite data were both available and (iii) representing the available daily SM data with a monthly time scale.
As expected in mountainous areas, the elevation is a major factor that determines the regional climate distribution, where for NE Portugal, the mean annual precipitation is positively correlated with elevation (r = 0.513), while the mean annual temperature and elevation are negatively correlated (r = −0.738), with the two climate elements varying inversely in the region (r = −0.559) [86]. Considering this, a 30 m resolution Digital Elevation Model (Advanced Land Observing Satellite (ALOS) World 3Dhttps://www.eorc.jaxa.jp/ALOS/en/aw3d/index_e.htm) was used to derive the elevation areal distribution in pixels, as depicted in Figure 1. The QSA elevation fell close to the pixel's median (P50), with a 42 m difference for an effective altitudinal range of 327 m (P95 − P5, with P referring to percentile). Therefore, it was assumed that the QSA ground climatic conditions provided a good representation of those prevailing in the SM pixel area.
The selected SM data series extended for 14 years (2003-2016), which was the longest continuous timespan with simultaneously available ground and satellite-borne data. In fact, earlier than 2003, the ESA CCI SM product imagery existing for the area contains a large set of missing data in the daily SM series.
For obtaining a monthly SM data series, which is required for a comparison with the monthly water balance that is computed with the ground data, three randomly selected years were analysed in detail. For these, the monthly SM averages were computed (i) with daily data before and after the accumulated period of precipitation and (ii) with data of three fixed days: the 1st of the month, the 15th of the month and the 1st of next month. It should be noted that missing data commonly occurred in these years, meaning that the monthly averages were computed with the available data. Non-significant differences were found between the two sets of monthly averages, which showed very similar individual values and were strongly correlated (r = 0.896). These results encouraged the application of the procedure for the entire study period.
Furthermore, the SM pattern of evolution throughout the year and its changes during the study period were analysed. For this purpose, the SM series was divided into two subsets, each one assembling the SM rise and decline periods, corresponding to the soil water recharge or wetting periods and to the soil water depletion or drying periods, respectively. Soil water depletion occurs from the highest SM value in the wet season to the lowest at the end of the dry season, whereas soil recharge occurs in the remaining months of the year. For the entire study period, the drying data subset was larger than the wetting subset.

Regression Models
Monthly SM and water balance data for the study period (2003-2016) were summed up to give a total of 168 values for each variable. Regression analysis was performed between the satellite-borne and the ground data series, taking S (mm) as the variable that was dependent on SM (m 3 m −3 , converted to a percentage, %) due to the prediction purposes being explored in this research.
Two regression models were tested (S* stands for the S model estimate): Model II (Logistic) : Water 2021, 13, 37 7 of 26 In model I, the dependent variable was bounded by the domain limits that were defined in the water balance computation (0 < S < Smax), meaning that all S estimates (S*) falling outside this range were set equal to the nearest domain limit (0 or Smax). The plots of SM against S suggested a non-linear relationship between these variables, and it was hypothesised that a logistic function (model II) would represent the relationship better than the linear function (model I). Furthermore, model II is asymptotically bounded to 0 and Smax, therefore providing a more realistic approach to the data range extremes than model I. Among other possible choices [87], the adopted logistic function [88] allowed for simple linearisation such that linear regression analysis was performed with its linearised form.
Regression analysis was performed with data from the entire study period (14 years), which involved testing the performance of the two models regarding statistically predicting S with SM data. As indicated earlier, water balance computations were repeated for a range of Smax values (25,50,75, 100 and 150 mm), which yielded the respective S series. Regression models were applied for the five soil water storage conditions in order to test the effect of Smax on the models' performances. Regression analysis was also performed separately with the wetting and the drying data subsets to appraise differences between the subsets in the SM-S relationship. In this case, only model II was applied with the S data issued from water balance computations for the five soil water storage conditions. Finally, regression analysis was performed independently for each year of the study period using both models I and II in order to assess the interannual variability in the models' performances and search for factors explaining this variability. This was done with the S data for all Smax conditions. The models' performances were evaluated using performance indicators, namely, the determination coefficient (r 2 ), efficiency coefficient (CE), root mean square error (RMSE, mm) and standard error of the estimate (SEE, mm), which were all applied to the observed and estimated S data series [89].

Soil Water Balance and Soil Moisture during the Study Period
During 14 years, the study period's annual P and ETo were, on average, 751.3 and 711.4 mm, respectively. The annual ETo ranged from 659.8 mm, in 2017, to 751.5 mm, in 2015. A far higher data dispersion was found for P ( Figure 2), with a coefficient of variation (standard deviation/mean) of 0.45, against 0.03 for ETo. The driest of the 14 years was 2012 (P = 359.7 mm), closely followed by 2005 (382.4 mm), which is a figure below the eighth percentile. The largest P record was 1584.6 mm in 2010, an exceptionally high value for Bragança, though not considered a statistical outlier of the annual P series.
The annual values of the soil water balance (SWB) output variables for the study period are also depicted in Figure 2. The average annual ETr, WD and WS were quite similar during the study period (383.7, 327.7 and 367.6 mm, respectively), though with very different data dispersions, with the highest being for WS.
The average monthly SWB in Bragança, which was computed using the P and ETo monthly averages for the study period, reflected the typical Mediterranean pattern prevailing in the study area ( Figure 3a). While the periods of soil water deficit (WD) and superavit (WS) were clearly identified, the soil water recharge was hardly noticeable in the September-October months, as the first autumn rainfalls filled the soil storage capacity up to its maximum, generating a superavit from October to April. The water deficit season lasted from May to September, with a peak deficit in July (Figure 3a).
Average annual values of the SWB variables for the study period showed small differences in comparison with the long-term averages, which were positive for WD and WS, negative for ETr and negligibly negative for ETo and P (Figure 3b). The long-term SWB was computed using the last published climatological normal from 1971-2000, namely, with the thirty-year monthly average temperature and precipitation [72]. Their monthly distribution throughout the year was similar in both calculation periods (Figure 3a), though with a higher P in October in the study period than in the 1971-2000 one. Despite the aforemen-tioned differences, the study period was assumed to provide an adequate representation of the normal climatic conditions prevailing in the case study area. 711.4 mm, respectively. The annual ETo ranged from 659.8 mm, in 2017, to 751.5 mm, in 2015. A far higher data dispersion was found for P ( Figure 2), with a coefficient of variation (standard deviation/mean) of 0.45, against 0.03 for ETo. The driest of the 14 years was 2012 (P = 359.7 mm), closely followed by 2005 (382.4 mm), which is a figure below the eighth percentile. The largest P record was 1584.6 mm in 2010, an exceptionally high value for Bragança, though not considered a statistical outlier of the annual P series.
The annual values of the soil water balance (SWB) output variables for the study period are also depicted in Figure 2. The average annual ETr, WD and WS were quite similar during the study period (383.7, 327.7 and 367.6 mm, respectively), though with very different data dispersions, with the highest being for WS. The average monthly SWB in Bragança, which was computed using the P and ETo monthly averages for the study period, reflected the typical Mediterranean pattern prevailing in the study area ( Figure 3a). While the periods of soil water deficit (WD) and superavit (WS) were clearly identified, the soil water recharge was hardly noticeable in the September-October months, as the first autumn rainfalls filled the soil storage capacity up to its maximum, generating a superavit from October to April. The water deficit season lasted from May to September, with a peak deficit in July ( Figure 3a).  Average annual values of the SWB variables for the study period showed small differences in comparison with the long-term averages, which were positive for WD and WS, negative for ETr and negligibly negative for ETo and P (Figure 3b). The long-term SWB was computed using the last published climatological normal from 1971-2000, namely, with the thirty-year monthly average temperature and precipitation [72]. Their monthly distribution throughout the year was similar in both calculation periods (Figure 3a), though with a higher P in October in the study period than in the 1971-2000 one. Despite the aforementioned differences, the study period was assumed to provide an adequate representation of the normal climatic conditions prevailing in the case study area.
During the study period, the extent of the dry season (number of months) and annual WD (mm) were positively correlated (r = 0.725). A large range in the dry season extent occurred throughout the series, with the dry season lasting from 4 to 9 months with a median duration of 7 months (Figure 4). The annual WD was more dependent on the  During the study period, the extent of the dry season (number of months) and annual WD (mm) were positively correlated (r = 0.725). A large range in the dry season extent occurred throughout the series, with the dry season lasting from 4 to 9 months with a median duration of 7 months ( Figure 4). The annual WD was more dependent on the starting than on the ending month of the dry season. The longest and the shortest recorded wet seasons lasted for 8 and 3 months, respectively, and a very strong correlation was found between the WS and P annual values (r = 0.983).

Relationship between the Soil Water Storage and the Satellite-Borne Soil Moisture Data
The overlap for the study period of the monthly ground SWB data (soil water storage, S) and satellite-borne data (SM as a volume percentage) is shown in Figure 5. Close matching of the two series was found for this case study and the data exploration proceeded to regression analysis. Indeed, the wet and dry seasons, as well as the soil water recharge and depletion periods, were clearly identified. However, the SM kept varying in the wetter months, while S flattened at Smax (100 mm in the  The monthly SM ranged from 0.085 to 0.311 m 3 /m 3 , while the monthly averages computed for the study period were the lowest in August (SM = 0.110) and the highest in January (SM = 0.269). The lowest annual SM average was found in 2012 (SM = 0.171), the year with the lowest P in the series. The average annual SM was negatively correlated with the extent of the dry season (r = −0.843).

Relationship between the Soil Water Storage and the Satellite-Borne Soil Moisture Data
The overlap for the study period of the monthly ground SWB data (soil water storage, S) and satellite-borne data (SM as a volume percentage) is shown in Figure 5. Close matching of the two series was found for this case study and the data exploration proceeded to regression analysis. Indeed, the wet and dry seasons, as well as the soil water recharge and depletion periods, were clearly identified. However, the SM kept varying in the wetter months, while S flattened at Smax (100 mm in the  The monthly SM ranged from 0.085 to 0.311 m 3 /m 3 , while the monthly averages computed for the study period were the lowest in August (SM = 0.110) and the highest in January (SM = 0.269). The lowest annual SM average was found in 2012 (SM = 0.171), the year with the lowest P in the series. The average annual SM was negatively correlated with the extent of the dry season (r = −0.843).

Relationship between the Soil Water Storage and the Satellite-Borne Soil Moisture Data
The overlap for the study period of the monthly ground SWB data (soil water storage, S) and satellite-borne data (SM as a volume percentage) is shown in Figure 5. Close matching of the two series was found for this case study and the data exploration proceeded to regression analysis. Indeed, the wet and dry seasons, as well as the soil water recharge and depletion periods, were clearly identified. However, the SM kept varying in the wetter months, while S flattened at Smax (100 mm in the The model performance indicators (r 2 , CE, RMSE and SEE), together with the parameters of the regression equations relating soil water storage (S, based on the ground data) with the soil moisture (SM, satellite-borne data) are presented in Table 2. The tested models, described in the Materials and Methods section, were linear bounded (I) and logistic (II). These were applied for Smax ranging from very low (25 mm) to high (150 mm). Figure 6 depicts an example of models' fits to the data for Smax 50 mm. gure 5. Monthly series of soil moisture (SM), an ESA satellite product, and soil water storage (S), an outcome of the ornthwaite-Matter water balance that was computed for Bragança (Smax = 100 mm).
The model performance indicators (r 2 , CE, RMSE and SEE), together with the parameters of the regression equations relating soil water storage (S, based on the ground data) with the soil moisture (SM, satellite-borne data) are presented in Table 2. The tested models, described in the Materials and Methods section, were linear bounded (I) and logistic (II). These were applied for Smax ranging from very low (25 mm) to high (150 mm). Figure  6 depicts an example of models' fits to the data for Smax 50 mm.  The determination coefficients (r 2 ) ranged from 0.690 (model II, logistic, with Smax 150 mm) to 0.886 (model II, logistic, with Smax 50 mm), meaning a minimum correlation coefficient (r) of 0.831 when all models and Smax's were considered. For the range of Smax's tested, the maximum r 2 was reached with 50 mm for both models (Figure 7a). Model II showed a higher r 2 than model I for Smax 25 to 75 mm, and a lower r 2 for the higher Smax's (100 and 150 mm). A similar model performance pattern of change was found for CE, with values slightly lower than r 2 . RMSE ranged similarly in both models from 5 mm (Smax 25 mm) to 29 and 34 mm (Smax 150 mm) in models I and II, respectively. A similar increase with the Smax increase was also found for SEE, though with lower values and more pronounced differences between the models (4-24 mm in model I and 5-31 mm in model II in the Smax range 25-150 mm). When expressed as a proportion of Smax, SEE fell into a much narrower range in model I (15-16%) and model II (15-21%, with 18% as the median) (Figure 7b). The determination coefficients (r 2 ) ranged from 0.690 (model II, logistic, with Smax 150 mm) to 0.886 (model II, logistic, with Smax 50 mm), meaning a minimum correlation coefficient (r) of 0.831 when all models and Smax's were considered. For the range of Smax's tested, the maximum r 2 was reached with 50 mm for both models (Figure 7a). Model II showed a higher r 2 than model I for Smax 25 to 75 mm, and a lower r 2 for the higher Smax's (100 and 150 mm). A similar model performance pattern of change was found for CE, with values slightly lower than r 2 . RMSE ranged similarly in both models from 5 mm (Smax 25 mm) to 29 and 34 mm (Smax 150 mm) in models I and II, respectively. A similar increase with the Smax increase was also found for SEE, though with lower values and more pronounced differences between the models (4-24 mm in model I and 5-31 mm in model II in the Smax range 25-150 mm). When expressed as a proportion of Smax, SEE fell into a much narrower range in model I (15-16%) and model II (15-21%, with 18% as the median) (Figure 7b). In both models, the b regression parameters increased with the Smax increase, while the opposite occurred with the a parameters; these changes expressed the effect of the S scale in its relationship with SM. The best-fitting functions relating a and b with Smax were obtained for both models with the following results: Model I-Linear bounded b = 0.822 + 0.049Smax, Model II-Logistic b = 4.260Smax 0.417 , a = 157.290Smax 0.585 , r 2 = 0.994.
In all four cases, the correlation with Smax was statistically significant and, except for the model I a parameter (Equation (4)), the r 2 of the regression functions were very high. As expected, non-linear functions of Smax best fit the a and b parameters of the lo- In both models, the b regression parameters increased with the Smax increase, while the opposite occurred with the a parameters; these changes expressed the effect of the S scale in its relationship with SM. The best-fitting functions relating a and b with Smax were obtained for both models with the following results: Model I-Linear bounded Model II-Logistic In all four cases, the correlation with Smax was statistically significant and, except for the model I a parameter (Equation (4)), the r 2 of the regression functions were very high. As expected, non-linear functions of Smax best fit the a and b parameters of the logistic model. The above equations allow for estimates of the S vs. SM regression parameters for any Smax, therefore providing, for each model, a single relationship S vs. SM that simultaneously accounts for Smax in an S estimation, as follows: Model I-Linear bounded S * = − 19.628 − 0.325Smax + 0.822SM + 0.049Smax × SM, r 2 = 0.856, CE = 0.856.
As shown in Figure 6 for Smax 50 mm, the logistic model (II) provided a visibly more adequate representation of the relationship between S and SM compared with the linear model (I). However, the linear model overcame the logistic model's performance (measured using r 2 and CE) for larger Smax's, depicting a lower SEE for all Smax's ( Figure 7). In the models given above with two variables (SM and Smax; Equations (7) and (8)), the linear and the logistic models had similar performances (as measured using r 2 ) in the whole range of Smax, although the linear model showed a slightly higher CE compared to the logistic model. It should be noted that Equations (7) and (8) describe the prediction surfaces of S responding to the SM and Smax input data, as calibrated for the Bragança area, which performed better in the first case but was better represented in the second one.

Soil Water Depletion and Recharge Phases
As stressed earlier, in a large number of wet months, from late autumn to early spring, Smax was reached, meaning a constant S value, while SM was not constant in those months. On the other hand, the soil water recharge during early autumn was much faster than the depletion during spring. This evidence encouraged deeper analysis of the S vs. SM coupled data series, with the regression analysis for each year of the study period and a chronological plot of the relationship SM vs. S shown in Figure 8.
Taking Smax 100 mm as an example, the quality of the relationship varied sharply throughout the study period, with r 2 ranging from 0.532 and 0.960 for the linear model and from 0.414 to 0.945 in the logistic model. Keeping Smax 100 mm as an example, Figure  8 depicts the large differences between years of the study period in terms of the monthly chronological changes of SM and S, with the four years depicted being ranked from low to high r 2 of the relationship. A large part of the recharge (wetting) and depletion (drying) phases presented a distinct S for the similar SM and vice versa. This was increasingly more evident as r 2 decreased, with clockwise progress from January to December being visible in 2008, where r 2 was as low as 0.600. As such, the model performance for the whole study period was seemingly affected by the different shape of the S vs. SM relationship in the soil water recharge and depletion phases. Following this, the previous data analysis steps were performed separately for the two data subsets in which the original series was segregated: soil water recharge months and depletion months. Table 3 shows the results of the segregated model (III) applied to the drying curve, wetting curve and drying + wetting curve. For each series and maximum storage tested, the drying curves presented a higher r 2 than the wetting curves, with SEE varying in the opposite way. The drying curve was fitted to a larger data set, 91 out of 168 months, while the wetting curve was fitted to the remaining 77 months. For the drying curves, r 2 ranged from 0.810 (Smax 150 mm) to 0.875 (Smax 50 mm), and the same trend was observed for CE. For the wetting curves, the best r 2 was obtained at Smax 75 (0.781), while a higher CE was associated with Smax 50 mm. For both curves and the drying + wetting curves, RMSE and SEE showed the same patterns found in models I and II, with both increasing as Smax increased. Furthermore, the percentage of SEE in Smax ranged from 12 to 18% (drying curve), 13 to 18% (wetting curve) and 14 to 19% (drying + wetting curve). Compared to the logistic approach applied to all data (model II, Table 2), data segregation in the drying and wetting periods (model III, Table 3) produced better global predictability of S from SM, presenting a higher r 2 and a lower SEE for each Smax tested (Figure 9). Improvements in the model's performance were more pronounced in larger Smax (from 75 to 150 mm), where segregation of the data series was visibly more effective in increasing r 2 and lowering SEE. Smax 50 mm remained the best-performing soil water storage capacity (r 2 = 0.892 in model III against r 2 = 0.886 in model II); Figure 10 depicts model III fitted to S vs. SM data with Smax 50 mm.  Table 3 shows the results of the segregated model (III) applied to the drying curve, wetting curve and drying + wetting curve. For each series and maximum storage tested, the drying curves presented a higher r 2 than the wetting curves, with SEE varying in the opposite way. The drying curve was fitted to a larger data set, 91 out of 168 months, while the wetting curve was fitted to the remaining 77 months. For the drying curves, r 2 ranged from 0.810 (Smax 150 mm) to 0.875 (Smax 50 mm), and the same trend was observed for CE. For the wetting curves, the best r 2 was obtained at Smax 75 (0.781), while a higher CE was associated with Smax 50 mm. For both curves and the drying + wetting curves, RMSE and SEE showed the same patterns found in models I and II, with both increasing as Smax increased. Furthermore, the percentage of SEE in Smax ranged from 12 to 18% (drying curve), 13 to 18% (wetting curve) and 14 to 19% (drying + wetting curve). Compared to the logistic approach applied to all data (model II, Table 2), data segregation in the drying and wetting periods (model III, Table 3) produced better global predictability of S from SM, presenting a higher r 2 and a lower SEE for each Smax tested (Figure 9). Improvements in the model's performance were more pronounced in larger Smax (from 75 to 150 mm), where segregation of the data series was visibly more effective in increasing r 2 and lowering SEE. Smax 50 mm remained the best-performing soil water storage capacity (r 2 = 0.892    Following the approach applied to models I and II, Equations (7) and (8) allow for S estimations from SM and Smax, considering the a and b parameters shown in Table 2. The r 2 's of the non-linear relationships between the model regression parameters and Smax were all very high. In fact, for the drying curve a and b parameters, r 2 was 0.986 and 0.957, respectively; for the wetting curve a and b parameters, r 2 was 0.989 and 0.991, respectively.
Model III-Segregated drying and wetting curve The combined application of the two above equations (Equations (9) and (10)) with the respective data subsets produced the global response (drying + wetting) of model III, with r 2 = 0.887 and CE = 0.868. Compared to the results obtained using Equation (8), the segregated two-variable logistic model that was used to account for soil water recharge and depletion periods showed better performance for the whole range of Smax tested, reaching a higher r 2 and CE. Again, these results highlight the effect of segregating data series in the wetting and drying subsets in improving the model's performance and predictive capacity in the estimation of S from SM.

Outliers' Effect on the Model Performance
The results presented correspond to an unscreened series of SM and S data. Yet, from the 168 monthly values, a small number fell far from the general trend defining the relationship between S and SM. No plausible explanation for such deviations was devised and; therefore, a screening exercise for identifying outliers was carried out.
As a first approach, the standardised residuals (SRs) of the trend lines fitted to S vs. SM data (for each one of the Smax tested) were used to screen outliers [90]. This approach implies defining an SR threshold as a screening criterion, which is a user-dependent decision that is supported by experience, expert assistance or the literature [89]. Absolute values of the standardised residuals (ABS SRs) of 3.0, 2.5 and 2.0 were adopted in the exercise, although no guiding values exist for this case. For ABS SR > 3.0, 2 outliers were detected, for ABS SR > 2.5, 7 outliers were detected, and for ABS SR > 2.0, 16 outliers were detected. Following the approach applied to models I and II, Equations (7) and (8) allow for S estimations from SM and Smax, considering the a and b parameters shown in Table 2. The r 2 's of the non-linear relationships between the model regression parameters and Smax were all very high. In fact, for the drying curve a and b parameters, r 2 was 0.986 and 0.957, respectively; for the wetting curve a and b parameters, r 2 was 0.989 and 0.991, respectively.
Model III-Segregated drying and wetting curve The combined application of the two above equations (Equations (9) and (10)) with the respective data subsets produced the global response (drying + wetting) of model III, with r 2 = 0.887 and CE = 0.868. Compared to the results obtained using Equation (8), the segregated two-variable logistic model that was used to account for soil water recharge and depletion periods showed better performance for the whole range of Smax tested, reaching a higher r 2 and CE. Again, these results highlight the effect of segregating data series in the wetting and drying subsets in improving the model's performance and pre- In a second approach, the interquartile range rule (IQR) [91] was applied to the residuals of the trend line fitted to S vs. SM data (for each one of the Smax tested). Residuals falling out of the range defined by the following limits corresponded to outliers in the relationship between S and SM (Q 1 and Q 3 are the first and third quartiles of the residuals' series): Eight outliers were detected using these rules. As a criterion that is not user-dependent, it was preferred to the ABS SR threshold approach. The performance of the regression models was tested using the data series after the outliers identified using the IQR were discarded (n = 160 for each Smax; Table 4). Furthermore, the IQR provided guidance for defining the ABS SR threshold value when the first approach was applied (in this case, the ABS SR was approximately 2.5). Table 4. Regression parameters and performance indicators of models I, II and III, which were applied with different soil water storage capacities, Smax: outlier-removed data series obtained using the interquartile range rule (SM, %, independent variable; S, mm, dependent variable; n = 160). The models' performances visibly improved after the outliers were removed from the data series, as r 2 ranged from 0.846 to 0.901 in the linear model and from 0.886 to 0.943 in the logistic model (compare with Table 1). Comparing Tables 1 and 2, the same pattern of change was observed for CE, while a decrease was found for RMSE and SEE, even though the SEE decrease was too small to determine changes when expressed as a proportion of Smax. In contrast, the model parameters did not visibly change when comparing results coming from the original data series with those obtained with the outlier-removed data series. This means the outliers' removal produced an increase in the models' performances without a visible change in the equations' parameters.

Model
These findings are confirmed in Figure 11, where r 2 , SEE (in % of Smax) and the a and b parameters of the regression functions for the original (independent variable) are plotted against those obtained with the outlier-removed data series (dependent variable). The slopes of the trend lines fitted to plotted data, with a null intercept, ranged from 0.904 (SEE, %Smax) to 1.061 (r 2 ), while for the a and b regression parameters, the trend line slopes were 0.975 and 1.005, respectively. The outlier removal increased the global r 2 by 6%, while for the global SEE, a decrease of ≈10% was observed. For all the models and storage capacities analysed, the removal of outliers did not substantially change the regression equations' parameters (a decreased by 2.5% and b increased by 0.5%). Furthermore, a and b presented the strongest r 2 's, with 0.9997 (a) and 0.99998 (b).
Water 2021, 13, x FOR PEER REVIEW 19 of 28 Figure 11. Linear fit and determination coefficient of the model performance indicators (r 2 and SEE, %Smax) and the regression parameters (a and b) obtained with all models and soil water storage capacities using the original data series and the outlier-removed data series. This is an important practical finding as it shows the robustness of the SM vs. S relationship, which was not sensitive to the presence of outliers in the data set; therefore, saving screening procedures prior to the model application. These findings also open the possibility of reducing the dataset to a minimum number of data necessary for model calibration. Thus, the equations obtained using the original data can be taken as a reference for prediction purposes. On the other hand, the model performance indicators (r 2 and CE, obtained using the outlier-removed data series) can be taken as a reference since the withdrawn data represent uncommon and unexplainable records that are meaningless for defining the relationship between SM and S.

Model Comparison
The models describing the relationship between S and SM performed well, showing fairly high-quality standards (r 2 > 0.8 and CE > 0.8, all models and Smax accounted for). However, the high r 2 values that were obtained do not entirely describe the models' per- This is an important practical finding as it shows the robustness of the SM vs. S relationship, which was not sensitive to the presence of outliers in the data set; therefore, saving screening procedures prior to the model application. These findings also open the possibility of reducing the dataset to a minimum number of data necessary for model calibration. Thus, the equations obtained using the original data can be taken as a reference for prediction purposes. On the other hand, the model performance indicators (r 2 and CE, obtained using the outlier-removed data series) can be taken as a reference since the withdrawn data represent uncommon and unexplainable records that are meaningless for defining the relationship between SM and S.

Model Comparison
The models describing the relationship between S and SM performed well, showing fairly high-quality standards (r 2 > 0.8 and CE > 0.8, all models and Smax accounted for). However, the high r 2 values that were obtained do not entirely describe the models' performances as large SEE and RMSE values were also found. Rather than compromising the models' descriptive capabilities, these results should be outlined as a limitation for model prediction purposes. Even though outlier removal led to a decrease in SEE, as indicated above, the reduction was about 10% globally, which is an improvement compared to the original series but too low of a gain in terms of model predictive power.
Dorigo et al. [43] compared the ESA CCI SM combined product with in-situ-based soil moisture observations recorded by the validation ground stations network. The median of the correlation coefficients obtained in each station for the latest merging period (July 2012 to December 2014) was 0.89, varying from 0.96 to 0.79 in the interquartile range (third and first quartiles of the series of correlation coefficients obtained in each station). These r values correspond to r 2 's of 0.79, 0.92 and 0.63, respectively. The model performance results in the present study fell in the third to second quartile ranges of r 2 values obtained in the network of stations where ground soil moisture data validated the SM combined product data [43].
Hydrological and land surface models predicting soil moisture and total soil water storage (TWS) achieved lower performances [60,61]. Tian et al. [60] tested data assimilation effects on the water balance model response, which was performed with soil moisture derived from SMOS and total water storage derived from GRACE. Data assimilation from these sources improved the model surface SM estimates when compared with in situ SM measurements, jointly performing better than the single assimilation (r = 0.47-0.86). A similar approach using a land surface model resulted in improved model estimates, which were matched with in situ measurements, of the surface SM (r = 0.86) and TWS (r = 0.95) [61]. In both studies, no direct relationship between remote and ground SM was provided.
The logistic model generally performed better than the linear model. However, the differences in r 2 for the logistic model compared to the linear model were small and, for the larger Smax, the latter performed better than the former. As such, in the case of constraints in the calculated means, the linear function can be adopted to convert SM into S with a minimal loss of predictive power. Apart from its practical simplicity, the linear model showed equally practical constraints, as it was necessarily bounded by upper (S = Smax) and lower limits (S = 0), meaning that the linear function only applies to the 0 to Smax range. This actually means that for SM < −a/b and SM > (Smax − a)/b, S is not predicted by the linear function but is kept constant and equal to 0 or Smax, respectively (a (intercept) and b (slope) are the linear regression parameters). Therefore, the linear model coupled three functions within the range of SM, with the thresholds for shifting from one to the other being variable with Smax.
In contrast, with the logistic model, a single function enabled describing the S vs. SM relationship for the whole range of SM, which is a practical advantage for calculation routines. Furthermore, the approach to S = 0 and S = Smax is asymptotical, meaning that, although the residuals of the regression line are higher in the logistic model compared to the linear model for the range of SM in which S = Smax and S = 0, the opposite occurs for the data not strictly matching this condition. In fact, high SM values not matching S = Smax or S = Smax values not matching high SM values generate lower the residuals of the regression line in the logistic model compared to the linear model, mutatis mutandis for the S = 0 range of SM data (see Figure 6). This seemingly contributed to explaining the better performance of the logistic model for most of the Smax values tested compared to that of the linear model. Furthermore, as shown in the example presented in Figure 6, the logistic model apparently fit the data better within the range 0 < S < Smax, and within this range, the logistic function depicted a quasi-linear shape.
Furthermore, the logistic model seemingly better represented S's behaviour in its lower range. In fact, the Thornthwaite-Mather method [83] applied in this study for computing the water balance uses an asymptotical approach to soil water depletion during the drying periods, which rarely results in null values of soil water storage. The linear model cannot represent this behaviour. Apart from the non-linearity in S, non-linearity was also reported in the SM remote appraisal. For instance, non-linearity was observed in the relationship between the in situ soil moisture measurements and the soil dielectric constant for low soil water contents, while for most of its range of application, the relationship was linear [69].
It should be noted that the shape of the S vs. SM relationship is not physically defined yet, meaning that only empirical interpretations can be derived in favour of the non-linearity. Nevertheless, and for the same reason, recommendations for using one of the models and discarding the other can only stand on their performance, assessed on a statistical basis. All in all, the logistic model better described the relationship between S and SM and, globally, was better fitted to data than the linear model in this case study, despite its poorer predictive capabilities, as represented by the large SEE.

Effect of Smax
The models were applied with the series of S monthly values obtained from the SWB computations using Smax 25 to 150 mm. The models' performances changed nonmonotonically with Smax, as higher r 2 and CE were obtained for Smax 50 and 75 mm when applying the logistic and linear models, respectively. The model performance indicators decreased either towards the lowest Smax tested or towards the highest one, where the poorest performances were obtained in both the linear and logistic models.
Smax depends on (i) the soil physical characteristics determining the maximum amount of water retained in soil pores that can be extracted by plant roots, and on (ii) the soil depth that is explored by plant roots [84]. The former includes the available water capacity (AWC) and bulk density (BD), with their effect requiring a correction for the presence of coarse particles (>2 mm) in case they were determined in fine soil earth (<2 mm). AWC and BD vary inversely with soil properties, such as texture and organic matter content (SOM). Hence, the product of AWC and BD depicts a narrow range of variation due to soil texture and SOM compared to the much wider range of variation that is associated with the soil depth. This means that Smax depends much more on the defined soil depth than on differences between soil types in terms of what concerns their water retention characteristics.
Despite the wide set of soil types found in the study area, with essential differences in AWC and BD, which are mainly associated with the parent material [79], the Smax per unit soil depth (mm cm −1 ) ranged from 0.12-0.20 (calculations based on data from [71]). Accordingly, for the Smax tested in this study, the soil depths ranged from 20 cm to 1 m.
The models' performances were expected to decrease with Smax. In fact, radar penetration does not exceed a 5 cm soil top layer [22]. SM's sensitivity to maximum storage capacity has already been reported, although not fully assessed for many soil types [62]. Soil water storage, S, is assumed to correspond to a constant moisture value throughout the soil profile. Although this is a reliable assumption for small soil depths, this might not represent the prevailing ground conditions for larger soil depths well. In the latter case, surface soil moisture, represented by SM, may shift from the profile-averaged soil moisture represented by S. Furthermore, this shifting may differ according to the soil wetting or drying phase during the year. Therefore, increasing the soil depth considered as a reference for soil water balance computations may result in a lower model performance, which was actually found in this study. The decline in model performance for the lower Smax tested does not fit in the above elaboration. An Smax as low as 25 mm resulted in a fast filling of soil water storage capacity during recharge by excess precipitation over the reference evapotranspiration. If Smax is surpassed, the SWB computes the surpassed amount as water superavit. As expected, in this study the mean annual superavit steadily decreased from Smax 25 to 150 mm (438 against 334 mm), as well as the mean annual number of months with superavit in the study period (6 against 4.5 months), which was equal to the number of months with S = Smax. Accordingly, in these months, S remained constant while SM varied at its higher range of values, reflecting water surplus accumulation over the ground surface. This may contribute to explaining the lower quality fit of both models to the data for the lowest Smax tested (25 mm).

Soil Water Recharge and Depletion
As shown, the segregation of the soil water recharge (wetting) from the depletion (drying) in the data series improved the logistic model's global performance. Similar S values were estimated from higher SM values in the wetting phase model than those estimated with the drying phase model. This finding is a reminder of the issue of moisture distribution in the soil profile, which is assumed uniform in S. SM reflects soil moisture at its shallow surface layer, which might not correspond to the profile-averaged moisture reflected by S. During the wetting phase in the study area, occurring during autumn and winter, the surface was exposed to frequent rainfalls, regardless of the redistribution rate of infiltrated water within the profile. In contrast, in the drying phase, corresponding to the Mediterranean spring and summer, the surface was exposed to very infrequent rainfalls and high evapotranspiration rates, while the soil storage in the profile declined at a much slower rate than that of water depletion as vapour transfers to the lower atmosphere via the surface soil layer.
The logistic model gave a good representation of this interpretation as an S estimate with the wetting and the drying functions fitted to the data for high and low SM values. The former corresponded to the full profile soil moistening after some months of excess rainfall over the evapotranspiration, while the latter corresponded to the end of the dry season (summer) after a large period of low rainfall amounts and high evapotranspiration rates. The segregation procedure of the two sets of data (wetting and drying phases) was not actually applied for the linear model. In fact, this would result in a non-realistic range of SM values that would produce a constant S estimation at S = 0 and at S = Smax.
The above description represents the hysteretic behaviour of the relationship between SM and S, which is not reported in the literature. In fact, regardless of the source of ground data used to be matched with the satellite-borne data, research has been targeted towards finding a single relationship between both data sets. However, hysteresis has been known for a long time in soil physics and soil hydrology and associated with water flowing through porous media [68,92]. Furthermore, hysteresis is found in surface hydrology, where sediment transport studies in streams show that sediment yields may depict such behaviour, as well as ionic species dissolved or adsorbed to the soiled phase [64,93,94].

A Note on the Models' Applicability
The methodological approach followed in this study was strictly statistical, prospecting direct applications in the agricultural water management field. For the case study presented, the approach aimed at converting a remote non-invasive monitoring soil moisture dataset into an actual volume of water in the soil per unit area, according to a specified soil water storage capacity. This is understood as a managerial improvement in this field. The regression functions obtained are the tools being provided to the potential end-users of these research results. The derived two-variable functions are especially highlighted in this context. As such, the expected context of the models' applications is practical and unsophisticated. Still, care is required in the models' applications for prediction purposes, given the wide range of possible results that the models' SEEs demonstrate.
The research reported on a case study located in NE Portugal, NW Iberian Peninsula. As such, the results obtained can be reliably used in applications in this region, not only for agricultural water management purposes but also for other purposes, including drought monitoring, which is an important issue in Mediterranean agri-environments. The ESA CCI SM combined product is developed on the basis of a worldwide network of ground stations providing in situ soil moisture data for the product's calibration [43]. Local ground data are, therefore, always necessary for improving the product's performance and local studies may add valuable information to that of the reference network of ground stations feeding SM calibration runs. It is important to stress that the central and north-western part of the Iberian Peninsula is covered by a single reference ground station of the SM calibration network [43]. On one hand, enlarging the SM product spatial coverage with increasing performance is a continuous challenge to which the present case study is intended to contribute. On the other hand, taking into account the spatial coverage of the calibration network of ground stations in the region, the present study's results can be considered for application in the NW Iberian Peninsula.
The present research applied meteorological data to derive soil water storage from soil water balance computations. As it is consistently related to SM combined product data, soil water storage is shown to have the potential for in situ SM calibration exercises, which was done indirectly in this case. This is a far-reaching prospect for the application of the approach tested in the present case study. In fact, the weather station networks are much larger and have a much denser spatial coverage than that of in situ soil moisture. The a priori definition of Smax for computing the water balance, as demonstrated, is a core element in the prospected developments of the present case study.

Conclusions
The ESA SM satellite product provides a remote and large-scale assessment of ground soil moisture at a daily time step. The product is a very interesting tool to be further explored in view of its wider applications in support of a range of research fields, as well as in support of environmental and economic sectors requiring such freely available information. For this, SM validation with ground data has to progress by increasing the number of ground stations and point datasets for comparison purposes to improve the reliability and coverage. Case studies exploring the relationship of SM data with point ground data contribute to such a purpose that the present case study entirely aligns with.
On the other hand, applications of SM data have seemingly not yet achieved the level of demand by end users that such a powerful product is expected to stimulate. The research that is focused on translating SM data to data that is useful and directly usable by the potential end users in several application sectors is still progressing, a research line for which the present study is also intended to contribute towards. Machine learning approaches to better define the relationship between S and SM to provide a more robust predictive tool compared to the conventional statistical methods used in this study is certainly a challenge to be undertaken in future research. As for the present, this study aimed at exploring the relationship between SM and soil water storage, as assessed through the soil water balance computation. In this context, it is worth noting that while SM expresses a remotely assessed volumetric soil water content, the soil water balance is expressed in terms of an equivalent water height in the soil, meaning a volume of water in the soil per unit ground area. This scientifically meaningless difference is, however, a significant asset when it comes to its application by end users.
A major conclusion of the research carried out is that the regression models representing the relationship between S and SM showed fairly high performances, with r 2 in the range of 0.8 to 0.9. This is an encouraging result in what concerns the use of strictly statistical approaches to deriving prediction tools for soil water storage from satellite-borne soil moisture data, namely the ESA SM product. However, as this study showed, the predictability was limited by the large data scatter around the fitted functions, as indicated by an SEE around 15% of the soil water storage capacity, which represents a still too large standard error when expressed in absolute terms (for instance, 15 mm or more SEE for Smax 100 mm).
A second major conclusion of the research was that the S vs. SM relationship was not unique, where more similar S's were found for larger SM values when the soil was recharg-ing (wetting phase) than when the soil was depleting water (drying phase). This actually represents a hysteresis in the S vs. SM relationship, which had not been reported earlier; this opens the challenging possibilities for improving the tools used to predict S from SM in a way that accounts for the pattern of evolution of the soil water balance variables along the year. The issue is remarkably important in the Mediterranean climate type, where the long dry season leads to non-uniform topsoil moisture profiles, which is meaningless for the water balance computations but biases SM due to the low penetration depth of satellite radar sensors.
Third, as an important outcome of the research, the parameters of the regression functions, which were obtained for each one of the soil water storage capacity tested, were highly correlated with Smax. Therefore, two-variable models (linear, logistic and segregated) were derived, which represented the soil water storage prediction surface in response to temporally variable SM and spatially variable Smax. Considering the extent of the study period and the level of performance reached (with an r 2 around 0.85), these models are assumed to be calibrated and, hence, available for validation exercises and use in the Bragança area.  Data Availability Statement: Publicly available ESA CCI SM product datasets were analyzed in this study. This data can be found here: [https://cophub.copernicus.eu/dhus/#/home]. The soil water storage data presented in this study was computed with meteorological data recorded at and managed by the Instituto Politécnico de Bragança, which are available on request from the corresponding author. The data are not publicly available because the data have not the format conditions to be released to public access.