Study area and data

The study area comprises the entire Portuguese national territory (mainland and islands). Data from 20 weather stations belonging to the main meteorological network of the Portuguese Institute for Sea and Atmosphere (IPMA, IP) were used. Stations are identified by acronyms and corresponding names on figure 1. The choice of these locations reflected the country's existing climate diversity.

Fig. 1 Study area (Portugal) and studied weather stations belonging to the main IPMA network. 

Table I shows the locations, geographic coordinates and distance to the sea of ​​each weather station. Except for two stations located in the Madeira Archipelago (at around 32-33( North latitude), all are located between 37(01'N (FAR) and 41(49'N (MTG). According to the NUTS II classification (Level 2 of territorial units for statistical purposes) applied to the Portuguese territory, five stations are located to the North, five in the Centro, five in Alentejo, one in Algarve and two in each of the archipelagos (Azores and Madeira). Only GRD and MTG, both northern locations, are located above 1000 meters of altitude. BGC, VRL, VIS, PTG, and CTD are located between 400m and 700m altitude whereas all others are below 400 m. STR is the only non-coastal station located at an altitude of less than 100m. In the mainland, six out of 16 are coastal locations (VCT, AVR, CC, ASL, SNS, and FAR).

Climate types and subtypes based on the Climate Classification of Köppen reflect the territory's climate diversity. For each of the 20 stations studied, table II presents normal values (1971-2000) of the most relevant climatic elements (mean temperature - 𝑇 𝑎 , air vapor pressure deficit - 𝑉𝑃𝐷 𝑎 , rainfall - 𝑅 𝑎 , wind speed at 2m height - 𝑉 𝑎 , relative humidity - 𝑅𝐻 𝑎 and daylight hours - 𝐷𝐿𝐻 𝑎 , followed by the corresponding climate types: Mediterranean climate with cool - Csb, or hot summer - Csa and semiarid climate - BS. Locations with the greatest 𝑇 𝑎 are found in the south and on the islands (only SNS has a 𝑇 𝑎 below 16(C) whereas the coldest stations are found in the north inland (10(C≤ 𝑇 𝑎 <14(C). The lowest 𝑅 𝑎 values were also found in the south, where rarely exceeding 700mm. The rainiest locations ( 𝑅 𝑎 > 1000mm) are in the nothern Portugal (MTG, VCT, VRL, and VIS) and in Azores (SCF and AGH). 𝑉 𝑎 exceed 10kmh^-1 in three out of four island stations (SCF, AGH, and PTS) and in three coastal stations (CC, SNS, and FAR). According to the United Nations Convention to Combat Desertification (^{Cherlet et al., 2018}) all but one location is in the Humid climate-zone (PTS, with BS climate-type, is the exception). Ten locations have Csb climate (four on the continental coast, five in the northern or central-northern interior, and one in the Azores archipelago), eight have a Csa climate (two coastal locations, five inland, mainly in the central and southern regions, and one in the Madeira archipelago), one with a Cfb climate (SCF/Azores) whereas PTS/Madeira with a dry climate (BS).

Table II Normal (1971-2000) annual values of temperature ( 𝑇 𝑎 ), air vapor pressure deficit ( 𝑉𝑃𝐷 𝑎 ), rainfall ( 𝑅 𝑎 ), wind speed at 2m height ( 𝑉 𝑎 ), relative humidity ( 𝑅𝐻 𝑎 ), daylight hours ( 𝐷𝐿𝐻 𝑎 ) climate classification (by Köppen) and for each of the studied stations. 

Reference evapotranspiration (ETo) estimation methods

Monthly and daily values of climatic elements measured at 20 selected meteorological stations were used to estimate Reference Evapotranspiration (ETo). Normal monthly values (1971-2000) were used from all the stations. Monthly and daily values for 2019 and 2020 were used from two locations only (AVR and EVR). All data relating to AVR were obtained from one weather station whereas records for EVR were provided by two different stations (one located in the city centre, operational until the beginning of this century and the other on the outskirts of the city, providing data afterward).

Both monthly and daily ETo were estimated by PM and HS methods. As recommended by FAO, the former is taken as the standard method (^{Allen et al., 1998}):

(1)

where R_n is the net radiation at crop surface (MJ m⁻² day⁻¹), G is the soil heat flux density (MJ m⁻² day⁻¹), T is air temperature at 2 m height (^◦C), u₂ is the wind speed at 2m height (m.s⁻¹), e_s-e_a is the saturation vapour pressure deficit (kPa) (saturation vapour pressure minus actual vapour pressure), ( is the slope of the vapour pressure curve (in kPa (C^-1) and ( the psychometric constant (kPa ◦C⁻¹). To solve Eq. (1), standard meteorological records of solar radiation or sunshine duration (daylight hours), minimum, maximum and mean air temperature, air humidity and wind speed were required. Whether for calculating different parameters or overcoming a lack of data, the recommendations proposed by ^{Allen et al. (1998}) were followed. All requirements for meteorological data approval have been met. Latitude, altitude and height at which the wind speed was measured were the metadata considered.

The ^{Hargreaves-Samani equation (Hargreaves & Samani, 1985}) was recommended by ^{Allen et al. (1998}) whenever data required in the PM equation (moisture, wind, and/or radiation) are missing:

(2)

where T, T_max and T_min (all expressed in (C) are, respectively, the mean, maximum and minimum temperatures of the time base used (daily or month), Ra is the extra-terrestrial radiation at the latitude, month and hemisphere considered (MJ m^-2 day^-1), ( is the latent heat of vaporization (MJ Kg^-1) for the mean air temperature T ((C), Krs ((C ^0.5) is the radiation coefficient (=0.17(C ^-0.5, by default) and 0.0135 is a factor for conversion from US to international units.

Analytical procedure

Annual, monthly and daily ETo values obtained by applying both PM and HS (with Krs=0.17) methods were compared. Linear regressions (y=b_o+b₁x, where y is ETo_(PM), x is ETo_(HS), b_o and b₁ are the intercept and slope, respectively) established between series of monthly or daily ETo values estimated by the two methods for each of the twenty locations. Goodness-of-fit was assessed using R² (coefficient of determination), RMSE (the root mean square error) and RE (relative error):

(3)

(4)

(5)

where ET_o(PM)I and ET_o(HS)I (i = 1, 2, . . . , n) represent pairs of values of ETo estimated using PM and HS equations, respectively, 𝐸𝑇 𝑜(𝑃𝑀) and 𝐸𝑇 𝑜(𝐻𝑆) are the respective mean values and n_o is the num_ber of days used in the assessment.

Other regressions with lines forced to pass through the origin of the axes (y=b₁x, where y, x and b₁ have the same meaning) were also established between the same data series. Parameters of goodness-of-fit (R², RMSE, and RE) were recalculated. Their values, naturally less favorable, are not comparable with those for the two-parameter lines. HS equation was then calibrated against PM equation. For this purpose, Krs coefficient was recalculated by multiplying the standard value (0.17) by the slope of the trend (straight) line obtained by regressions ETo_(PM) = b₁ ETo_(HS). This procedure was suggested by ^{Rodrigues and Braga (2021}). The use of new Krs values should allow the HS method to replace the PM, whenever the scarcity of meteorological or climatic data is an obstacle and the correlation between the values obtained by the two methods is high. The Analysis ToolPak (Excel) software was used to analyze data related to the regressions performed.

The results obtained by comparing the two methods (differences, correlation, etc.) were evaluated and discussed based on the climatic factors that mostly influence the climate in Portugal: proximity to the sea (here measured by mean annual temperature ranges referred to 30 consecutive years) and elevation. Wind speed was also considered as a factor, either due to the importance it has in ETo calculations using the PM method, or due to the relevance given to several studies found in the literature, as mentioned above.

Annual and monthly reference evapotranspiration (ETo)

1.1. Annual values of evapotranspiration (ET)

The annual ETo values obtained by applying the two methods from normal values (1970-2000) are shown in table III. Values estimated by PM method for the 20 stations considered ranged from 781mm (CC) to 1191mm (EVR) whereas those estimated by HS method ranged from 667mm (CC) to 1256mm (ASL). Unsurprisingly, the greatest values for ETo (above 1000mm per year) obtained by both methods were generally found in regions with hot summers (Csa) whereas the lowest (less than 950mm) were recorded where summer is cool (Csb and Cfb). Even so, this trend is more visible in the results obtained by the PM method than by the HS method.

In half of the locations studied (MTR, BGC, VCT, MDL, VRL, VIS, STR, ASL, CTD, FNC) ETo_(HS)>ETo_(PM). ETo_(HS)<ETo_(PM) in the others. Absolute differences were relevant (> 100mm) in VRL, VIS, CC, ASL, CTD, STR, SNS, and PTS, and only residual (< 5mm) in MTR and FNC. These differences seem to be relatively independent of the climate type/subtype defined by Köppen for each of the 20 locations, that is, the maximum monthly temperature in the summer period (Csb vs. Csa) does not seem to determine the sign of the difference in the results of each method. In fact, ETo_(HS)>ETo_(PM) in half of the stations with Csb climate and in five out of eight locations with a Csa climate. Neither precipitation regimes nor annual precipitation seem to influence the direction of such differences: e.g., in both the wettest station (SCF, with Cfb climate) and in the driest (PTS, with BS climate), ETo_(HS)> ETo_(PM).

Figure 2 highlights the combined influence of mean annual thermal amplitude ( 𝑇𝐴 𝑎 ), mean annual wind speed ( 𝑉 𝑎 ) and elevation on the differences between the annual ETo values obtained by the two methods. In general, ETO_(HS)< ETO_(PM) whenever 𝑇𝐴 𝑎 are smallest (usually in coastal or insular areas) and 𝑉 𝑎 are greatest (≥ 2-2.5ms^-1). Otherwise, the values obtained by the HS method overestimated those found by the PM method whenever 𝑇𝐴 𝑎 are greatest (generally > 10(C) and 𝑉 𝑎 are smallest (< 2 ms^-1 in most cases). Only in three coastal locations (FNC, ASL, and VCT), all with 𝑉 𝑎 < 2ms^-1, ETo_(PM) was not underestimated. In ASL and VCT (with 𝑇𝐴 𝑎 similar to those observed in most inland locations) there was even a clear overestimate of the ET_O calculated by the PM method. The only inland locations where ETo was underestimated, were those with 𝑉 𝑎 > than 2-2.5ms^-1 (GRD, EVR, and PTG). Finally, elevation does not seem to have the same influence as each of the factors above does (locations at approximate altitudes show different trends). When elevation was evaluated as a factor associated only with atmospheric pressure variation (and not associated with temperature variation), ETo estimates using the PM method were based on sea level (( = 0,067 kPa◦C⁻¹). In this context, the natural and consequent increase of ETo would range from about 15-17mm at MDL, VRL, VIS, CTD, and EVR (all located between 200 and 500m) to about 40mm at GRD and MTG (both located at about 1000m). If in these last two locations, this would mean an underestimate of ETo (PM) obtained by HS method and a relevant approximation between values in the cases of BGC and PTG.

Fig. 2 ETo(HS) - ETo(PM) values as a function of wind speed and thermal amplitudes (mean annual values) for the 20 meteorological stations. 

Figure 3 shows the influence of each of these factors (thermal amplitude, wind speed and elevation) on annual ETo_(HS)-ETo_(PM). For larger 𝑇𝐴 𝑎 reference ETo is overestimated whereas smaller ones underestimated it (fig. 3a). Despite the level of significance of the regression (p-value <0.05), the goodness-of-fit is just fair (R²=0.68). For thermal amplitude values close to 8(C, the differences were minimal. The results obtained by the two methods are very different for locations near the sea level, tending to be closer as the station elevation increases (except for VIS) (fig. 3b). Differences did not exceed 50 mm for stations located at more than 500m (4 out of 20). At elevations < 100mm, | 𝐸𝑇 𝑜 𝐻𝑆 − 𝐸𝑇 𝑜(𝑃𝑀) | ranged from 5mm to 161mm. The relationship between| ??𝑇 𝑜 𝐻𝑆 − 𝐸𝑇 𝑜(𝑃𝑀) |, and 𝑉 𝑎 was well-described (R²= 0.782) by a straight line (with a negative slope). The statistical relevance of this trend is also significant (p-value < 0.1) (fig. 3c). This means that the HS method tends to underestimate the values obtained by PM under windier conditions and overestimate them under less windy conditions. The 𝑉 𝑎 corresponding to a zero-difference was 2.61 ms^-1 (a reasonable criterion to distinguish, for this purpose, windy locations from non-windy locations).

Fig. 3 ETo_(HS)- ETo_(PM) vs.: (a) Mean annual thermal amplitude ( 𝑇𝐴 𝑎 ); (b) elevation and (c) mean annual wind speed ( 𝑉 𝑎 ) for the 20 weather stations studied. 

The annual ETo_(HS) values ​​obtained by both methods from data recorded in AVR and EVR during 2019 and 2020 were also different (data not shown). As with normal data, those obtained in AVR by ET_o(HS) underestimated (by approximately 15-20mm) those found by the PM method. On the contrary, ET_o(HS) > ET_o(PM) in both years, with differences greater than 100mm. The relocation of the weather station in EVR (to a less windy location) could explain this fact, reinforcing the importance of wind speed in the differences between the ETo obtained by the two methods.

1.2. Monthly values of Evapotranspiration (ETo)

The variations in ETo values throughout the year estimated by both methods largely reflect the central influence of temperature on the evapotranspiration process: minimum values in December/January and maximum values in July/August at these latitudes. Winter minimum values of 20-30 mm are common in most places (ASL, FAR, and insular locations were exceptions). In summer, ETo often reached monthly values close to 200mm in the inland, mainly south, whereas on the islands and on two coastal locations (CC and SNS) they never exceeded 120mm. The differences ETo_(HS)- ETo_(PM) also vary throughout the year (fig. 4), generally reaching the greatest magnitudes during the hottest period. The exceptions were PTG, FAR, and SCF, where the differences were more noticeable in the autumn-winter period, and AVR, where they never exceeded five mm.

Monthly ETo_(PM) values were always lower than ETo_(HS) in VCT, VRL, VIS, ASL, CTD, and STR (the six least windy locations) and always greater at CC, SNS, FAR, and PTS (four of the windiest locations) (eg., fig. 4a and 4b, respectively) . In the other locations, the differences are negative in one part of the year and positive in another (eg., fig. 4c). ETo _(HS) - ETo _(PM) exceeded 20mm in at least one month at MDL, VRL, VIS, STR, ASL, and SNS (eg., fig. 4d). At MTG, AVR, GRD, PTG, FAR, AGH, FNC, the differences between monthly values never exceeded 10 mm (eg., fig. 4e).

Fig. 4 Differences between the monthly average ETo (normal) values obtained by the HS and PM methods, for four locations. 

2. Correlations between monthly or daily values obtained by HS and PM methods

2.1 From normal monthly values

Monthly values of ETo_(PM) were plotted against those of ETo_(HS). The regression equations and goodness-of-fit parameters are shown in table IV. The linear correlations were very significant in all cases (p-value <0.01) and the goodness-of-fits were very high (R² was around 0.98 and 0.97 in the island locations and SNS, and ≥0.99 in the others). Although unsurprisingly the RMSE and RE values ​​in most cases increased significantly when the line was forced to cross the origin, the R² values ​​remained quite high (only in SCF R²<0.97). Decreases were greater than 50% in 10 locations (VCT, VRL, VIS, STR, ASL, CTD, SNS, MDL, CC, and PST) and only residual (<5%) in 5 (GRD, SCF, PTG, AGH, and FNC).

When in the equation forced to pass through the origin the b₁ values were > 1 (VCT, MTL, VRL, MDL, BGC, VIS, STR, ASL, CTD, FNC, and especially in VIS), the Krs used (0.17(C^-0.5) proved to be excessive; when b₁ <1 (PTG, EVO, FAR, PTS, SCF, AGH, and especially in SNS, and CC) the Krs values were naturally deficient. In AVR and in GRD, the Krs value proved to be adequate. The HS method can then be calibrated by applying the new Krs values to the HS equation. For the places farthest from the sea, Krs values ranged from 0.155(C ^-0.5 (MDL) to 0.177 (C^-0.5 (EVR) whereas on the mainland coast, the range of values was wider (from 0.149(C ^-0.5 at ASL to 0.201(C ^-0.5 at CC). At island stations, corrected Krs values ranged from 0.169 (C ^-0.5 (FNC) to 0.188 (C^-0.5 (PTS). The lowest value was registered in VIS (0.14(C ^-0.5) (fig. 5).

A significant (p-value <0.05) and an acceptable (R²=0.69) linear correlation was found (fig. 6a) when the influence of 𝑇𝐴 𝑎 on the variation of the Krs coefficient was tested. The quality of regression fit was considerably improved when only mainland locations were considered (R²=0.81). In a hypothetical isothermal environment (not verifiable in this climate-type), Krs would be between 0.21 and 0.23, depending on the established regressions. As thermal amplitudes increase (increase supposedly associated with increasing distance from the sea) the value of Krs decreases. Adjusted Krs values in insular territories are visibly lower than those obtained in coastal sites with similar thermal amplitudes. In addition, W-E transects without a coherent trend (fig. 3a) suggest that other factors may also be required to better calibrate Krs.

Fig. 5 Recalibration of the Hargreaves-Samani equation (modified Krs values). 

Also a significant (p-value <0.05) and acceptable (R²=0.72) linear correlation was found when Krs values are plotted against 𝑉 𝑎 (fig. 6b). When only mainland locations were considered the goodness-of-fit also increased (R²= 0.87). The best fit between Krs and wind speed is not surprising, as this already explained the variation of ET_O(HS)-ET_O(PM) (annual values) better than any other parameter. For mainland locations, in the absence of wind, Krs=0,123(C^-0.5 (0.113 < Krs < 0.133, with 95% of confidence), a useful value for predicting Krs as a function of wind speed (T-test for (=0.05). Finally, when the effect of atmospheric pressure is removed, the Krs values increase very slightly (only at the two highest locations do the increments exceed 0.005(C^-0.5).

If a multiple linear regression is used, with 𝑉 𝑎 and 𝑇𝐴 𝑎 as independent variables (predictors), more significant (p-value <0.05) and better correlations (R² equal to 0.79 or to 0.91 when all locations or only mainland locations were considered, respectively) were found (fig. 7). P-values for both independent variables were < 0.05together then. In short, Krs is better explained by both variables jointly. Table V presents Krs values (approximated to 0.005(C^-0.5) for different combined values of 𝑉 𝑎 and 𝑇𝐴 𝑎 (both approximated to units).

Fig. 6 (a) Krs-values = f ( 𝑇𝐴 𝑎 ) and (b) Krs-values = f ( 𝑉 𝑎 ) for 20 weather stations in Portuguese territory. R² - coefficient of determination. 

Fig. 7 Krs-values estimated (by multiple regression) as a function of both (normal) annual thermal amplitude and wind speed values, for 20 weather stations in Portuguese territory. R² - coefficient of determination. 

Table V Krs values (in (C ^-0.5; approximation: 0.005(C ^-0.5) fitted for different combination of mean annual thermal amplitude (TA 𝑎 ) and mean annual wind speed ( 𝑉 𝑎 ), using the equation: Krs = 0.00892 𝑉 𝑎 -0.00253 𝑇𝐴 𝑎 +0.16785. 

2.2. From monthly values recorded in 2019 and 2020

Correlations between the monthly ETo values estimated by the two methods for 2019 or 2020 were also good, whether the trend lines found has been forced to pass through the origin or not (R² was always greater than 0.97) (fig. 8). Correlations between the average values for AVR referred to these two years were not significantly different (p≥0.05) from those found when normal values were used, but those for EVR were significant (p<0.01). Consequently, similar Krs values were obtained (0.17-0.171 in any case) in AVR, but not in EVR (Krs ( 0.16 in the two years studied and around 0.18 for the 30 years considered).

2.3 From daily values recorded in 2019 and 2020

Correlations between the daily values for to AVR and EVR obtained by both methods were also significant (p-value <0.01). Goodness-of-fit was also high (R²(0.90) for any of the years considered (2019 or 2020), allowing us to use the HS method instead of PM on a daily basis as well (fig. 9). RMSE for AVR ranged from 0.37 mm day^-1 (2019) to 0.41 mm day^-1 (2019) whereas for EVR it ranged from 0.77 mm day^-1 (2020) to 0.80 mm day^-1 (2019). When the regression line was forced to pass through the origin (b_o=0), R² did not drop substantially (also around 0.90). The decrease in RMSE resulting from using this relationship to estimate ETo_(PM) from ETo_(HS) was around 10% in EVR for 2020. The default value of Krs used by default (=0.17) was overestimated in EVR ((0.16) whereas in AVR it seemed to be adequate. Estimates for AVR were similar to those found for monthly estimates (from both normal data and the data for 2019 and 2020), whereas for EVR they were clearly smaller than those obtained when monthly data were used.

When comparing the daily values separately for each of the two semesters (April-September and October-March) the R² values decreased substantially in both locations, to around 0.78 in both semesters in the case of EVR, and to 0.84 and 0.59, respectively in the hottest and coldest semester, in the case of AVR (data not shown) The Krs for both semesters were similar in AVR ((0.17), but not in EVR (0.164 for the warmer semester and 0.139 for the other).

Fig. 8 ETo_(PM) vs. ETo_(HS) (linear regressions) with monthly values for 2019 (() and 2020 (() in (a) AVR and (b) EVR. Only the straight lines y=b₁x+b_o are shown: 2019 ( () and 2020 ( ^{_ _ _} ). R² - coefficient of determination. 

Fig. 9 ETo_(PM) vs. ETo_(HS) (linear regressions) with daily values for 2019 (() and 2020 (() in (a) AVR and (b) EVR. Only the straight lines y=b₁x+b_o are shown: 2019 ( () and 2020 ( ^{_ _ _} ). R² - coefficient of determination.