Abstract
Regions with a large percentage of precipitation occurring near freezing experience high percentages (>10%) of misclassified precipitation events (rain versus snow) and necessitate efforts to improve precipitation phase determination schemes through the use of more accurate surface air temperature thresholds (Trs). Meteorological data from 169 sites in Scandinavia were used to test the applicability of using physiographic categories to determine Trs. Three classification methods involving varying degrees of automation were evaluated. The two automated methods tested did not perform as well as when tested on a smaller region, showing only 0.16% and 0.20% reduction in error. A semi-manual method produced the largest average reduction in misclassified precipitation (0.53%) across all sites. Further refinement of classification criteria for mountain and hill stations showed that at mesoscales (>5 km), maximum elevation is a better predictor of Trs (0.89% average reduction in error) than terrain relief (0.22%), but that relief becomes increasingly important at microscales (0.90%). A new method for categorizing mountainous stations based on upslope or downslope air movement increased the average reduction in error up to 0.53%. These results provide a framework for future landscape classification methods and confirm the importance of microscale topography for determining Trs in alpine regions.
INTRODUCTION
The precipitation phase determination scheme (PPDS) is one of the most important parameters in a snow model (Kongoli & Bland 2000), yet remains one of the most difficult tasks for hydrologists and meteorologists (Lackmann et al. 2002) in temperatures near freezing. PPDS in conceptual hydrological models often uses a set surface air temperature (Ta) threshold (Trs), assigning all precipitation events with Ta > Trs as rain and all other events as snow. This is a simplistic approach (Daly et al. 2000; Harpold et al. 2017a) that does not account for the influence of atmosphere or landscape variables such as topography (e.g. Harpold et al. 2017b), warm and cold air-mass boundaries (Feiccabrino et al. 2012), or ocean temperatures (e.g. Dai 2008). As a result, mid- to high-latitude, topographically complex regions, such as Scandinavia, are associated with high rates (10–40%) of misclassified precipitation events occurring between −3 and 5 °C. The previous work has shown that deriving Trs from groups of physiographically similar sites reduces misclassified precipitation and holds promise as a low-cost method for improving PPDS (Feiccabrino & Grigg 2016). This study builds upon this work by (1) assessing whether this approach can effectively be applied to a wider latitudinal range of meteorological stations and (2) examining the relative importance of elevation and terrain relief at different spatial scales in assigning Trs in mountainous regions.
Correct classification of precipitation into solid and liquid phases is of paramount importance due to vastly different atmospheric, hydrological, and ecological responses to rain and snow (Ye et al. 2013). Near freezing precipitation events can be modeled differently depending on the snow fraction (SF) assigned by a PPDS, which impacts snowpack properties such as snow density, albedo, snowpack layering, and water retention capacity (Loth et al. 1993). Precipitation mass corrected for gauge undercatch of 2–14% rain and 5–80% snow (Kokkonen et al. 2006) is also affected by changes in modeled SF. PPDSs can be based on hydrological, meteorological, or combined approaches. This paper focuses on conceptual hydrologic models, typically using surface temperature and precipitation mass as the meteorological forcing inputs and widely used across many disciplines for their simplicity, data availability, and low computational requirements.
The most commonly used surface temperature measurements for Trs are air temperature (Ta), dew-point temperature (Td), wet-bulb temperature (Tw), or a combination of relative humidity (RH) and Ta (Ye et al. 2013). However, Td, Tw, and RH measurements are much less available than Ta, and Trs is typically calibrated to the air temperature resulting in the least misclassified precipitation. This calibration is often conducted over large areas irrespective of terrain, ocean, or seasonal influences (Jennings et al. 2018). Trs are known to vary with time and location; however, the practice of validation and calibration of Trs adjusted over different land surfaces is not often applied in models or research (Harpold et al. 2017a). The use of broadly established set Trs indirectly forces two notable assumptions: (1) that near surface air is coupled to the atmosphere above, without substantial differences due to physiographic or biophysical processes on the Earth's surface (Aalto et al. 2018) and (2) that atmospheric conditions and energy exchanges from precipitation microphysics are invariant. These assumptions are incorrect when the atmospheric lapse rate (the rate of air temperature decrease with height) is greater (more unstable) than normal which has been shown to occur over open ocean water (Dai 2008) due to conductive heat transfer from the water to the lower atmosphere. In mountainous terrain, these assumptions are also often invalidated because of mechanical lifting and cooling of air at the dry adiabatic lapse rate (9.8 °C/km), which is more unstable than the average atmospheric lapse rate (6.5 °C/km) (Feiccabrino et al. 2015). These and other landscape-driven changes in environmental lapse rates represent an opportunity to improve high rates of misclassified precipitation (error) rates (>10%) associated with PPDS in near freezing temperatures.
Open-water conductively warms Ta at the ground–atmosphere interface, while the frozen, snow-covered ground has the opposite effect. In the winter, warming of surface air temperatures over ice-free water near Iceland has been shown by Ólafsson & Haraldsdóttir (2003) to increase the environmental lapse rate compared to a model assumed constant lapse rate (e.g. 7.5 °C in CHRM (Fang et al. 2013)). The higher atmospheric lapse rate over water gives snow a better chance of reaching the ground at a given Ta than over land, which results in a wide range of Trs in Iceland from 0.5 °C in inland areas to 2.1 °C along the northern coast (Ólafsson & Haraldsdóttir 2003). Dai (2008), using global 3-h data, found average land Trs 1.2 °C and ocean Trs 1.9 °C, a 0.7 °C warm bias over oceans. Similar findings at the regional scale were found in Scandinavia where ocean stations were shown to have a warmer Trs than land stations by ≈ 0.5 °C (Feiccabrino & Grigg 2016).
Precipitation patterns are strongly affected by geographic barriers (hills/mountains) at multiple geographic scales. Air forced to rise over geographic barriers can cool to saturation, allowing water vapor to condense and cause enhanced precipitation on the windward side (Roth et al. 2018). This orographically enhanced precipitation causes a thicker melting layer, a lower 0 °C isotherm, and a decrease in the snow elevation compared to upwind areas (Minder et al. 2011). Descending air on the lee side of mountains dries resulting in lower atmospheric RH, more sublimation, and decreased precipitation totals compared to the windward side (Jennings et al. 2018; Roth et al. 2018). Geographic barriers also cause large variance in precipitation mass at local scales (Henn et al. 2018) and increased PPDS errors when fixed lapse rates or Trs are used (e.g. Harpold et al. 2017b). The impacts of topography on Trs are evident in data from Northern Hemisphere land stations which show the coolest Trs in lowland and maritime climates and the warmest Trs in continental mountain climates with a maximum of 4.5 °C on the Tibetan Plateau (Jennings et al. 2018). At local to regional scales, Trs from Scandinavia show a similar trend with lowland stations in Scandinavia having cooler Trs (1.0–1.1 °C) than hill and mountain stations Trs (1.2–1.6 °C) (Feiccabrino & Grigg 2016).
This study aims to verify and improve upon the methods developed by Feiccabrino & Grigg (2016) with the end goal of decreasing precipitation phase uncertainty in conceptual hydrological models. The previous study developed a geographic information system (GIS)-based landscape classification method for a 7° latitude wide cross-section in Scandinavia from the North Sea over the Scandinavian mountains to the Bay of Bothnia and resulted in a reduction in error of 0.59% and 1.26%, using Ta and Tw-based Trs, respectively. The limited latitudinal range of the Feiccabrino & Grigg's (2016) study enabled the simplification of windward versus leeward stations based on the prevailing westerlies but restricted the future applicability of this method as a broader country-wide approach to assigning Trs. In this study, an expanded data set of 169 meteorological stations from across Norway and Sweden was used, and the simplified windward versus leeward designation was replaced with the GIS-derived classification of the upslope or downslope movement of air that can be modified based on the regional prevailing wind direction. This study also introduces a semi-manual classification scheme that considers the more varied wind sources and substrates of the larger geographic area.
Another aspect of the previous study's results which is addressed is the difficulty in assigning Trs for mountain and hill stations. The mountain and hill categories had both the largest Trs variability between stations and the highest percent of misclassified precipitation. This study tests which elevation variables, in addition to or in replace of relief, are most predictive of Trs at multiple spatial scales. The original study used a 15 km buffer surrounding each station in order to assign landscape categories, which may have been too coarse to capture microscale (<10 km) differences in topography. Harpold et al. (2017a) suggest that in mountain and coastal environments, extreme changes in atmospheric dynamics over short distances may drive microscale differences in optimal Trs.
STUDY AREA
A meteorological observation data set for the years 1995–2011 from 84 Norwegian stations and 85 Swedish stations was selected for the analysis (Figure 1(a)). This data set was obtained freely from the websites of the Swedish Meteorological and Hydrological Institute (SMHI 2018) and the Norwegian Meteorological Institute (NMI 2018). Data collected after 2011 were not used because of multiple changes to data collection methods made between 2011 and 2014. The Scandinavian peninsula is an ideal setting to develop new physiographic-based PPDS because of its latitude, proximity to the ocean on three sides, prevailing wind direction, and topographic diversity over a relatively small region.
Much of the weather impacting Scandinavia travels over the North Atlantic Ocean and the Norwegian Sea and is subject to an increase in the environmental lapse rate causing higher Trs over water. The prevailing westerly winds then experience rising/cooling followed by sinking/warming, as they travel over the Scandinavian mountains which run SSE–NNW along the Norway/Sweden border. These orographic effects change the moisture content and precipitation intensity both of which impact Trs. The stations range in latitude from ≈ 55°N in southern Sweden to 77°N in Svalbard and in elevation from the sea level to 723 m. 20% to 60% of total precipitation in Scandinavia occurs between −3 and 5 °C, and the average rate of misclassified precipitation at these temperatures is 11.00% (Figures 1(b) and 1(c)).
METHODS
Original-automated classification – all sites
The original-automated (GIS)-classification method (Table 1) by Feiccabrino & Grigg (2016) was tested, but this time on an expanded and more geographically complex set of meteorological stations. Basic physiographic categories were determined by user-defined criteria to classify sites as the ocean, coast, fjord, rolling, hill, or mountain (Table 1; Figure 1(a)). Key aspects of this method are the differentiation of ocean-influenced versus land sites within a 15 km radius of the station, the differentiation of land sites based on maximum relief within a 15 km station radius, and the assumption that prevailing westerlies caused all Norwegian stations to be windward and all Swedish stations to be leeward regardless of local terrain. Further details on the method can be found in Feiccabrino & Grigg (2016).
Physiographic category . | % Water within 15 km . | Elevation change within 15 km . |
---|---|---|
Ocean | 90–100 | N/A |
Coast | 60–90 | |
Fjord | 40–10 | |
Rolling | <10% | 0–499 m |
Hill | 500–999 m | |
Mountain | Above 1,000 m |
Physiographic category . | % Water within 15 km . | Elevation change within 15 km . |
---|---|---|
Ocean | 90–100 | N/A |
Coast | 60–90 | |
Fjord | 40–10 | |
Rolling | <10% | 0–499 m |
Hill | 500–999 m | |
Mountain | Above 1,000 m |
New-automated method – all sites
As in the previous study, ArcGIS was used to develop a ‘new-automated’ method that attempted to improve on the simplistic country-based classification of windward and leeward sites. This approach uses a measured upslope or downslope classification based on the location of the maximum elevation within each 15 km station buffer relative to the station (Supplemental Table S-1 in Supplementary Materials). The orientation of the maximum elevation point relative to the station was determined using the Euclidean Direction tool within the Spatial Analyst extension. Sites with a maximum elevation located to the west of the station were classified as downslope sites, while sites with a maximum elevation located to the east were classified as upslope. Elevation data were obtained from a 30-arc second European digital elevation model (GTOP30; USGS 1996). This upslope or downslope distinction was then combined with the 15 km maximum relief-based categorization of land sites from the original-automated method (Table 1; Figure 1(a)).
Semi-manual method – all sites
With the expansion of the study area (Figure 1), the number of impacting climatic variables increased. Sites from southern Scandinavia could have weather coming from Europe over the Baltic Sea which would not experience the same orographic effect as sites further north, while Svalbard stations are much further north, have frozen ground and little relief. To address this variability, a semi-manual classification (Supplemental Table S-1) was done using the original-automated method to first determine basic physiographic categories (Table 1). Further categorizations of the land sites were then determined manually based on the relative location of the station to mountains or high-ground and the potential for different air mass source region effects on advecting air. Ocean, coastal, and fjord stations were similarly classified based on the direction of the expected mean wind flow and substrate (i.e. frozen ground, deep or shallow ocean; see Supplemental Table S-2).
Relief versus elevation – mountain and hill sites only
In Scandinavia, relief of the terrain is believed to have a stronger effect than the elevation on changing atmospheric conditions due to the windward ascent of air and the leeward decent of air during precipitation events. However, since the elevation is readily available and often used in hydrological models (Lehning et al. 2011; Henn et al. 2018; Roth et al. 2018), this study tests the use of station, maximum and average elevation from within a 15 km station radius as classification criteria (Supplemental Table S-3). These trials were conducted only on the 37 stations originally classified as mountain or hill sites, where relief was large enough to potentially impact Trs. The methods followed the new-automated classification scheme except maximum and average elevation within the 15 km radius or station elevation replaced maximum relief.
Station radius size – mountain and hill sites only
A final variation in the classification methodology used 5 and 1 km station radii to determine the different GIS-derived elevation parameters for mountain and hill sites (Table 2) in order to better characterize site-specific microscale topography. This approach used the new-automated methodology at 5 and 1 km radii to determine relief and upslope versus downslope. Further tests using 5 and 1 km radii were performed for the different elevation parameters described in the previous section.
Station radius size . | Relief or elevation parameter . | % Reduction in misclassified precipitation . |
---|---|---|
15 | Maximum relief | 0.22 |
5 | 0.64 | |
1 | 0.90 | |
15 | Maximum elevation | 0.89 |
5 | 0.72 | |
1 | 0.89 | |
15 | Average elevation | 0.36 |
5 | 0.34 | |
1 | 1.09 |
Station radius size . | Relief or elevation parameter . | % Reduction in misclassified precipitation . |
---|---|---|
15 | Maximum relief | 0.22 |
5 | 0.64 | |
1 | 0.90 | |
15 | Maximum elevation | 0.89 |
5 | 0.72 | |
1 | 0.89 | |
15 | Average elevation | 0.36 |
5 | 0.34 | |
1 | 1.09 |
Bold values highlight the radius with the largest percent reduction in error for each elevation parameter.
Calculations
Using the same method as Feiccabrino & Grigg (2016), non-precipitation observations for each station were removed along with precipitation events occurring in air temperatures warmer than 5 °C and cooler than −3 °C. Mixed phase observations were removed, as were freezing rain/drizzle observations, which can be characterized as rain or snow and consisted of less than 1% of the remaining precipitation observations. The total sum of misclassified precipitation events for each 0.1 °C interval of possible Trs between −3 and 5 °C was calculated as the total of PPDS assigned snow events observed as rain and PPDS assigned rain events observed as snow divided by the total number of precipitation events in the observation data set. See Feiccabrino & Grigg's (2016) study for further details on calculating misclassified precipitation.
Optimum Trs values assigned for physiographic groups in a PPDS were set at the air temperature with the lowest group average of station percent misclassified precipitation. This is not necessarily the Trs that would result in the lowest number of misclassified events, as the sample size was different for each station. However, this allowed stations with small and large sample sizes to be given equal weighting toward the Trs for the physiographic groups within a PPDS.
RESULTS
All sites – comparison of three methods
Optimal Trs assigned by all three landscape-based PPDS classification methods resulted in slight reductions in the average percent of misclassified precipitation relative to the country-wide thresholds (Figure 2; Supplemental Table S-1). The semi-manual classification method showed the lowest average percent of misclassified precipitation at 10.49% which amounted to a 0.53% average reduction in error compared to the single country-wide threshold. The results from the two automated methods were very similar with the new-automated method performing slightly better than the original-automated method, 0.20% versus 0.16% average reduction in misclassified precipitation, respectively. In all three methods, the largest average reduction in error was seen in Norway, with the semi-manual method showing the largest difference between countries, 0.88% for Norway, and 0.19% for Sweden.
An examination of the spatial distribution of the results from all three methods grouped by the major physiographic category shows that the average percent reduction in error was greatest for sites categorized as the ocean (0.68%) and mountain (0.67%). However, these results also show that mountain and hill sites have the highest average percent of misclassified precipitation, at 12.02% and 11.81%, respectively. The rolling and coastal categories showed the least improvement at 0.17% and 0.11% but along with the fjord sites have the lowest percentages of misclassified precipitation (9.60–10.73%).
Mountain and hill sites – relief versus elevation
The higher rates of misclassified precipitation events from all three methods for mountain and hill categories were the motivation to further refine the classification of topographically complex sites. The first set of calculations compared the use of maximum relief versus various elevation parameters from within a 15 km station radius (Supplemental Table S-3). When no additional classification of sites based on the prevailing wind direction was used, station elevation produced the greatest (0.75%) reduction in error of misclassified precipitation, while average elevation produced the smallest reduction in error (0.19%). However, when the categories were expanded to include upslope and downslope wind designations, relief and all elevation parameters showed some increase in the reduction of error (Figure 3). The use of maximum elevation within the 15 km station combined with upslope versus downslope yielded the largest reduction in error at 0.89%. All the elevation parameters outperformed relief when using the 15 km station radius, with the exception of average elevation alone (no prevailing wind direction) (Figure 3, Supplemental Table S-3).
Mountain and hill sites – station radius size
The comparison of using 15, 5, or 1 km station radii to determine maximum relief showed a steady increase in the reduction in error of misclassified precipitation with decreasing station radii (Figure 3, Table 2). The use of maximum relief within a 1 km radius resulted in the largest reduction in error at 0.90%. When the smaller station radii were used to determine maximum and average elevation, the results showed no improvement between 15 and 5 km radii. However, the use of a 1 km radius improved these results, with the largest reduction in misclassified precipitation (1.09%) occurring with average elevation calculated for a 1 km radius.
DISCUSSION
Automated versus semi-manual methods
The three basic landscape classification methods (original-automated, new-automated, and semi-manual) tested in this study all resulted in an average percent reduction in misclassified precipitation events relative to the use of a single country-wide threshold (Supplemental Table S-1, Figure 2). However, both automated methods resulted in less improvement (0.16% and 0.20% versus 0.59% reduction) than the Feiccabrino & Grigg (2016) study, which used the original-automated method on a smaller, geographically constrained set of stations. The difference in the performance of the same method between this and the previous study can be explained by the larger and more climatically heterogeneous set of stations used in the current study. The semi-manual method addresses the expected shortcomings of the larger data set and results in a similar reduction in error (0.53%) as the original study (0.59%). One of the key variables accounted for by the semi-manual method and not the automated methods is the changing precipitation source area from north to south along the Scandinavian peninsula (Supplemental Table S-2). For example, coastal stations in northern Norway often receive snow advected from the colder Barents Sea which is less likely to melt while falling to the ground, regardless of surface temperatures, and results in a higher Trs (1.5–1.4 °C) than those for central Norway coastal sites (1.2–1.3 °C), which receive weather from the Gulf Stream-warmed Norwegian Sea. Likewise, southern coastal sites have the lowest Trs (0.6 °C) because of the impact of a modified continental air mass from the European mainland with less maritime influence. A future approach to objectively account for diverse precipitation source areas is to use wind and storm track data to first delineate regions with similar precipitation source areas that can then be classified separately by other landscape variables.
Although the new-automated method with GIS-derived upslope or downslope designations showed only a slight improvement relative to the original-automated method (0.16% versus 0.20%), these results would likely be improved when used on a smaller subset of stations or modified to accommodate known geographic changes in the prevailing winds. This method was most effective relative to the other methods at reducing average % error for the hill stations (0.5%) and was least effective (0.1%) at Norway mountain stations (Figure 2). These results suggest that small and previously unrecognized changes in upslope and downslope conditions in regions of low-to-moderate relief may impact Trs, but that in areas with multiple mountain ridges and passes such as the southern Norwegian mountains, other physiographic factors are important. In intermountain areas, the air loses moisture after each successive mountain pass, and the atmospheric changes caused by upslope and downslope winds become less effective than the first coastal mountains. Additionally, dense cold air can get trapped in the valleys between mountains causing warm air to pass over the mountains without experiencing orographic lifting.
PPDSs have broad societal applications across multiple disciplines including weather forecasting, transportation safety, and snowpack modeling. Conceptual PPDSs rely on a best estimate of Trs that can be applied across both spatially and temporally changing atmospheric conditions. The methods examined in this study provide the basis for an automated approach to deriving Trs that accounts for spatial, landscape-driven atmospheric changes averaged through time using widely available meteorological data that are a better starting point for the development of conceptual PPDS than broadly standardized Trs.
Refinement of classification for topographically complex regions
Higher rates (>11%) of misclassified precipitation at mountain stations across all three methods were the motivation for the further refinement of landscape classification for mountain and hill stations. Results show that the use of elevation parameters (station elevation, maximum elevation, and average elevation) instead of relief as a classification criteria yielded better results when using the 15 km station radius (Figure 3). The largest reduction (average: 0.89%) occurred using maximum elevation and upslope/downslope wind direction (Supplemental Table S-3). These results indicate that atmospheric conditions within 15 km are impacted by the up- or downslope movement of air caused by topographic highs. Additionally, Jennings et al. (2018) suggest that at higher elevations, less dense air may allow snow to persist at higher temperatures, and thus local high points could warm Trs and may also explain the better performance of maximum elevation.
A reduction in the station radii from 15 km to 5 km and to 1 km improved the performance of relief and average elevation as criteria for classification (Figure 3; Table 2). Maximum elevation did not improve between the 15 and 1 km radii, which is consistent with the explanation of lower air pressure at higher elevation impacting Trs (e.g. Jennings et al. 2018) because altitudinal changes in air pressure are laterally consistent. The increased importance of relief and average elevation at microscales suggests that orographically enhanced precipitation and lowering of both the zero-degree isotherm and snow lines (Minder et al. 2011) on the windward side of a mountain begin to affect Trs within 5 km of the station.
The findings from the mountain and hill data set contribute to other efforts to better predict snowfall in topographically complex regions (e.g. Marks et al. 2013). This study confirms the importance of microscale landscape parameters and could be combined with the microscale snow hydrology work by Cristea et al. (2017) which uses the relative elevation position and aspect to better characterize the snow cover in alpine regions. Future efforts of landscape classification should include the use of maximum elevation at macroscales (>10 km) to develop a first order of classification. Mountain and hill sites (maximum elevations >499 m) can then be further classified using a microscale (<5 km) analysis of terrain relief and/or average elevation. In both stages of classification, the interaction between topography and the prevailing wind direction should be accounted for by an analysis of upslope versus downslope conditions.
CONCLUSIONS
This study contributes to an expanding body of work on the use of physiographic variables to derive more accurate Trs at global, hemispheric, and regional scales (Ólafsson & Haraldsdóttir 2003; Dai 2008; Harpold et al. 2017a; Jennings et al. 2018). These efforts are crucial to improving conceptual PPDS, which are widely used across meteorology, hydrology, and ecology. The methods tested reveal both the potential for and limitation of automated, GIS-based landscape classification at regional scales. The original classification method when applied to a climatically diverse region did not perform as well as it did when applied to the smaller latitudinal range used by Feiccabrino & Grigg (2016), although it still resulted in a small improvement (0.16% reduction in error) in misclassified precipitation relative to the set country-wide Trs. The modification of the original method to include the GIS-derived upslope versus downslope air movement showed a 0.20% reduction in error and also fell short of improvements made in the original study. Another modification tested was the manual sub-classification of stations following the initial automated landscape classification which resulted in a 0.53% average reduction in error and highlights the need to better integrate precipitation source area and intermountain physiographic effects into future automated methods.
A closer examination of mountain and hill stations indicates that with the 15 km station buffer, the greatest average reduction in misclassified precipitation (0.89%) occurred when using maximum elevation and the upslope and downslope movement of air. When smaller station buffers were used, maximum elevation results remained unchanged, but the use of terrain relief and average elevation improved average reduction in error rates up to 0.73%.
Based on the results of this study, several recommendations can be made for future work on landscape classification. The first is that the use of GIS-based, automated methods in climatically diverse regions could be improved by an initial classification of the precipitation source area based on the analysis of wind direction data during precipitation events. Once similar precipitation source regions are established, stations can be further classified using a 15 km station buffer to determine: (1) proximity to the ocean, (2) maximum elevation, (3) and whether the station is up- or downslope relative to the precipitation source. For stations with maximum elevations >499 m, a third tier of classification that uses terrain relief and/or average elevation within a 1 km station buffer is recommended. This and future work aimed at improving the landscape classification of Trs provide needed and practical approaches to decreasing the error currently associated with conceptual PPDS in mid- and high-latitude regions.
SUPPLEMENTARY MATERIAL
The Supplementary Material for this paper is available online at https://dx.doi.org/10.2166/nh.2020.081.