Random forest predictive model development with uncertainty analysis capability for the estimation of evapotranspiration in an arid oasis region

Evapotranspiration (ET) is the process of transfer of water from the surface of the earth to the atmosphere including evaporation and transpiration (Shiri et al. 2014; Nourani et al. 2019), and often used to estimate actual evapotranspiration in water balance studies and water resources management (Tao et al. 2015). In arid oasis conditions, crops are a material basis on which human beings depend for their survival as well as being an ecological protection barrier in such areas. Knowledge of crop-water demands is an important practical consideration for improved water-use efficiency (Benli et al. 2006). This is because ET is a primary source of water loss, so its accurate evaluation can provide valuable information for water balance, irrigation system design, and water resources management (Torres et al. 2011; Wen et al. 2015). This is especially true for arid regions, such as the northwest region in China, where population growth, expansion of agriculture, and other socio-economic activities are significantly constraining the available water resources.

Due to the lack of observation data, the precise estimation of ET has produced the need for another comprehensive concept called reference evapotranspiration (ET₀) (Abdullah et al. 2014). ET₀ can be measured directly using lysimeters which are characterized by providing accurate measurement results; however, the application of the methods is limited due to their cost and complexity (Ferreira et al. 2019), which increases the requirements of employing data-based methods to predict ET₀. Several conventionally empirical models like Hargreaves equation, Priestley–Taylor equation, and Ritchie equation have been developed to estimate ET₀ using meteorological data. Because the PMF-56 equation takes into account moisture availability, mass transfer, and required energy for the process (Granata 2019), it has been recommended for the computation of ET₀ by the Food and Agricultural Organization of the United Nations (FAO) as the only standard equation which is usually applied to validate other models and has been accepted in many regions across the world. PMF-56 equation can be broadly applied in various environments and climate conditions due to its good precision and stability (Huang et al. 2019). However, some restrictions still exist in the application of PMF-56 equation, for example, it is difficult to obtain all meteorological data required in the estimation process, particularly in a developing country, where the number of meteorological stations is limited and weather data records could be scarce (Abdullah et al. 2015). Within this context, an alternative data-driven model which requires easily available input variables is necessary and significant.

As the ET₀ depends on several interacting meteorological factors, such as temperature, humidity, wind speed, and radiation, it is difficult for the ordinary formula to express all the related physical processes (Yassin et al. 2016; Yin et al. 2016). In this context, artificial intelligence or data-driven models are considered as efficient tools to deal with non-linear relationships between independent and dependent variables. In the past few decades, artificial intelligence models, including artificial neural network (ANN), extreme learning machine (ELM), support vector machine (SVM), and so on, have been extensively used in the area of predicting and forecasting (Kisi & Cimen 2009; Yoon et al. 2011; Tabari et al. 2012; Acharya et al. 2013; He et al. 2014; Deo & Şahin 2015). In terms of ET₀ prediction, Traore et al. (2010) assessed the performance of feed forward backpropagation neural network (BPNN) algorithm (a type of ANN) based on different inputs in estimating ET₀ in the Bobo-Dioulasso region. The results showed that the BPNN algorithm had a better performance than conventional Hargreaves equation and that wind was found to be the most effective variable significantly required for modeling with high accuracy when added into inputs. Huo et al. (2012) compared the performance of ANN models with multiple linear regressions, the Penman equation, and two empirical equations for calculation of ET₀ in northwest China, concluding that ANN models exhibited higher accuracy than the others, and they also concluded that temperature, R_h, was the most important input affecting ET₀. Abdullah et al. (2015) proved that ELM was efficient, simple in application, of high speed, and had a very good generalization performance at predicting Penman–Monteith (P-M) ET₀ using four different complete and incomplete meteorological input combinations in Iraq. Patil & Deka (2016) developed the ELM model utilizing three different input combinations to calculate ET₀ in the Thar Desert, India, and Hargreaves equation, ANN and least-square support vector machine (LS-SVM) models were used for a contrast. The results revealed that ELM is a simple yet efficient algorithm and superior to the other two methods. Tabari et al. (2012) estimated the performances of SVM, adaptive neuro-fuzzy inference system (ANFIS), multiple linear regression (MLR), and multiple non-linear regression (MNLR) for estimating ET₀ using six input vectors of climatic data in a semi-arid highland environment in Iran. The results displayed that the capability of SVM and ANFIS models for ET₀ prediction was better than those achieved using the regression and climate-based models. Kisi & Cimen (2009) used the SVM approach for modeling ET₀ in three stations in central California. The results were compared with empirical models and ANN model and revealed that the SVM method could be employed successfully in simulating the ET₀ process. These models have demonstrated promising prediction ability of ET₀ in many parts of the world, but some deficiencies exist. ANN models become easily stuck in a local minimum, and the optimization process is effortlessly influenced by initial point selection. SVM and numerous ELM models are machine learning methods based on kernel function, and generalization abilities depend largely on the choice of the kernel function.

Random forest (RF) is another emerging machine learning technique and a natural non-linear modeling tool, the superiority of which is good tolerance for outliers and noise, difficulty in producing an over-fitting phenomenon. As well, it can overcome the ‘black-box’ limitations of ANN and provides evaluation of the importance degree of input variables (Rodriguez-Galiano et al. 2014). RF with its merits has been widely used in classification and prediction (Gislason et al. 2006; Cutler et al. 2007; Heung et al. 2014; Gong et al. 2018). Wang et al. (2015) proposed the RF model to evaluate flood hazard risk and implemented the method in Dongjiang River Basin, China; consequently, the capacity of the RF model was similar to the SVM model with a correlation coefficient of 0.916, but the RF method had a better performance with its advantages including providing credible assessment consequences of importance degree of input variables. Dong et al. (2013) classified whether rockburst will happen and the intensity of rockburst in underground rock projects utilizing RF method, and selected some main control factors of rockburst, including the values of in-situ stresses, uniaxial compressive strength and tensile strength of rock, and the elastic energy index of rock to analysis. The results indicated that the RF model exhibited high classification accuracy compared with the ANN and SVM approach with misjudgment ratios of 0, 10%, and 20%, respectively. In RF modeling of ET₀, Fukuda et al. (2013) accessed the applicability of RF model for estimating mango fruit yields using 10-day rainfall and irrigation data in response to water supply under different irrigation regimes. The RF models accurately estimated the maximum and mean values of mango fruit yields, and the results displayed the applicability of RF in the field of agricultural engineering. Feng et al. (2017) proposed RF and generalized regression neural networks (GRNN) models for daily ET₀ estimation in southwest China, and the result revealed that the RF model was slightly better than GRNN model for estimating daily ET₀. Although the RF model demonstrated significant potential in many studies, the use of RF model for evaluating ET₀ has been rarely recorded by research, especially in the arid environment of northwest China. It is, thus, important to predict ET₀ using the RF model to provide a reliable method in data-limited areas.

Despite these advantages, there is a deficiency in the application of the RF model for ET₀ predictions. Almost all the artificial intelligence models are stochastic algorithms, the RF approach is no exception, and running the model will not reproduce the same result even in an identical situation. Uncertainty analysis is an indispensable procedure for getting reliable results in model simulations. For uncertainty analysis, two primarily different aspects of uncertainty include uncertain input variables, model parameters, and model structure. By means of its general applicability, the Monte Carlo simulation technique is a widely used method for uncertainty analysis in hydrological modeling (Shrestha et al. 2009; Antanasijević et al. 2014). However, one remarkable issue is that the uncertainty of the model in estimation is usually ignored by most studies, and no such studies have been reported adding uncertainty analysis in predicting ET₀ so far. In this condition, uncertainty analysis is conducted in the paper for assessing the precision of the RF model.

The present study was carried out in an arid oasis area of the middle reaches of the Heihe River Basin, northwest China (Figure 1), where water resources play an important role in the sustainable development of the ecological environment. Besides, the study area is a typical irrigated agricultural area as well as an important commodity grain base of Heihe River. Agriculture consumes most water, plus water resources are in severely short supply in this region. As a vital component to describe the hydrological cycle, estimate water balance, and schedule irrigation (Rawat et al. 2019), ET₀ determination with reliable accuracy is significant in such a water-scarce region (Nourani et al. 2020). Hence, the objective of this paper was to investigate the precision of the RF model by using different variables' combination of meteorological data including T_max, T_min, U₂, R_h, S_un in an arid region, northwest China. The results obtained from the RF modes for various input combinations are compared to each other, and subsequently determine the effects of different meteorological arguments on ET₀ according to the importance degree of variables. Moreover, the uncertainty analysis is performed for the RF model by Monte Carlo simulations for the purpose of a better accurate result applying to arid areas.

The RF was developed by Breiman (2001) based on a CART decision tree model, including regression (RFR) and classification (RFC) algorithm. The basic idea based on statistical theory is that extracting repeatedly and randomly K samples from the original training sample set N for generating a new set of training samples through the bootstrap resampling method, then producing K decision trees and comprising random forest according to the bootstrap sample set. In terms of the classification model, the classified results of new data depend on the number of votes obtained by classification tree votes, and for the regression model, all the averages of the predictive value of decision trees are regarded as final prediction outcomes (Figure 2).

This paper uses the regression algorithm whose calculation processes are as follows.

First, randomly generate k training samples (Θ₁, Θ₂, …,Θ_k) from the total training sample using the bootstrap sampling method, corresponding to K decision trees can be constructed.

Second, at each node of the decision tree, the m features are randomly selected from the M features as the splitting features set of the current nodes, then selecting one node from the m features to split according to the principle of node purity minimum, each decision tree is grown to the largest extent possible, no pruning.

In addition, index importance assessment is a prominent advantage of the RF algorithm, the purpose of which lies in evaluating the effect of each variable on the accuracy of the RF model. IncNodePurity index adopted in this research was used to assess the importance of each parameter, and compare that by calculating the reduced values of impurity of the nodes of all tree variables. That higher index importance measurement can intuitively reflect the main factors affecting estimated ET₀. Besides, in this research, we applied the randomForest package to train data and access variable importance in the R environment.

The weather data for this study were obtained from two sites in Zhangye (100°17′E, 39°05′N) and Gaotai (99°50′E, 39°22′N), as shown in Figure 1. In this study, five years of meteorological data was sourced from the National Climatic Centre of the China Meteorological Administration. The duration of the data is from 2013 to 2017 at daily timescales, which includes T_max, T_min, U₂, R_h, S_un. There were 1,826 records and these were divided into two parts: the training part composed of 1,461 daily records which account for about 80% of the total data set, and the testing part, the remaining 365 records, which accounts for about 20% of the total data set. The statistical characteristics of daily weather data and the PMF-56 ET₀ for each station are shown in Table 1. In terms of the skewness values, T_max, T_min, and R_h showed lower skewed distribution than other variables. Also, it can be seen that U₂ shows a higher skewed feature than the other variables (1.06 and 1.21 for the two sites, respectively). T_max, T_min, and PMF-56 ET₀ demonstrate a strong variability and the CV values exceed 0.62, which principally resulted from seasonal changes. The variation intensity of the rest of parameters is intermediate (CV is bounded between 0.25 and 0.75). Additionally, there are no significant differences for other data between the two weather stations.

The selection of appropriate input variables has a direct impact on the performance of the model; moreover, finding suitable inputs can provide an efficient way of estimating ET₀ for many regions where weather data are not always available. For the development of the RF model, this study selected different combinations of various daily climatic data as input, and ET₀ computed by daily PMF-56 equation as output for training and testing the models. Eight different combinations were considered in the present study and are referred to in the short form as shown in Table 2. Temperature is the most influential variable on ET₀ and predominant physical factor in the evaporation process (Jain et al. 2008; Wen et al. 2015). Thus, combination 1, as the base inputs, consists of T_max and T_min and the other combinations are formed by integrating S_un; R_h and U₂; R_h, U₂ and S_un; R_h, U₂, and S_un into combination 1, respectively. Each combination was trained and tested by the RF model.

The performance of RF model for PMF-56 ET₀ applied to the studied stations for the training and testing periods are summarized in Tables 3 and 4, which demonstrate the precision of the proposed RF model by the formulae of r, RMSE, MAE, and NS. As can be seen, there were no significant changes in respect to all of the metrics of these models in training as well as testing periods. This research selected the criteria during the testing phase to compare the capabilities of these models in the prediction of PMF-56 ET₀ and all of the following analyses were performed in the testing period.

Considering all models for the two stations, it can be observed that the RF8 model outperforms all of the other models in the four-estimation norm, with the highest r and NS values as well as the lowest RMSE and MAE values. We also clearly see all of the r and NS values surpass 0.8 for both stations, indicating the performances of RF models in PMF-56 ET₀ prediction were encouraging. Thus, it was selected as the best-fit model for estimating the PMF-56 ET₀ at the two stations. In the remaining models, RF5, RF6, and RF7 models including four parameters as inputs had a higher r and NS values, lower RMSE and MAE values and were found to be better than RF2, RF3, and RF4 models with three input parameters at each individual site. Alternatively, RF1 with only T_max and T_min as inputs had the biggest errors rates compared to other models. This demonstrated that the performance of the models relied on the number of input parameters. However, weather factors were usually incomplete in data-limited regions, especially in arid environments, such as northwest China. The selection of the model should be decided according to the available meteorological parameters. The models whose input comprised T_max and T_min are needed and can be used in this study for practical application.

Concretely, in terms of Zhangye station, the corporation of R_h can significantly improve the performance of RF models. Adding R_h to temperature-based inputs, the RF4 model improved r, RMSE, MAE, and NS by 6.2%, 37.4%, 39.4%, and 13.1%, respectively. Likewise, the RF6 and RF7 models introducing R_h as input variable achieved higher simulation precision (with the higher r and NS values, the lower RMSE and MAE values) than the RF5 model with the absence of R_h. From these results, it was shown that the addition of R_h was more sensitive to output relative to the S_un and U₂. The RF2 model performed the second best in ET₀ estimation among the RF2, RF3, and RF4 models. Note that RF3 including U₂ on the basis of RF1 improved r, RMSE, MAE, and NS value by 4.1%, 22%, 23.7%, and 8.5%, respectively. It was an objective fact that RF6 was superior to RF5 according to four evaluation criteria. This result showed that inserting R_h was more effective than U₂ to the estimation of ET₀. The RF3 model, whose inputs included T_max, T_min, and U₂ were found to be worse than RF2 and RF4 models among the three models. As a result, the ET₀ is most easily affected by R_h, followed by S_un and U₂. This conclusion is in disagreement with the findings of many studies (Dai et al. 2009; Petkovic et al. 2015; Tao et al. 2015; Xing et al. 2016), namely, the S_un is considered as the most effective parameter for simulating PMF-56 ET₀. Generally, the results depend on the selected geographical location and climate type of the study area.

For the case of Gaotai station, there were different results compared with those of Zhangye station. It was shown that the RF3 model performed slightly better than RF2 and RF4 models in terms of four statistical indicators, and it can be stated that PMF-56 ET₀ was easily influenced by U₂. This was also confirmed by RF5 and RF7 models with the insertion of U₂ into the inputs presented in Tables 3 and 4. The RF5 and RF7 models remarkably increased the r and NS values of 0.6% and 1.1%, and 0.3% and 0.5%, respectively, and decreased the RMSE and MAE values of 8.4% and 6.5%, and 3.9% and 6.5%, respectively, relative to the RF6 model, exhibiting the superiority of RF5 and RF7 models to the RF6 model significantly. The results of this comparison revealed that integrating U₂ improved the accuracy of the model significantly. Accordingly, adding U₂ is found to be more influential than S_un and R_h on ET₀ simulation, which is the same outcome obtained by Traore et al. (2010) and Karimaldini et al. (2012). It is observed that the input scenarios listed in Table 2 have a distinct performance for the two stations due to the different geographical locations.

To compare the performance of the temperature-based models and the other models with the absence of temperature, the Supplementary material lists performance statistics of four input combinations including: (1) S_un and U₂; (2) S_un and R_h; (3) U₂ and R_h; (4) S_un, U₂ and R_h, and the four inputs are expressed as RF9, RF10, RF11, and RF12, respectively. It is apparent that all the RF models produced higher RMSE (more than 1.38) and MAE (more than 1.08) as well as lower r (less than 0.81) and NS (less than 0.66), and were inferior to combinations 1–8 inserting T_max and T_min into inputs for PMF-56 ET₀ forecasting (Supplementary material, Tables A1 and A2). As the best fitting models, RF12 had r values of 0.809 and 0.790, RMSE values of 1.427 and 1.382, MAE values of 1.097 and 1.082, and NS values of 0.654 and 0.619 for Zhangye and Gaotai stations, respectively, which cannot meet the prediction standards of PMF-56 ET₀. In such circumstances, RF9–RF12 models should not be selected as techniques to estimate PMF-56 ET₀. Therefore, the following does not elaborate on the four models, but mainly focuses on RF1–RF8 models.

Figures 3–6 exhibit the hydrograph and scatter plots of the ET₀ values computed by the PMF-56 equation and the values estimated by different combinations of the RF model of the validation period for the two stations. A total of eight combinations of RF model displayed a good prediction of ET₀. In addition, it is obviously seen that the ET₀ values estimated by the RF8 model were closer to the PMF-56 ET₀ values and followed the same trend than the other models while the RF1 model performed the worst in this area. From the fit line with the form of y = ax + b, the coefficients a and b of the RF8 model were closer to 1 and 0 than the other models, because the lowest values of b (equal to 0), and the highest values of the slope (equal to 1) denote the best fit of models. These were confirmed by r, RSME, MAE, and NS values shown in Tables 3 and 4. As well, it was observed that the fitting performance of the maximum and minimum PMF-56 ET₀ was not very good, especially that of peaks of the first few models.

Due to the importance of PMF-56 ET₀ in irrigation and agricultural water use, water resources planning and management, the estimation of total PMF-56 ET₀ obtained by different combinations of RF model was also considered in this paper. The total ET₀ amounts calculated by PMF-56 and RF models in the testing phase are given in Table 5. It is remarkable that all models had a quite good estimation of total PMF-56 ET₀ value since there was a smaller relative error (all values less than 3.5%) for both sites, especially the RF1 model, whose input parameters were only T_max and T_min at Zhangye station, with a relative error of −0.2%. In addition, noting the fact that the RF8 model with all the variables as inputs did not have the lowest relative error (−1.6%) among all the models at Gaotai station, it still performed well and its value fell within the reliable range. Although generally, reliable weather data sets such as R_h, U₂ and S_un are limited in the arid regions, combining the above calculation results and demands of practical use, the RF model can be employed to predict PMF-56 ET₀ where restricted data are available.

Index importance assessment is an advantage of the RF model which can directly obtain an order of all of the weather parameters. As shown in Figure 7, temperature is the most relevant variable connected with the estimation of PMF-56 ET₀ in the two stations, with IncNodePurity values individually accounting for 69% and 71%, which illustrates that T_max and T_min can be employed to predict PMF-56 ET₀ combined with Tables 3 and 4; on the contrary, the RF model cannot simulate ET₀ with higher accuracy when temperature is missing. In addition, for the importance degree of the other three factors there existed a similarity at the two sites; R_h can be considered as an influential index due to higher IncNodePurity values at Zhangye and Gaotai stations. It is relevant to note that this result is inconsistent with the consequence obtained at Gaotai station (described by Table 4). The aforementioned outcomes are achieved by different combinations among all parameters, indicating that inserting U₂ into inputs has higher precision compared to adding other variables for Gaotai site. Nevertheless, IncNodePurity value of random forest explains each index's contribution to ET₀, therefore, in terms of importance of variables, there is no doubt that R_h is the most relavant factor affecting ET₀ in this area. Similar results were also carried out in Shiyang River Basin, northwest China by Huo et al. (2012), where R_h has a large effect on daily PMF-56 ET₀ except for air temperature in an arid region, northwest China.

The techniques of Monte Carlo simulations were used to corroborate the applicability of RF models in modeling PMF-56 ET₀, an important hydro-meteorological parameter for agriculture, ecosystems, and several other socio-economic activities. The method proposed here has been used to quantify the uncertainty by predicting the confidence intervals of the simulation results. Figures 8 and 9 illustrate 95% confidence intervals for the estimates of daily PMF-56 ET₀ applying the RF model for Zhangye and Gaotai stations during the testing period. From these two figures, we find that there was a good match between 95% confidence intervals and results obtained by the RF model, and most of the observed ET₀ data lay within the confidence intervals at the two stations. Results of Monte Carlo analysis of the RF model for the two stations are given in the upper left corner of the panels. In general, satisfactory results indicate that more observed data were bracketed within the 95PPU (all values are over 60%), while a lower d-factor value can be obtained (d-factor values less than 1 are considered appropriate). Remarkably, the RF8 model produced acceptable d-factor values, noting that the d-factor values were 0.26 and 0.27 at Zhangye and Gaotai station, respectively. Also, the PMF-56 ET₀ values bracketed by 95PPU were more significant for both stations; it was observed that 91% and 87% of the PMF-56 ET₀ data were bracketed by the 95PPU at Zhangye and Gaotai station, respectively. In addition, 95% confidence intervals of RF3 to RF8 models were found to be a relatively good fit to the peak of PMF-56 ET₀ for the two sites. Although the RF1 model has wider 95% confidence and higher d-factor values than other combinations for the two stations, the uncertainty is still within acceptable limits. As can be seen from these figures, the trend of 95% confidence intervals calculated by the RF1 model is basically close to that of PMF-56 ET₀. The values bracketed by 95PPU were more than 62%, and the d-factor values were less than 0.58, indicating that the RF1 model is able to predict daily PMF-56 ET₀ with smaller uncertainties. Considering the purpose of this paper and the above discussion, the prediction uncertainties of the RF1 model were determined to further illustrate the application of this model and therefore find a reliable model in an arid area with a lack of sufficient meteorological data. Besides, it is worth noting that the maximum and minimum ET₀ values cannot be simulated perfectly by models, consistent with the previous conclusion (as shown in Figures 3–6), which also further illustrates the differences between the total ET₀ amounts computed by PMF-56 and RF models as displayed by Table 5. In spite of some errors, taking the discussion of the section on evaluation of the importance of variables into account, we find that the RF model with only T_max and T_min as inputs is still considered as an appropriate technique to simulate daily PMF-56 ET₀ in arid conditions.

Water resources play an essential role in arid environments, so new modeling and water assessment methods are crucial for maintaining sustainability of water resources, strategies for water quality and usage. ET₀ provides a vital parameter of water resources calculation, regional water resources management, and irrigation plan development. This research discussed the performance of the RF model to predict PMF-56 ET₀ using different combinations of daily climatic data, including maximum air temperature (T_max), minimum air temperature (T_min), sunshine duration (S_un), wind speed (U₂), and relative humidity (R_h) for Zhangye and Gaotai stations, in an arid region, northwest China. It was found that the precision of the models was respectively improved when adding S_un, U₂, and R_h into the temperature-based model. Moreover, the importance evaluation of indices indicated PMF-56 ET₀ was more readily influenced by R_h with the exception of air temperature in this region. The best performance was achieved by the RF8 model with all the meteorological arguments as inputs. Although the precision of the model depends on the number of input climatic variables, all of the combinations of RF model turned out to be capable of producing reliable precision in ET₀ modeling, as mentioned above. Thus, the RF model should be the recommended model for PMF-56 ET₀ modeling in arid regions where weather data are limited. The Monte Carlo simulation technique was also employed for quantifying RF model uncertainty. The results of uncertainty analysis indicated that the PMF-56 ET₀ values bracketed by 95% confidence interval (95PPU) were larger, namely, most of the PMF-56 ET₀ values fell within 95PPU, and d-factor values calculated by upper and lower limits of the confidence interval were smaller in the studied area. In summary, the RF model is considered as an appropriate way of forecasting PMF-56 ET₀ and can provide an alternative tool under a minimal amount of climate data, as well as a reference for water resources management.

It should be noted that there have been extensive studies comparing the performance of RF and other artificial intelligence models for simulating ET₀, regarding the fact that almost all studies confirmed that RF model achieved higher simulation precision than others. Under these circumstances, this article did not conduct comparative research. In addition, although the RF model provides significant potential for more accurate estimation of the ET₀ with a lack of appropriate weather data in arid regions, certain drawbacks still persist. In this investigation, the selected sites are insufficient and the amounts of data size used to develop the model are smaller; besides, maximum and minimum ET₀ values simulated by RF model cannot accurately reflect the observed data. Therefore, further study can potentially focus on choosing more studied points with different climate types and combining RF method and another algorithm, such as Kalman filtering and wavelet transform techniques, for obtaining more reliable and practical results.

This research was funded by the National Key Research and Development Program of China (2017YFC0404305), Key Research Program of Frontier Sciences, CAS (QYZDJ-SSW-DQC031), Natural Science Foundation of Gansu province, China (18JR4RA002,18JR3RA393), CAS ‘Light of West China’ Program and the USQ-CAS collaborative research agreement USQ943692018 (2019–2021). The authors thank all the anonymous reviewers for their constructive comments.

The Supplementary Material for this paper is available online at https://dx.doi.org/10.2166/nh.2020.012.