Effects of mountain pine beetle-killed forests on source water contributions to streamflow in headwater streams of the Colorado Rocky Mountains

Abstract

Natural or human-influenced disturbances are important to the health and diversity of forests, which in turn, are important to the water quantity and quality exported from a catchment. However, human-induced disturbances (prescribed fire and harvesting) have been decreasing, and natural disturbances (fires and insects) have been increasing in frequency and severity. One such natural disturbance is the mountain pine beetle (MPB), (Dendroctonus ponderosae) an endemic species. A recent epidemic resulted in the mortality of millions of hectares of lodgepole pine (Pinus contorta) forests in Colorado, USA. Beetle-induced tree mortality brings about changes to the hydrologic cycle, including decreased transpiration and interception with the loss of canopy cover. This study examined the effect of the mountain pine beetle kill on source water contributions to streamflow in snowmelt-dominated headwater catchments using stable isotopes (²H and ¹⁸O) as tracers. Study catchments with varying level of beetle-killed forest area (6% to 97%) were sampled for groundwater, surface water, and precipitation. Streams were sampled to assess whether beetle-killed forests have altered source water contributions to streamflow. Groundwater contributions increased with increasing beetle-killed forest area (p=0.008). Both rain and snow contributions were negatively correlated with beetle-killed forest area (p=0.035 and p=0.011, respectively). As the beetle-killed forest area increases, so does fractional groundwater contribution to streamflow.

Graphical abstract

Keywords

mountain pine beetle / isotope tracers / streamflow generation / headwaters

1 Introduction

The extreme precipitation process of landfalling tropical cyclones (LTCs) often severely threatens life and property in the coastal area of the Northwest Pacific Ocean, and the accompanied very heavy rainfall (>100 mm/24 h) is the main cause of flood disasters (Chen et al., 2019; Yu and Chen, 2019). China experiences the most LTCs in the world, with an average of 7–8 such events per year. Although LTCs continues to weaken after landfall, they often remain as very strong rainstorm systems, and can still bring extreme rainfall of more than 1000 mm in 24 h to China (Chen and Xu, 2017).

The precipitation process of LTCs is very complex. Previous studies have shown that in addition to the important influence of large-scale environmental field on LTC precipitation (Wang and Holland, 1996; Rogers et al., 2003; Lonfat et al., 2004; Chen et al., 2006; Lonfat et al., 2007; Wingo and Cecil, 2010; Jiang and Ramirez 2013; Reasor et al., 2013; Yu et al., 2015), complex underlying surface forcing (Yu et al., 2010, 2017) such as water, terrain and mesoscale convective activities all play a significant role in LTC precipitation (Yu et al., 2010). Wang et al. (2012) showed that the Threat Score (TS) value of 25 mm-rainfall of 24 h forecasts in LTCs over China was only 0.2. Although the TC track forecast has made great progress in the past decades, the TC precipitation forecast ability is still far behind the track forecast (Yu et al., 2020). Hence, accurate typhoon rainfall forecast is still full of challenges.

To improve the typhoon rainfall forecast ability, it is necessary to know the LTC precipitation forecast performance of numerical models in a targeted specific way. At present, some traditional verification methods such as Equitable Threat Score (ETS), Probability of Detection (POD) and False Alarm Ratio (FAR), have been widely used to evaluate quantitative precipitation forecast (QPF). However, although these methods have some application value, they cannot provide the detailed error information to their users.

In recent years, a “non-traditional” object-oriented verification method has been developed. The “precipitation field” is a continuous area where the precipitation rate or accumulated precipitation exceeds a certain threshold value defined as the “object” of the verification. A few object-oriented methods have been developed in the past 20 years, including CRA (contiguous rain area, Ebert and McBride, 2000), MODE (Method for Object-based Diagnostic Evaluation, Davis et al., 2006, 2009), and SAL (Structure – Amplify – Location, Wernli et al., 2008). These methods have been further applied to ensemble precipitation forecasts, mesoscale convective systems, and other evaluation community studies (Ebert and Gallus, 2009; Gilleland et al., 2009; Gallus, 2010; Demaria et al., 2011; Moise and Delage, 2011; Zacharov et al., 2013; Clark et al., 2014). As early as 2007, Marchok et al. (2007) tried to analyze the LTC precipitation characteristics through precipitation patterns, averages and extremes in the United States. Since then, the researches on QPF evaluation by object-based method have been less developed, especially for LTCs in China.

As one of the earliest proposed object-oriented verification methods, CRA has been widely used in different types of rainfall verification such as mesoscale convective rainfall, seasonal precipitation, monsoonal rain and ensemble precipitation forecasts (Gallus 2010; Demaria et al., 2011; Moise and Delage, 2011; Dube et al., 2014; Ashrit et al., 2015; Murphy et al., 2015; Sharma et al., 2015, 2017). However, the application of CRA verification to TC rainfall forecasts is still very limited (Ebert et al., 2005, 2011; Chen et al., 2018; Yu et al., 2020). The original CRA method proposed by Ebert and McBride (2000) mainly considered three error components: displacement, volume and pattern. Moise and Delage (2011) modified the CRA method by considering the rotation of the whole rainfall structure. This technique has demonstrated its usefulness in evaluation of the climate model metrics. Chen et al. (2018) evaluated TC rainfall forecasts over the ocean using a modified CRA method. Rotation of TC rainfall distributions was introduced as a new error source. In this new version of CRA method, the total TC rainfall error can be decomposed into four components, which are respectively from the rotation, displacement, volume and pattern.

Yu et al. (2020) adopted the CRA verification method to discuss the forecast ability of the ACCESS_TC (the TC version of the Australian Community Climate and Earth System Simulator) regional model of the Australian Meteorological Bureau for the precipitation of the landfalling typhoons in China. It was the first systematic study based on years of numerical forecasts to carry out QPF verification of LTCs by using object-oriented verification method of CRA. Thus it could provide an example to find the error sources and to compare the performance of LTC rainfall forecasts for different weather numerical models. The comparison among models including both regional models and global models is especially necessary. It helps to continually record the progress and error sources of QPFs in those numerical models, so that different ways to improve LTC rainfall forecasts can be explored.

Although the above TC rainfall forecast verification studies have explored the forecast ability for different lead-time and rainfall thresholds, little attention has been paid to the rainfall forecast ability during a TC life cycle. It is hard to analyze the physical error sources only from the lead-time dependent forecast abilities, especially for TCs with long durations. One lead forecast time generally covers different life stages of a TC, so that the error analysis based on lead time would combine the effects of many factors. Therefore, the investigation from the view of TC life cycle can help to reveal the major physical factors in different life stages. Lekima (2019) is a good LTC case for the study in this paper. It made landfall as a super typhoon, maintained for several days over land and produced very heavy rainfall during moving inland. Since the rainfall process in a landfalling typhoon is complicated, which is influenced by the changing environments and underlying surfaces (Yu and Wang, 2018), the quantitative analysis of rainfall forecast abilities at different life stages of a LTC can illustrate both the variable rainfall forecast ability during TC landfall process and the related causes of the different forecast uncertainties. We believe that the findings of this study can be applied in analogous LTCs characterized by super landfall intensity, long duration and deep northerly inland movement.

This paper will focus on the LTC Lekima (2019) rainfall forecast errors in different landfall stages by using the CRA verification method. The rainfall forecast of an operational regional TC numerical model, GRAPES-TCM (the TC version Model of Global and Regional Assimilation PrEdiction System) developed by Shanghai Typhoon Institute of China Meteorological Administration (STI/CMA) is used in this study. The remainder of this paper is organized as follows. Section 2 introduces the data and the analysis method. Section 3 gives a general introduction of LTC Lekima (2019). The verification results and the diagnostic results related to the rainfall forecast errors are presented in Section 4 and Section 5, respectively. Conclusions and discussions are given in Section 6.

2 Data and analysis methods

2.1 GRAPES-TCM model

Forecasts from the regional TC numerical model of GRAPES-TCM are used in this study. The GRAPES-TCM has been adapted and updated for operational forecasting applications on TCs for several years (Huang et al., 2007; Tan et al., 2021). The vertical direction adopts the height-based terrain-following coordinate (Klemp coordinate), and it is divided into 51 layers. The top model layer is set as 35 km. The physical parameterization schemes in the model employ the package of the NCEP Global Forecast System (GFS). The model domain is nested with a finest resolution of 0.1°. It makes 72 h forecasts with a new run every 6 h. The TC vortex initialization is based on the estimates of TC location, intensity and storm size. The initial and lateral information comes from the global forecast model of GFS. To evaluate the rainfall forecast ability of GRAPES-TCM during the period of typhoon Lekima affecting the coastal areas of China, the model prediction data from 24 h before landfall to the end of its life after landfall are selected. We will focus on the rainfall forecast ability of GRAPES-TCM during different life stages of the LTC Lekima (2019). The rainfall thresholds of 30/50/100/250 mm are chosen in the CRA verification.

2.2 TC best track data

Intensities and locations of Typhoon Lekima are extracted from the best track data officially released by the STI/CMA.The TC position data at landfall time is obtained by data interpolation from the nearest analysis time.

2.3 Gauge-satellite merged rain data

Gauge-satellite merged rain data over China during recent years has been produced based on the CMORPH (Climate Precipitation Center Morphing) satellite rain product by the National Meteorological Information Center of the CMA (Shen et al., 2010, 2014), which is the reference data in this study. It is known that satellite rain data, including the CMORPH precipitation product, could present reasonable rainfall distributions but tend to overestimate light rainfall and underestimate heavy rainfall (Yu et al., 2009; Chen et al., 2013). The gauge-satellite merged rain data can improve the accuracy by using the rain gauge observations over China (Shen et al., 2010). The merged rain product features improved quality, smaller bias over land, and finer spatio-temporal resolutions (0.1° and 1 h) (Shen et al., 2014). Severe weather systems could be better captured with hourly rainfall description. Therefore, gauge-satellite merged rain data could be used as a suitable reference to evaluate the GRAPES-TCM model rainfall forecasts during the landfall of Lekima.

2.4 Initial analysis fields of ECMWF forecast products

The initial analysis fields of European Centre for Medium-Range Weather Forecasts (ECMWF) are used in this paper to compare with the GRAPES-TCM forecasts, which could help find the possible forecast errors of some related physical processes. The horizontal resolution of the high-level field is 0.25°. The physical quantities to be examined include environmental VWS, relative vorticity, frontogenesis, and advection of equivalent potential temperature.

2.5 CRA verification method

Similar to Yu et al. (2020), the CRA verification method is used in this study. The main difference between the CRA and other object-based verification methods, including the MODE (Davis et al., 2006, 2009) and the SAL (Wernli et al., 2008), is that the CRA decomposes the error and identifies its major contributors. Particularly, important attributes such as the first (or secondary) error source could be obtained by both numerical model developers and other users.

The errors of the forecasted rain field are decomposed into displacement (D), rotation (R), volume (V) and pattern (P) error components, which can be written as follows,

(1)${\mathrm{MSE}}_{\mathrm{total}}={\mathrm{MSE}}_{\mathrm{displacement}}+{\mathrm{MSE}}_{\mathrm{rotation}}+{\mathrm{MSE}}_{\mathrm{volume}}+{\mathrm{MSE}}_{\mathrm{pattern}}.$

Based on the minimum squared error best-fit criterion, the error components can be calculated. Finally the forecasted rain field can be adjusted via shifting and rotating to best fit the observed rainfall. By moving (M) displacement and adjusting rotation (R), the forecasted rain field can be corrected. More detailed information is available in Chen et al. (2018) and Yu et al. (2020).

2.6 Traditional verification methods

Some traditional verification measures including the ETS, POD and FAR, are also used in this study. Their definitions are shown as below.

(2)$\mathrm{E}\mathrm{T}\mathrm{S}={\displaystyle \frac{H-CH}{H+F+M-CH}},$

(3)$\mathrm{FAR}=F/(H+F),$

(4)$\mathrm{POD}=H/(H+M),$

where $CH=(H+M)\times (H+F)/(Z+H+M+F)$, and some thresholds (30, 50, 100, and 250 mm) are used to define the rainfall events with different intensities. Then it is scored as falling under one of the four categories of correct no-rain estimate, false alarms, misses, or hits (Z, F, M, or H). More detailed information could be found in Yu et al. (2020).

3 Overview of TC Lekima

The TC Lekima (2019) formed over the north-western Pacific (Fig. 1(a)) at 0600 UTC (Coordinated Universal Time) on August 4, 2019, and then slowly moved north-westward, with its intensity gradually increasing (Fig. 1(b)). After strengthening to the grade of typhoon (the maximum wind speed (MWS) = 32.7 m·s⁻¹) at 2100 UTC on August 6, its motion speed gradually increased. After a rapid intensification (RI), it became a super typhoon (MWS= 51 m·s⁻¹) at 1500 UTC on August 7, which maintained for 60 h. Lekima reached its maximum intensity of wind-force scale 17 (MWS= 62 m·s⁻¹) at 1800 UTC on August 8.

Fig.1 (a) The best track of typhoon Lekima and accumulated rainfall from 1200 UTC on August 8 to 0000 UTC on August 13, 2019 (unit: mm; shaded); (b) MSLP (unit: hPa; observation: solid black line, forecasts with different initial times: black dashed lines) in the typhoon center at 6-h intervals from 1200 UTC on August 6 to 0000 UTC on August 13, 2019. The vertical dashed line indicates the landfall time of Lekima. The dashed curves represent all precipitation forecasts in this study. The red one represents the case used for later detailed analysis (Fig. 7 and Fig. 8).

Full size|PPT slide

At 1745 UTC on August 9, Lekima made landfall in Wenling City of Zhejiang Province of China, with a MWS of 52 m·s⁻¹ near its center. Afterwards, it experienced a rapid decay (RD) process, decreasing to 35 m·s⁻¹ within 3 h, and weakened into a severe tropical storm (MWS= 30 m·s⁻¹) at 0100 UTC on August 10. At approximately 1300 UTC on August 10, Lekima passed through the northern region of Zhejiang Province and entered Jiangsu Province, with its intensity further weakening to a tropical storm (MWS= 23 m·s⁻¹), then moved northward, and entered the Yellow Sea at around 0400 UTC on August 11. At 1200 UTC on August 11 Lekima made its second landfall at Qingdao City of Shandong Province. Afterwards, it entered the Shandong Peninsula and meandered in the Laizhou Bay. At 0000 UTC on August 13, it became a tropical depression.

Therefore, the TC Lekima experienced both RI and RD processes in its life cycle. It was still a super typhoon when it made landfall, and became the third most intense landfall typhoon in Zhejiang according to China meteorological records. After landfall, it still moved northward (Fig. 1(a)) and lingered for up to 3 days (Fig. 1(b)), bringing heavy rainfall on coastal provinces in East China. According to the distribution of accumulated rainfall from 1200 UTC on August 8 to 0000 UTC on August 13 (Fig. 1(a)), there are two major centers of rainfall above 300 mm. One is located in the central and northern Zhejiang near the typhoon landing point, and the other is located in the central and northern areas of Shandong.

4 Model forecast verification

4.1 Rainfall forecasts with different lead times

Figure 2 shows conventional verification results of the Lekima rainfall forecasts, including ETS, POD and FAR (Yu et al., 2020), which are used to assess the model forecast performance (Figs. 2(a), 2(c) and (e)) for five different forecast lead times (0–6, 0–24, 24–48, 48–72, and 0–72 h) and four different rainfall thresholds (30, 50, 100, and 250 mm). The verification results of rainfall forecasts after the displacement adjustment of the CRA method are also illustrated (Figs. 2(b), 2(d) and 2(f)). For the 30 mm rainfall forecasts, the average ETS is not high (only 0.2) for short lead-time (6 h) forecasts (Fig. 2(a)), while that of the 24 h rainfall forecasts is the highest (0.4). This may be related to the model initialization. A certain dynamic imbalance in the physical quantity fields such as typhoon intensity, structure and environmental flow field may lead to a relatively long model “spin-up” process, and consequently could affect the 6-h rainfall forecasts. However, the model dynamic balance process has completed in 24 h. At this time, the model has relatively small forecast errors for the track and intensity of the typhoon, which also makes the score of 24 h rainfall-forecast the highest. With the increase of forecast lead time, the ETS of 48/72-h forecast rainfall gradually decreases to about 0.2.

Fig.2 (a), (c), and (e) The verification measures of ETS, POD and FAR for rainfall forecasts with different lead times. (b), (d), and (f) The same as (a), (c), and (e), but for the rainfall after the CRA displacement adjustment. The x-axis represents different forecast lead times. The solid, dashed, dotted lines, and the asterisks represent the results for 30 mm, 50 mm, 100 mm, and 250 mm rainfall, respectively.

Full size|PPT slide

Besides, the 0–72 h ETS also reaches 0.4, which is comparable to that of the 24 h forecast, indicating that the model still has a higher ETS for the accumulated precipitation in a longer period. The verification results above are consistent with those of the ACCESS_TC model in Yu et al. (2020), indicating that for the current regional numerical models, the 0–6 h 30 mm rainfall nowcasting still presents challenges. As for 50 mm rainfall forecasts, the variation of the ETS is similar to that of 30 mm forecasts, but with slight decreases. When the 24 h rainfall thresholds reach 100 mm and 250 mm, the ETS decreases significantly to 0.14 and 0.05, respectively.

To further analyze the sources of rainfall forecast errors in the GRAPES-TCM during the period influenced by the Typhoon Lekima, as introduced in Section 2.5, the mean square error (MSE) of rainfall forecasts in the verification area is decomposed by the CRA method into four error components: rotation (R), displacement (D), volume (V) and pattern (P). Their relative contributions to the total rainfall errors under different forecast lead times and rainfall thresholds are compared, and then the major error sources are determined. It is found that for all forecast lead times and rainfall thresholds (Fig. 3), the rainfall forecast errors of the model are generally dominated by P and D. The former is related to the distribution structure, while the latter is related to the TC track. Comparatively, the components V and R contribute less, and especially for R, the contribution is nearly negligible.

Fig.3 Average contribution rates of the four error sources to the total rainfall error for different forecast lead times and different rainfall thresholds. (a) 0–6 h, (b) 0–24 h, (c) 24–48 h, (d) 48–72 h and (e) 0–72 h. R, D, V, and P represent the rotation, displacement, volume and pattern error components, respectively. The x-axis is the rainfall threshold (mm) and the y-axis is the error contribution (%).

Full size|PPT slide

Specifically, for 6 h rainfall forecasts (Fig. 3(a)), the main errors come from rainfall patterns (~50%–60%). The D maintains 20%–30% for different rainfall thresholds, and the V also accounts for a certain proportion (~10%–20%, especially for 50 mm rainfall forecasts). For 24 h rainfall forecasts, the proportion of D increases with the increase of rain amount, and exceeds 50% for rainfall amount of 250 mm. The D of 48 h forecasts occupies a higher proportion (even close to 80%). In contrast with the other lead times except the 6 h lead time, the 72 h accumulated rainfall has less proportion of D and more proportion of P. This shows that the rainfall pattern error is the major source of error in this case. For the 250 mm rainfall, the proportion of rainfall volume error increases significantly, which shows that in terms of heavy rainfall, the model rainfall-extreme forecast also has obvious deficiencies.

Although the improvement in ETS for 6–24 h forecasts is not obvious after the CRA correction by adjusting displacement and rotation (Fig. 2(b)), the ETS for longer (48–72 h) lead-time rainfall forecasts has increased significantly. The displacement error of the CRA method, reflecting the horizontal shift between the observed and forecasted rain centers, is directly related to the error of forecasted typhoon track, which increases with lead time. Therefore, the ETS for 48–72 h lead-time rainfall forecasts could be largely improved after the CRA shift, which further indicates that the typhoon track error is an important reason for the decrease of ETS for longer-time rainfall forecasts. The ascending trend of ETS for heavy rainfall amounts (>100 mm) after CRA shift adjustment is more obvious than that of lower rainfall amounts. The variation trend of POD with different forecast lead times and rainfall thresholds is generally consistent with that of ETS (Figs. 2(c) and 2(d)). Similarly, after the displacement adjustment of the CRA method, the FAR also decreases significantly, especially for large rainfall thresholds (Figs. 2(e) and 2(f)).

4.2 Rainfall forecasts in different lifetime of Lekima during landfall

To more intuitively express the life cycle of landfall typhoons, relative time (totally 4 days) will be used for the time axis (Fig. 4). Taking the landfall time of Lekima (at 1800 UTC on August 9, 2019) as the reference time (0 h, as shown in Fig. 4), 24 h before the landfall is −24 h, and 24/48/72 h after the landfall is 24/48/72 h. The numbers of TC forecasts from −24 h before landfall to 72 h after landfall are generally comparable except for the 250 mm rainfall (figure not shown).

Fig.4 (a), (c), and (e) The rainfall forecast verification measures of ETS, POD, and FAR during different stages of the typhoon life cycle (from 24 h before the landfall to being weakened to a tropical depression). (b), (d), and (f) The same as (a), (c), and (e), but for the results after the displacement adjustment based on the CRA method. “−24 h,” “0 h,” “24 h,” “48 h”, and “72 h” shown on the x-axis represent the time of 24 h before the landfall ( −24 h), at the time of landfall (0 h), and 24/48/72 h after the landfall (24/48/72 h), respectively. The solid, dashed, dotted lines and the squares represent 30 mm, 50 mm, 100 mm, and 250 mm rainfall, respectively.

Full size|PPT slide

As can be seen in Fig. 4(a), for 30 mm rainfall forecasts, the variation of ETS before and after the landfall is not significant, which basically remains at about 0.3. Moreover, the ETS even rises slightly 48 h after the landfall, and then decreases significantly. For 50 mm rainfall forecasts, the variation trend of ETS is similar to that of 30 mm rainfall, but the value decreases slightly. However, for 100 mm rainfall, the ETS declines very quickly to 0.1, 0.05 and 0 at 24, 48, 72 h after the landfall, respectively. Moreover, the ETS is even lower for the 250 mm rainfall. As shown in Fig. 1, Lekima continued to move northward to the inland after landfall, and produced heavy rainfall in Shandong region at 24–72 h after landfall. Therefore, for the current model, the forecast ability for persistent heavy rainfall is limited during the stage when the typhoon moved further inland.

From contribution ratios of various errors for all rainfall thresholds (Fig. 5) the R always contributes less than 10% of the total error, further decreasing after landfall and becoming negligible. At landfall (0 h), the major error comes from both P and D, but the D significantly decreases and the P increases when moving inland (24–48 h). For the 30–50 mm rainfall, the P contributes more than 50% of the total error. However, for the rainfall larger than 100 mm, the proportion of V also increases significantly, which may be related to the typhoon vortex structure and the interaction with westerlies during this period (Fig. 6). More detailed analysis will be provided in section 5. At the end of the TC lifecycle (72 h after the landfall), the proportion of V continues to increase for 30–50 mm rainfall forecasts, indicating that the model still has a significant deviation in the rainfall volume forecasts after the typhoon encountered the westerlies. For the 250 mm rainfall forecasts, the proportion of D exceeds 60%, which might be attributed to the fact that small-scale rainfall center errors will have a significant impact on rainfall errors.

Fig.5 The same as Fig. 3, but for the rainfall during different landfall stages. “−24 h,” “0 h,” “24 h,” “48 h”, and “72 h” shown on the x-axis represent the time of 24 h before the landfall ( −24 h), at the time of landfall (0 h), and 24/48/72 h after the landfall (24/48/72 h), respectively.

Full size|PPT slide

Fig.6 The 24 h accumulated rainfall observations and forecasts of typhoon Lekima during (a) 2019080900–2019081000 UTC and (b) 2019081000–2019081100 UTC, respectively. The merged rainfall from the observation (the 1st column), the forecasted rainfall from the model (the 2nd column), the forecasted rainfall after the displacement adjustment based on the CRA method (the 3rd column), and the merged and shifted model rainfall based on the CRA method (the 4th column) are shown.

Full size|PPT slide

Figures 4(b), 4(d), and 4(f) show the performance of the modified rainfall forecast results by the displacement adjustment of the CRA method. For all rainfall thresholds, the ETSs of the TC rainfall forecasts during different periods of its lifecycle have increased by various degrees, particularly for large thresholds. However, except for the 250 mm rainfall, the ETSs of other rainfall thresholds almost keep the original variation trend, which shows that the decrease of ETS after landfall could not simply be related to the TC track errors, but might correlate with other factors. In addition, it shows similar results for the POD and FAR.

To illustrate the forecast errors in TC heavy rainfall clearly, the 24 h 50 mm rainfall forecasts are taken as an example (considering that the 24 h rainfall forecasts have the best performance compared with other lead times), and the differences of rainfall forecasts in stages of near-landfall and inland are shown in Table 1 and Fig. 6. Near landfall for 50 mm rainfall forecast (Fig. 6(a)), the number of observed grid points with precipitation is more than that of the forecast, whereas the average daily precipitation rate is close (Table 1). The observed maximum daily precipitation rate is 124 mm·d⁻¹, while the forecast is 159 mm·d⁻¹, indicating that the range of forecasted rainfall is relatively small, while the total rainfall amount is relatively consistent, but with a higher rainfall intensity. From the statistical results of error decomposition, the errors of this rainfall forecast are mainly D and P, accounting for nearly 50% of the total error. The correlation coefficient between the model forecast field and the observation field is 0.41, with a root mean square (RMS) error of 61 mm·d⁻¹. After the displacement adjustment based on the CRA method, the correlation coefficient between the two rain fields rises to 0.7, and the RMS error also drops to 43.6 mm·d⁻¹, which indicates that after the displacement correction for the model rainfall, the rainfall forecast errors caused by the maximum rain center errors can be effectively reduced.

Tab.1 CRA verification for 24 h accumulated rainfall forecasts (Fcst) with 24 h lead time

	2019080900–2019081000				2019081000–2019081100
Statistics on CRA grid	Obs		Fcst		Obs		Fcst
Girdpoints>50 mm·d−1	3357		2912		2798		2416
Average rain rate/(mm·d−1)	73		73		77		51
Maximum rainrate/(mm·d−1)	457		275		312		158
90% rain rate/(mm·d−1)	124		159		162		82
Rain volume/km3	38		38		35		24
	Fcst	Move		M+R	Fcst	Move		M+R
Correlation coefficient	0.41	0.7		0.7	0.06	0.25		0.28
RMSE/(mm·d−1)	61	43.5		43.5	69.8	64.3		63.6
	2019080900–2019081000				2019081000–2019081100
Skill scores	Fcst		Obs		Fcst		Obs
Probability of detection (POD)	0.59		0.64		0.53		0.52
False alarm ratio (FAR)	0.34		0.26		0.4		0.4
Equitable threat score (ETS)	0.36		0.44		0.32		0.31
Error decomposition	2019080900–2019081000				2019081000–2019081100
Rotation error	0%				2%
Displacement error	49%				11%
Volume error	0%				13%
Pattern error	51%				74%

For inland rainfall, the forecasted rainfall could reflect the two observed rain centers in Shandong and northern Zhejiang (Fig. 6(b)), with the proportion of D being only 11% (Table 1). However, the model forecasted rainfall volume is obviously weaker, especially for the heavy rainfall in Shandong. The average observed precipitation rate is 77 mm·d⁻¹, while the forecast value is 51 mm·d⁻¹, and the rainfall V accounts for 13%. The observed maximum precipitation rate and 90 percentile precipitation rate are 312 mm·d⁻¹ and 162 mm·d⁻¹, respectively, which are both nearly twice the forecasted values. This indicates a weaker rainfall intensity and narrower distribution range in the forecasts. As seen from the scatter diagram, the corresponding forecasted rainfall errors increase significantly compared to the observations. P accounts for 74% of the total error, indicating significant errors in forecasted rainfall distribution.

5 Synoptic diagnostic analysis of forecast errors

To understand some physical reasons for the rainfall forecast errors, synoptic diagnoses are carried out for the near-land and inland rainstorm processes (located in Zhejiang and Shandong, respectively; Fig. 1(a)).

5.1 Precipitation process near landfall

Figure 7 shows the observed and GRAPES-TCM forecasted precipitation processes just before and after landfall. Before landfall (2019080900–2019080912 UTC), the forecasted rainfall intensity is significantly higher than the observation (Fig. 7(b)), and the forecasted heavy rainfall follows the TC center, showing a clear symmetric ring structure. However, the observed rainfall is distributed more asymmetrically, with the rainfall maximum located on the east side of the typhoon center. After landfall (2019080918–2019081000 UTC), the observed heavy rain is mainly located in the coastal area, the central and eastern Zhejiang. Due to the influence of coastline and terrain, the observed heavy rain becomes more asymmetric, but the forecasted heavy rain maintains a relatively symmetric distribution, resulting in excessive concentration of precipitation in the TC center region. The asymmetric distribution could not be predicted, leading to weak forecasted or even missing rainfall in the northern Zhejiang. Therefore, the observed rainfall intensity and asymmetric distribution near landfall cannot be predicted by the current model.

Fig.7 (a) Observed and (b) forecasted 6-h rainfall (unit: mm·h−1; initial forecast time: 1200 UTC on August 8, 2019) during the landfall (2019080818–2019081000 UTC). The origin of the figure is the typhoon center. The x- and y-axis shows the grid points away from the typhoon center. The grid distance is 0.1°.

Full size|PPT slide

Previous studies have shown that the asymmetric distribution of typhoon rainfall is related to many factors, such as environmental vertical wind shear (VWS), asymmetric convergence caused by the boundary layer friction, and the distribution of asymmetric environmental water vapor or convection (Shapiro, 1983; Chen and Yau, 2001; Yu et al., 2015, 2017). In general, the environmental VWS, TC movement and TC intensity are closely connected with the rainfall pattern of landfalling TC, so it is worth analyzing the influence of these three factors on the distribution of the forecast rainfall.

Through diagnostic analysis on environmental VWS from both GRAPES-TCM forecast and ECMWF analysis fields (Fig. 8(a)), it is found that near landfall (2019080900–2019081000UTC), the forecasted VWS shows large differences from the analysis in both magnitude and direction, compared to the inland period (2019081012–2019081112UTC). At the initial time of 2019080812 UTC of GRAPES-TCM, the forecasted VWS is north-westerly with the magnitude greater than 5 m/s, and it is similar to the analysis. The observed rainfall in 2019080812–2019080818 UTC distributes asymmetrically with the maximum rainfall in the down shear-left of the analysis VWS. The distribution of the forecasted rainfall is analogous to the observation but with higher rainfall intensity and denser structure. Yu et al. (2017) investigated the relationship between rainfall intensity and distribution in LTCs over China, and found that the rainfall axisymmetry is closely related to the TC intensity. In general, stronger TCs have higher averaged rain rate and higher averaged amplitude of axisymmetric rainfall. Comparing the forecasted MSLP (mean sea level pressure) of typhoon center and its vortex intensity with the observations before landfall (Figs. 1(b), 8(b), and 8(c)), we found that the forecasted TC intensity is much stronger at all times, resulting in more axisymmetric rainfall. This indicates that the initialization scheme adopted in the typhoon model could not adequately represent the typhoon intensity, adversely affecting the asymmetric rainfall distribution pattern. At 2019080900 UTC, the VWS magnitude of the analysis is larger than the forecasted value. An asymmetric rainfall structure is also evident, as shown in Fig. 7(a). Soon after TC landfall, the magnitude of the analysis VWS becomes sufficiently weak (less than 5 m/s) so as not to influence the rainfall distribution, while the forecast VWS retains its strength. Thereafter when the TC moves further inland and interacts with the westerly trough (Figure not shown), the magnitude of VWS gradually increases, as does the forecasted value (Fig. 8(a)).

Fig.8 Comparison between the ECMWF initial analysis and the GRAPES-TCM forecast (initial forecast time: 1200 UTC on August 8, 2019) during the landfall. (a) Temporal change of environmental VWS between 200 and 850 hPa (unit: m·s−1; red: ECMWF analysis; blue: the GRAPES-TCM forecast; arrow: the VWS direction; black vertical dashed line: the landfall time of typhoon). (b) and (c) are the 850 hPa relative vorticity (unit: 10−5·s−1) from the ECMWF analysis and GRAPES-TCM forecast, respectively. The origin of the figure is the typhoon center. The x- and y-axis show the grid points away from the typhoon center. The grid distance in (b) and (c) is 0.25° and 0.1°, respectively.

Full size|PPT slide

5.2 Precipitation process during moving inland and northward

After the northward movement of the typhoon, Shandong and the Bohai Bay became the inland heavy rainfall center (Fig. 9(a)). After 2019081012 UTC, the typhoon center left the northern Zhejiang and entered the southern Jiangsu. However, the large-scale heavy precipitation occurred in the central and southern Shandong in the typhoon periphery, with the average rainfall rate generally exceeding 9 mm·h⁻¹. Meanwhile, the maximum value at the center reached nearly 30 mm·h⁻¹, which persisted for nearly 18 h and then gradually declined. After 2019081118 UTC, the main rainfall center moved eastward to the east of the Bohai Bay. The model has predicted heavy rainfall in the typhoon periphery (Fig. 9(b)), but the forecasted rain intensity is obviously weak, and the rain center is slightly displaced from the observation. The overall intensity of model rainfall forecast during 2019081100–2019081106 UTC is weak. After 2019081106 UTC, the heavy rainfall area (>9 mm·h⁻¹) became more scattered, mainly distributed in the coastal areas of the Bohai Bay. The model generally reflects the rainfall distribution in these periods, but deviations in the rainfall pattern and volume still exist, which can also be seen from the contribution rates of different errors in the late stage of the typhoon life cycle (Fig. 5).

Fig.9 The same as Fig. 7, but for (a) observed (top two panels) and (b) forecasted rainfall (bottom two panels) at an interval of 6 h during the TC weakening stage when the typhoon moved further inland during 2019081018–2019081200 UTC (unit: mm·h−1 , initial forecast time: 100 UTC on August 10, 2019). The dashed and solid lines (shown in Fig. 9(a) and Fig. 9(b), respectively) represent the observed and forecasted TC tracks for the every 6-h rainfall period.

Full size|PPT slide

According to the ECMWF analysis data, the TC peripheral circulation interacts with the westerlies during the northward movement of TC (Fig. 10(a)). The baroclinicity during this period becomes an important physical parameter of rainfall, which is closely related to frontogenesis and equivalent potential temperature advection. In this paper the frontogenesis function is defined following Srock and Bosart (2009) as$F_{\mathrm{d}\mathrm{i}\mathrm{v}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}})\left|{\nabla }_h{\theta }_e\right|$, where $u$ and $v$ are respectively zonal and meridional wind components and ${\theta }_e$ is equivalent potential temperature. The equivalent potential temperature advection is defined as, which is a scalar quantity with negative (positive) value representing cold and dry (warm and moist) advection.

At 2019081012 UTC, a large frontogenesis area with the magnitude of the frontogenesis function greater than 0.1 × 10⁻⁸ K·m⁻¹·s⁻¹ distributed in Shandong (Fig. 10(a)). This should be related to the active interaction between two airflows with different${\theta }_e$. At 2019081100 UTC, a strong and narrow front in the south-north direction passed across the middle part of Shandong (Fig. 10(a)). It corresponds well with the rainfall center greater than 10 mm·h⁻¹during 2019081018–2019081100 UTC. At the same time, the intensities of cold-dry and warm-moist air advections reached their peaks (Fig. 10(c)). Though the GRAPES-TCM forecasted advection has a similar pattern to that of the ECMWF analysis (Fig. 10(d)), it is much weaker, especially in the northern heavy rain area. In other words, the forecasted front is also weaker than the real one. Another problem with the forecasted front involves its influence range. It is found that there exist major discrepancies in the rainfall patterns. This maybe directly caused by the forecasted front location.

Fig.10 (a) and (b) The 850 hPa frontogenesis (unit: 10−8 K·m−1·s−1) from the ECMWF analysis and the GRAPES-TCM forecast (initial forecast time: 1000 UTC on August 10, 2019), respectively. (c) and (d) The 850hPa equivalent potential temperature advection (unit: 10−1 K·s−1) from the ECMWF analysis and the GRAPES-TCM forecast, respectively.

Full size|PPT slide

6 Discussion and conclusions

The deep understanding of the numerical rainfall forecast ability and the main error source is essential to the operational precipitation forecast and the improvement of TC models. Based on the object-oriented verification method of CRA, the rainfall forecast performance of a regional numerical model, the GRAPES-TCM developed by the STI/CMA, was analyzed for super typhoon Lekima in 2019 during its landfall life cycle (from 24 h before landfall to being weakened to a tropical depression).

After landfall, Lekima maintained for several days over land, and produced very heavy rainfall during its northward movement in China. The heavy rainfall process had two stages, near-landfall and moving inland. The sources and differences of the model rainfall forecast errors in the two stages were comparatively analyzed, and a comparison with the ECMWF initial analysis was also performed. The major sources of rainfall forecast errors in different landfall stages were confirmed. On this basis, existing problems in forecasting the thermodynamic evolution of the model were preliminarily revealed. The main conclusions are as follows.

Analysis of the variable rainfall forecast ability during different landfall stages has shown that the TC heavy rain forecast ability for the inland (24–72 h after landfall) is much weaker than that just near landfall, and the reasons for the numerical model forecast errors are also different. Though the forecast ability for the threshold of heavy rainfall (>100 mm) has been obviously improved with the displacement adjustment based on the CRA method, the inland rainfall forecast ability is still largely lower than that in the near landfall. It indicates that the decreasing trend of rainfall forecast skill scores after landfall has little correlation with the typhoon track errors, but might be related to other factors, such as the interactions with external environmental fields.

For the inland rainfall, the northward-moving typhoon interacted with the westerlies. The forecast rainfall intensity is much smaller. The results show that the forecasted frontogenesis process is much weaker than the observation. Also, the different influence ranges of the cold-dry and warm-moist air flows could cause a discrepancy in the thermodynamics of the air advections. The rainfall forecast errors could be related to the bias in physical processes and in the interactions between the typhoon and westerlies.

However, near landfall the observed heavy rainfall shows an obvious asymmetric distribution, with the maxima located east of the TC center. By contrast, in the forecasted rainfall field, the heavy rainfall shows a notably symmetric distribution. From the VWS and MSLP analysis before the landfall of Lekima, it is found that the forecasted TC intensity is much stronger than the observation, accompanied by lower VWS effects, thus resulting in a more axisymmetric rainfall structure. The model has not effectively simulated the asymmetric rainfall pattern near landfall, which might be related to the selection of the initial assimilation scheme or the description of relevant physical processes in the typhoon model.

For TC Lekima rainfall forecasts with different lead times, the 24 h rainfall forecast has the highest scores, and the decline of the 48–72 h rainfall skill score is directly related to the typhoon track errors, which can be seen from the rainfall forecast improvement by the CRA shift adjustments. The variation trend of the 30–50 mm rainfall forecast scores is relatively consistent for different rainfall thresholds. When the rainfall threshold exceeds 100 mm, the verification score drops significantly.

The rainfall error analysis based on the CRA method shows that for all forecast lead times and rainfall thresholds, the forecast errors are mainly from pattern (P) and displacement (D), and the proportion of D also gradually increases with the increase of rainfall threshold. However, after moving inland, the proportion of D significantly decreases while that of volume (V) increases, especially for rainfall above 100 mm.

It is worth noting that, currently in the field of model precipitation verification, traditional methods of classified verification based on grid points or stations such as ETS, POD, FAR, or continuity verification methods such as deviation, RMS error, and correlation, are still widely used. Although these verification methods can evaluate the forecast ability of the model from different perspectives, they cannot determine the sources of errors, not to mention quantitative analysis on the proportions of various errors. Therefore, the CRA method was used in this paper to verify the rainfall forecast performance of regional model of GRAPES-TCM when typhoon Lekima struck at eastern coastal China, and to confirm the main sources of forecast errors in different landfall stages. The results in this study will be helpful for forecasters and model developers to carry out researches on model precipitation forecasting, or to adjust the assimilation schemes or physical processes in the model in a targeted way.

References

[1] Alpert P (1986). Mesoscale indexing of the distribution of orographic precipitation over high mountains. J Appl Meteorol Climatol, 25(4): 532–545

[2] Baker W L, Veblen T T (1990). Spruce beetles and fires in the nineteenth-century subalpine forests of Western Colorado, U.S.A. Arct Alp Res, 22(1): 65–80

[3] Bearup L A, Maxwell R M, Clow D W, McCray J E (2014). Hydrological effects of forest transpiration loss in bark beetle-impacted watersheds. Nat Clim Chang, 4(6): 481–486

[4] Bentz B J, Regniere J, Fettig C J, Hansen E M, Hayes J L, Hicke J A, Kelsey R G, Negron J F, Seybold S J (2010). Climate change and bark beetles of the western United States and Canada: direct and indirect effects. Bioscience, 60(8): 602–613

[5] Bergen J D (1971). Vertical profiles of wind speed in a pine stand. For Sci, 17(3): 314–321

[6] Bethlahmy N (1974). More streamflow after a bark beetle epidemic. J Hydrol (Amst), 23(3 – 4): 185–189

[7] Bethlahmy N (1975). A Colorado episode: beetle epidemic, ghost forests, more streamflow. Northwest Sci, 49(2): 95–105

[8] Biederman J A, Brooks P D, Harpold A A, Gochis D J, Gutmann E, Reed D E, Pendall E, Ewers B E (2014a). Multiscale observations of snow accumulation and peak snowpack following widespread, insect-induced lodgepole pine mortality. Ecohydrology, 7(1): 150–162

[9] Biederman J A, Harpold A A, Gochis D J, Ewers B E, Reed D E, Papuga S A, Brooks P D (2014b). Increased evaporation following widespread tree mortality limits streamflow response. Water Resour Res, 50(7): 5395–5409

[10] Biederman J A, Somor A J, Harpold A A, Gutmann E D, Breshears D D, Troch P A, Gochis D J, Scott R L, Meddens A J H, Brooks P D (2015). Recent tree die-off has little effect on streamflow in contrast to expected increases from historical studies. Water Resour Res, 51(12): 9775–9789

[11] Boon S (2007). Snow accumulation and ablation in a beetle-killed pine stand in Northern Interior British Columbia. BC J Ecosyst Manag, 8(3): 1–13

[12] Boon S (2009). Snow ablation energy balance in a dead forest stand. Hydrological Processes, 23(18): 2600–2610

[13] Boon S (2012). Snow accumulation following forest disturbance. Ecohydrology, 5(3): 279–285

[14] British Columbia Ministry of Forests (1995). Bark beetle management guidebook. Forest practices code. Victoria, BC. Forest Practices Branch.

[15] Carpenter S R, Mooney H A, Agard J, Capistrano D, DeFries R S, Diaz S, Dietz T, Duraiappah A K, Oteng-Yeboah A, Pereira H M, Perrings C, Reid W V, Sarukhan J, Scholes R J, Whyte A (2009). Science for managing ecosystem services: beyond the Millennium Ecosystem Assessment. Proc Natl Acad Sci USA, 106(5): 1305–1312

[16] Chen F, Zhang G, Barlage M, Zhang Y, Hicke J A, Meddens A, Zhou G, Massman W J, Frank J (2015). An observational and modeling study of impacts of bark beetle-caused tree mortality on surface energy and hydrological cycles. J Hydrometeorol, 16(2): 744–761

[17] Cheng J D (1989). Streamflow changes after clear-cut logging of a pine beetle-infested watershed in southern British Columbia, Canada. Water Resour Res, 25(3): 449–456

[18] Clow D W, Rhoades C, Briggs J, Caldwell M, Lewis W M Jr (2011). Responses of soil and water chemistry to mountain pine beetle induced tree mortality in Grand County, Colorado, USA. Appl Geochem, 26: S174–S178

[19] Colorado Water Conservation Board (CWCB) (2010). Colorado’s Water Supply Future: State of Colorado 2050 Municipal and Industrial Water Use Projections.

[20] Eisenhart K, Veblen T T (2000). Dendrochronological detection of spruce bark beetle outbreaks in northwestern Colorado. Can J Res, 30(11): 1788–1798

[21] Fassnacht S R, Dressler K A, Bales R C (2003). Snow water equivalent interpolation for the Colorado River Basin from snow telemetry (SNOTEL) data. Water Resour Res, 39(8), doi: 10.1029/2002WR001512

[22] Feng X, Taylor S, Renchaw C E, Kirchner J W (2002). Isotopic evolution of snowmelt 1. A physically based one-dimensional model. Water Resour Res, 38(10): 1–8

[23] Fettig C J, Klepzig K D, Billings R F, Munson S A, Nebeker T E, Negron J F, Nowak J T (2007). The effectiveness of vegetation management practices for prevention and control of bark beetle infestations in coniferous forests of the western and southern United States. For Ecol Manage, 238(1 – 3): 24–53

[24] Gat J R (2010). Isotope Hydrology: A Study of the Water Cycle. Series on Environmental Science and Management, 1–69

[25] Hazlett S, Hansen N (2007). Current status of the mountain pine beetle in the west—A literature review. USDA Forest Service, P.A.09.2009. Washington, D.C. 1–4

[26] Hermann A, Stichler W (1981). Isotope input into runoff systems from melting snow covers. Nord Hydrol, 12: 309–318

[27] Ingraham N L (1998). Isotopic Variations in Precipitaiton. In: Kendall C, McDonnell J J, eds. Isotope Tracers in Catchment Hydrology. Amsterdam: Elsevier Science B.V., 87–118

[28] International Atomic Energy Agency (IAEA) (1997). Technical procedure for cumulative monthly sampling of precipitation for isotopic analysis.

[29] Leatherman D A (2008). The Health of Colorado’s Forests: High Elevation Forests. Edited by Colorado Department of Natural Resources. Colorado State University, Fort Collins, CO. p. 17

[30] Livneh B, Deems J S, Buma B, Barsuguli J J, Schneider D, Molotch N P, Wolter K, Wessman C A (2015). Catchment response to bark beetle outbreak and dust-on-snow in the Colorado Rocky Mountains. J Hydrol (Amst), 523: 196–210

[31] Love L D (1955). The effects on stream flow of the killing of spruce and pine by the Engelmann Spruce beetle. Trans Am Geophys Union, 36(1): 113–118

[32] Maggart A (2014). Effects of mountain pine beetle caused tree mortality on streamflow generation mechanism in Colorado. MS Thesis Catchment Science, Colorado State University, Fort Collins, CO

[33] Man G (2012). Major Forest Insect and Disease Conditions in the United States: 2011. USDA Forest Service, Report No: FS-1000

[34] Maxwell J (2009). Forest headwaters-saving our water source. USDA Forest Service Regional Office Denver, CO. Unpublished memo. 6pp

[35] Menger A L (2015). Response of streamflow and stream chemistry to pine beetle induced tree mortality across northern Colorado. MS Thesis Watershed Science, Colorado State University. Fort Collins, CO

[36] Mikkelson K M, Bearup L A, Maxwell R M, Stednick J D, McCray J E, Sharp J O (2013a). Bark beetle infestation impacts on nutrient cycling, water quality, and interdependent hydrological effects. Biogeochemistry, 115(1 – 3): 1–21

[37] Mikkelson K M, Maxwell R M, Ferguson I, Stednick J D, McCray J E, Sharp J O (2013b). Mountain pine beetle infestation impacts: modeling water and energy budgets at the hill-slope scale. Ecohydrology, 6(1): 64–72

[38] Morris J L, Cottrell S, Fettig C J, Hansen W D, Sherriff R L, Carter V A, Clear J L, Clement J, DeRose R J, Hicke J A, Higuera P E, Mattor K M, Seddon A W R, Seppӓ H T, Stednick J D, Seybold S J (2016). Managing bark beetle impacts on social-ecological systems: priority questions to inspire future research. J Appl Ecol, doi: 10.1111/1365-2664.12782

[39] Parnell A, Jackson A (2015). SIAR: Stable Isotope Analysis in R (Version 4.2). [Software].

[40] Potts D F (1984). Hydrologic impacts of a large-scale mountain pine beetle (Dendroctonus ponderosae) epidemic. Water Resour Bull, 20(3): 373–377

[41] Pugh E, Gordon E (2013). A conceptual model of water yield impacts from beetle-induced tree death in snow-dominated lodgepole pine forests. Hydrol Processes, 27(14): 2048–2060

[42] Pugh E, Small E (2012). The impact of pine beetle infestation on snow accumulation and melt in the headwaters of the Colorado River. Ecohydrology, 5(4): 467–477

[43] Raffa K F, Aukema B H, Bentz B J, Carroll A L, Hicke J A, Turner M G, Romme W H (2008). Cross-scale drivers of natural disturbances prone to anthropogenic amplification: the dynamics of bark beetle eruptions. Bioscience, 58(6): 501–517

[44] Rhoades C C, McCutchan J H, Cooper L A, Clow D, Detmer T M, Briggs J S, Stednick J D, Veblen T T, Ertz R M, Likens G E, Lewis W M (2013). Biogeochemistry of beetle-killed forests: explaining a weak nitrate response. Proc Natl Acad Sci USA, 110(5): 1756–1760

[45] Seidl R, Schelhaas M J, Lexer M J (2011). Unraveling the drivers of intensifying forest disturbance regimes in Europe. Glob Change Biol, 17(9): 2842–2852

[46] Spittlehouse D (2007). Influence of the mountain pine beetle on the site water balance of lodgepole pine forests. In: Proceedings Mountain Pine Beetle and Catchment Hydrology Workshop: Preliminary Results of Research from BC, Alberta and Colorado. 25–26

[47] Stednick J D (1996). Monitoring the effects of timber harvest on annual water yields. J Hydrol (Amst), 176(1 – 4): 79–95

[48] Stednick J D, Troendle C A (2016). Chapter 12. Hydrologic effects of forest management practices on water resources. In: Forest Hydrology: Processes, Management, and Assessment. CABI Publishing, 192–203

[49] Taylor S, Feng X, Renshaw C, Kirchner E (2001). Isotopic evolution of snowmelt: validation and parameterization of a one-dimensional model using controlled snowmelt experiments. Water Resour Res, 37: 759–769

[50] Taylor S, Feng X, Williams M, McNamara J (2002). How isotopic fractionation ofsnowmelt affects hydrograph separation. Hydrol Processes, 16(18): 3683–3690

[51] Teti P (2007). Solar Radiation and Snow Ablation in Natural and Managed Pine Stands. In: Proceedings Mountain Pine Beetle and Catchment Hydrology Workshop: Preliminary Results of Research from BC, Alberta and Colorado. 17–18

[52] Thomson A J, Shrimpton D M (1984). Weather associated with the start of mountain pine-beetle outbreaks. Can J Res, 14(2): 255–258

[53] USDA (2014). Aerial Detection Survey [Data file].

[54] USGS Land Cover Institute (2005). Conifer Forest Layer

[55] USGS National Elevation Dataset (2013). The National Map

[56] Veblen T T, Hadley K S, Nel E M, Kitzberger T, Reid M, Villalba R (1994). Disturbance regime and disturbance interactions in a Rocky Mountain subalpine forest. J Ecol, 82(1): 125–135

[57] Weed A S, Ayres M P, Hicke J A (2013). Consequences of climate change for biotic disturbances in North American forests. Ecol Monogr, 83(4): 441–470

[58] Wehner C E (2016). Effect of mountain pine beetle kill on streamflow generation mechanisms. MS Thesis Watershed Science, Colorado State University. Fort Collins, CO

[59] Winkler R D (2007). Snow accumulation and melt in southern interior lodgepole pine forests. In: Proceedings Mountain Pine Beetle and Catchment Hydrology Workshop: Preliminary Results of Research from BC, Alberta and Colorado. 19–20

Acknowledgements

Support for this project was made available by the National Science Foundation grant 1204460: Water Quality and Supply Impacts from Climate-Induced Insect Tree Mortality. The work expands upon work done by Ariann Maggart and Ashley Menger. The authors also acknowledge the useful comments provided by the anonymous reviewers.

RIGHTS & PERMISSIONS

2017 Higher Education Press and Springer-Verlag Berlin Heidelberg