2.1 Winter wheat yields

Our primary data consists of 1,024 county-level administrative regions from 145 cities in 21 provinces in China, covering the period from 2001 to 2016. We measure the county-level average winter wheat yields using the satellite imagery from Iizumi and Sakai (Reference Iizumi and Sakai2020) based on the Chinese county-level administrative regions. The satellite imagery of crop yields has been used in multiple analyses of yields (see, e.g., Zhang et al., Reference Zhang, Niu and Yu2021; Sun et al., Reference Sun, Zhao, Zhang, Chen, Ying, Lv, Che and Wang2022). Figure 1 presents the geographical locations of winter wheat growing regions in 2001, 2006, 2011 and 2016. As can been seen in figure 1, winter wheat is mostly cultivated in Northern China, such as in Shanxi Province, Shandong Province, Hebei Province and Henan Province, showing a good geographical coverage over a large territory in China.

Figure 1.

Winter wheat yields in China, (a) Winter wheat yield in 2001, (b) Winter wheat yield in 2006, (c) Winter wheat yield in 2011, (d) Winter wheat yield in 2016. Notes: Solid blue lines represent boundaries of province-level administrative regions. Dashed black lines represent boundaries of county-level administrative regions. Red pixels represent high volume of winter wheat yields while green pixels represent low volume of winter wheat yields.

Table 1 presents variable definitions and summary statistics. As shown in panel A, the mean and standard deviation of winter wheat yieldsFootnote ³ are 4,386.26 kg/hectare and 3,065.71 kg/hectare, respectively. This indicates a high yields variability in our sample, which is consistent with large territory coverage in figure 1. In addition to winter wheat yields, we match the county-level winter wheat yields with county-level climatic characteristics, fertilizers and soils. The rest of this section includes our data-construction methods for climatic variables as well as fertilizer and soil texture.

Table 1.

Summary statistics

2.2 Climatic characteristics

We use daily temperature to construct the FDDFootnote ⁴ to investigate the effect of freezing stress occurrence on winter wheat yields (Zhu et al., Reference Zhu, Kim, Jin, Lin, Wang, Ciais, Mueller, Aghakouchak, Huang, Mulla and Makowksi2022). Additionally, we also use daily temperature to construct the GDDs in spring (March to May), autumn (September to November) and winter (December to February) following equation (1) from Zhu et al. (Reference Zhu, Kim, Jin, Lin, Wang, Ciais, Mueller, Aghakouchak, Huang, Mulla and Makowksi2022) to control for the thermal effect during the winter wheat growing seasons. As is shown in panel B of table 1, the mean and standard deviation of FDD are 225.97$^\circ$C and 372.88$^\circ$C, respectively. For GDDs in autumn, winter and spring, we find that the mean and standard deviation for high GDD terms across three seasons are lower compared to low and medium GDD terms.Footnote ⁵

To investigate the effect of snowpack insulation expected to reduce the freezing stress on winter wheat yields, we construct the SCF using FDD as well as daily snow depth (SND) as a ratio of days when both FDD>0 and SND$>2$ cm, to days when FDD$>0$, following equation (2) from Zhu et al. (Reference Zhu, Kim, Jin, Lin, Wang, Ciais, Mueller, Aghakouchak, Huang, Mulla and Makowksi2022).Footnote ⁶ In addition to the insulation effect, snow is one of the important water storage sources that provides soil moisture. In order to control for the effect of water supply from snowmelt, we partitioned precipitation into rainfall and snowfall to represent water supply from different sources on the basis of daily temperature and daily precipitationFootnote ⁷ following equation (3) from Zhu et al. (Reference Zhu, Kim, Jin, Lin, Wang, Ciais, Mueller, Aghakouchak, Huang, Mulla and Makowksi2022). Panel B in table 1 shows that the mean and standard deviation of SCF are 0.18 and 0.25, respectively. This indicates that on average 18 per cent of days with freezing stress occurrence (FDD$>0$) have SND greater than 2 cm. We also find that average snowfall in winter is 19.84  mm, which is higher than counterparts in spring (4.48 mm) and autumn (4.63 mm). Differently from snowfall, average rainfall in spring (178.40 mm) and autumn (159.8 mm) is higher than that in winter (49.97 mm). The difference between rainfall and snowfall in winter is consistent with intuition that snow is heavier in winter, leading to the fact that a higher proportion of precipitation is from snowfall instead of rainfall in winter than in autumn and spring.

2.3 Fertilizer and soil

As has been shown in previous literature, winter wheat yields are significantly affected by nitrogen fertilizer rates (Litke et al., Reference Litke, Gaile and Rua2018; Tabak et al., Reference Tabak, Lepiarczyk, Filipek-Mazur and Lisowska2020) and soil textures (Thaler et al., Reference Thaler, Eitzinger, Trnka and Dubrovsky2012; Johnen et al., Reference Johnen, Boettcher and Kage2014). We include the county-level nitrogen fertilizer ratesFootnote ⁸ and proportions of soil textures,Footnote ⁹ such as clay particles, sand particles and soil organic carbon components, to control for the effect of nitrogen fertilizer use and soil components on winter wheat yields (Zhu et al., Reference Zhu, Kim, Jin, Lin, Wang, Ciais, Mueller, Aghakouchak, Huang, Mulla and Makowksi2022). As shown in panel C of table 1, the mean and standard deviation of the nitrogen fertilizer rate are 193.11 kg/hectare and 80.21 kg/hectare, respectively. The mean and standard deviation of clay particles, sand particles, and organic carbon components are 1,101.91 g/kg and 311.73 g/kg, 1,326.94 g/kg and 549.28 g/kg, and 508.62  dg/kg and 305.02 dg/kg, respectively.Footnote ¹⁰

3.1 Conceptual framework

Our primary empirical strategy relies on quantile regression for panel data (QRPD), which is a special case of the generalized quantile estimator in which the estimator addresses a foundational problem of including fixed effects (Powell, Reference Powell2014a, Reference Powell2020, Reference Powell2022). Our choice of an empirical framework stems from the nature of the research objective and variables of interest. Our primary objective is to explore the heterogeneous impact of snowpack insulation on winter wheat yields. This indicates a potentially varying effect at different points in the conditional distribution of the outcome variable. Furthermore, the snowpack exhibits significant inter- and intra-annual variability (Mankin and Diffenbaugh, Reference Mankin and Diffenbaugh2015). Therefore, the variability is a crucial aspect of evaluating the impact of snowpack insulation on crop yields (Huning and AghaKouchak, Reference Huning and AghaKouchak2020; Qin et al., Reference Qin, Abatzoglou, Siebert, Huning, AghaKouchak, Mankin, Hong, Tong, Davis and Mueller2020).

Previous literature has shown that heterogeneous effects are useful information, but usually missed by mean regression techniques (Bitler et al., Reference Bitler, Gelbach and Hoynes2006; Powell, Reference Powell2014a; Powell and Wagner, Reference Powell and Wagner2014). Moreover, in recent advancements, there has been an impetus to combine the advantages of panel data regression as a tool to account for the unobserved heterogeneity with the ability of quantile regression to estimate robust estimates accounting for outliers (Koenker, Reference Koenker2004; Canay, Reference Canay2011; Galvao Jr, Reference Galvao2011; Ponomareva, Reference Ponomareva2011; Kato et al., Reference Kato, Galvao and Montes-Rojas2012). However, combining the quantile and panel regression, these methods exploit additive fixed effects in the quantile function at the cost of loss of generality in the quantile regression inference and interpretation. In a novel approach, Powell (Reference Powell2022) extends this approach to build non-additive fixed effects with conventional quantile regression advantages such as a nonseparable disturbance term.Powell (Reference Powell2022) provides an estimator in an instrumental variable framework for fixed effects that are consistent with the small time period and are feasible to estimate in standard econometric software.Footnote ¹¹ Powell (Reference Powell2022) presented an estimator which utilizes fixed effects for identification purposes, while the results are interpreted as traditional cross-sectional quantile regression estimates. Fixed effects enable arbitrary correlation with the level of variables of interest, accounting for unobserved county-level differences in non-parametric specification.

In summary, our proposed empirical strategy has two key advantages. Firstly, our dataset shows that the winter wheat yields exhibit positive skewness, which indicates a potentially nonlinear relationship.Footnote ¹² Previous researchers have argued that climatic variables such as temperature and rainfall have nonlinear relationships with crop yields (Schlenker and Roberts, Reference Schlenker and Roberts2009; Burke et al., Reference Burke, Hsiang and Miguel2015; Proctor, Reference Proctor2021). Therefore, we believe that the investigation of the winter wheat yields should take explicitly into account the distribution of outcome variables. However, several applied economic research studies regarding the effect of snowpack insulation on winter wheat yields do not explicitly model nonlinear relationships. We suggest that the quantile relationship is better suited for modeling climatic variables such as snowpack insulation and winter wheat yields. Furthermore, our proposed approach is robust to accounting for outliers in the data (Cameron and Trivedi, Reference Cameron and Trivedi2022; Powell, Reference Powell2022).

Secondly, we estimate the heterogeneous relationship between winter wheat yields and snowpack insulation with an explicit quantile function. The approach is semiparametric in terms of not having parametric distribution of regression errors and its departure from the normality assumptions (Cameron and Trivedi, Reference Cameron and Trivedi2022). Moreover, the nonlinear relationship between winter wheat yields and snowpack insulation has varying effects throughout the outcome distribution, and previous authors have found that snowpack insulation heterogeneously impacts crop yields (Li et al., Reference Li, Qi, Smith, Pattey and Qian2022; Zhu et al., Reference Zhu, Kim, Jin, Lin, Wang, Ciais, Mueller, Aghakouchak, Huang, Mulla and Makowksi2022). To our knowledge, we are one of the first to estimate the effect of snowpack insulation using a distribution of crop yields.

3.2 Model specification

We begin with an estimation of the two-way fixed effects (TWFE) model to estimate the relationship between the winter wheat yields and snowpack insulation effect as well as other climatic characteristics alongside the agro-economic variables such as fertilizer, and soil attributes. Our main estimating equation is specified as follows:

(1)\begin{equation} \begin{aligned} {\mathrm{Wheat Yields}}_{i,t} & = \alpha _1{\mathrm{FDD}}_{i,t} + \alpha _2{\mathrm{FDD}}_{i,t} \times {\mathrm{SCF}}_{i,t} + \mathop {\sum }\limits_{j = 1}^3 \alpha _{3,j}{\mathrm{Rainfall}}_{i,t}^j \\ & \quad+\mathop {\sum }\limits_{j = 1}^3 \alpha _{4,j}{\mathrm{Snowfall}}_{i,t}^j+ \mathop {\sum }\limits_{j = 1}^3 \alpha _{5,j}{\mathrm{GDD}}_{i,t}^j + \alpha _{11}{\mathrm{Fert}}_i + \alpha_{12}{\mathrm{Sand}}_i \\ & \quad+ \alpha _{13}{\mathrm{Clay}}_i + \alpha _{14}{\mathrm{OM}}_i + {\mathrm{procitydis}}_i + {\mathrm{time}}_{t} + e_{i,t}, \end{aligned} \end{equation}

where $\mathrm {Wheat Yields}$ is the winter wheat yields. The ${\mathrm {FDD}}_{i,t}$ is the FDD in county $i$ at year $t$. The ${\mathrm {FDD}}_{i,t} \times {\mathrm {SCF}}_{i,t}$ is the interaction term between FDD and SCF in county $i$ at year $t$. The ${\mathrm {Snowfall}}_{i,t}^j$, ${\mathrm {Rainfall}}_{i,t}^j$ and ${\mathrm {GDD}}_{i,t}^j$ are the snowfall, rainfall and GDD in county $i$ at year $t$ for season $j$, respectively. The ${\mathrm {Fert}}_i$, ${\mathrm {Sand}}_i$, ${\mathrm {Clay}}_i$ and ${\mathrm {OM}}_i$ are the nitrogen fertilizer rates, clay particles, sand particles and soil organic carbon components in county $i$, respectively. The ${\mathrm {procitydis}}_i$ is county FEs, controlling for any time-invariant unobserved heterogeneity at the county level. The ${\mathrm {time}}_{t}$ denotes the year FEs.Footnote ¹³ We primarily focus on the estimates of $\alpha _1$ and $\alpha _2$ which represent the relationship between winter wheat yield with FDD and FDD*SCF.

In TWFE, the estimates reflect the average relationship, which is further refined in quantile regression as introduced by Koenker and Bassett Jr (Reference Koenker and Bassett1978). The following equation illustrates how the impact of the variable of interest changes based on total distribution:

\[ y_i=\mathbf{d}_i^{\prime} \beta\left(u_i^*\right). \]

To build on the framework from Powell and Wagner (Reference Powell and Wagner2014), the structural quantile function (SQF) can be written as

\[ S_{wheat yield_i}=\mathbf{d}_i^{\prime} \beta(\tau), \quad \tau \in(0,1). \]

The SQF provides the $\tau$ th quantile of $y_{\mathbf {d}}=\mathbf {d}_i^{\prime } \beta (u^*)$ for a fixed $\mathbf {d}$ and a randomly-selected $u^* \sim U(0,\,1)$. In our context, the SQFs provide the distributions of the outcome variable winter wheat yield. More specifically, we aim to account for covariates and fixed effects as shown in equation (1) as follows:

(2)\begin{equation} \begin{aligned} S_{wheat yield_{it}} & =\alpha _1{\mathrm{FDD}}_{i,t} + \alpha _2{\mathrm{FDD}}_{i,t} \times {\mathrm{SCF}}_{i,t} + \mathop {\sum }\limits_{j = 1}^3 \alpha _{3,j}{\mathrm{Rainfall}}_{i,t}^j \\ & \quad+ \mathop {\sum }\limits_{j = 1}^3 \alpha _{4,j}{\mathrm{Snowfall}}_{i,t}^j+ \mathop {\sum }\limits_{j = 1}^3 \alpha _{5,j}{\mathrm{GDD}}_{i,t}^j + \alpha _{11}{\mathrm{Fert}}_i + \alpha _{12}{\mathrm{Sand}}_i \\ & \quad+ \alpha _{13}{\mathrm{Clay}}_i + \alpha _{14}{\mathrm{OM}}_i + {\mathrm{procitydis}}_i + {\mathrm{time}}_{t} + e_{i,t} , \quad \tau \in(0,1). \end{aligned} \end{equation}

We present equation (2) building on Powell (Reference Powell2014a) with the aim of estimating the $\alpha _1$ and $\alpha _2$, which represents the relationship between the winter wheat yield with FDD and FDD*SCF.

5.1 Alternative snow cover fraction

In the previous section, we constructed SCF as a ratio of days when both FDD$>0$ and SND$>2$ cm. To investigate whether threshold selection is a concern, we report estimates from alternative thresholds of SND, such as 5 cm, 10 cm, 15 cm and 20 cm. Figure A6 (online appendix) shows estimates for FDD and FDD*SCF with thresholds of SND for 5  cm, 10 cm, 15 cm and 20 cm. As can be seen in figure A6, the estimates of FDD in panels (a), (b), (c) and (d) are significantly negative, which is consistent with the results from figure 2. As SND thresholds are higher, the estimates for FDD*SCF are predominantly significantly positive. Though the estimates for FDD*SCF in figure A6 are different from the results in figure 2, the snowpack insulation effect still persists to offset the negative effect of FDD. The difference is due to a higher threshold of SND, which makes the effect of snowpack insulation stronger.

5.2 Alternative growing degree days

In the previous section, we constructed season-specific GDDs following thresholds used in Zhu et al. (Reference Zhu, Kim, Jin, Lin, Wang, Ciais, Mueller, Aghakouchak, Huang, Mulla and Makowksi2022). To alleviate the concern about threshold selection, we construct new season-specific GDDs following thresholds used in Tack et al. (Reference Tack, Barkley and Nalley2015). The thresholds to categorize low, medium and high GDDs are 0$^\circ$ C/18$^\circ$ C/23$^\circ$C for autumn, 0$^\circ$ C/5$^\circ$ C/10$^\circ$C for winter, and 0$^\circ$ C/17$^\circ$ C/28$^\circ$C for spring. Furthermore, we also construct new non-season-specific GDDs following the thresholds of 0$^\circ$C/26$^\circ$C/35$^\circ$C used in Zhu et al. (Reference Zhu, Kim, Jin, Lin, Wang, Ciais, Mueller, Aghakouchak, Huang, Mulla and Makowksi2022). In figure A7 in the online appendix, panel (a) shows estimates for FDD and FDD*SCF using new season-specific GDDs and panel (b) shows estimates for FDD and FDD*SCF using new non-season-specific GDDs. Results from figure A7 show that FDD is estimated to have a statistically significant negative effect, which is consistent with findings from figure 2. From panel (a) in figure A7, FDD*SCF is estimated to have a statistically significant positive effect at higher quantiles of the winter wheat yields distribution (from the 75th to 95th quantiles), which is aligned with the results from figure 2. From panel (b) in figure A7, FDD*SCF is estimated to have a predominantly statistically significant positive effect, showing the existence of the snowpack insulation effect.

6.1 Economic development

Previous literature has shown that nightlight is a good proxy for economic development (see, e.g., Chen and Nordhaus, Reference Chen and Nordhaus2011; Nordhaus and Chen, Reference Nordhaus and Chen2015; Gibson et al., Reference Gibson, Olivia, Boe-Gibson and Li2021). We construct the annual average of nightlight intensity at county-level administrative regions from 2001 to 2016Footnote ¹⁹ and categorize regions as low and high nightlight regions based on the median of nightlight intensity. Figure 5 shows QRPD estimation results for FDD and FDD*SCF for low and high nightlight regions. Panel (a) shows that FDD is estimated to have a statistically significant negative effect on winter wheat yields for regions with low nightlight intensity but a positive effect for regions with high nightlight intensity throughout the entire distribution of winter wheat yields. Panel (b) shows that FDD*SCF is estimated to have a statistically significant negative effect at lower and middle quantiles of the winter wheat yields distribution (from the 25th to 65th quantiles), and a positive effect at higher quantiles of the distribution (from the 70th to 95th quantiles). Comparing the results in figure 5 to the results in figure 2, we find that the snowpack insulation effect persists at higher quantiles of winter wheat yields distribution for regions with low nightlight intensity (FDD$<0$ and FDD*SCF$>0$), indicating that a heterogeneous effect of snowpack insulation exists. One of the possible explanations is that agricultural activities are less likely to be invested in regions with low economic development (i.e., low nightlight intensity), preserving the original natural effect of snowpack insulation.

Figure 5.

Heterogeneity analysis: nightlight, (a) QRPD estimation results: FDD, (b) QRPD estimation results: FDD*SCF. Notes: The dependent variable in all regressions is the mean of winter wheat yields. The top two panels show QRPD estimation results of FDD for low and high nightlight intensity. The bottom two panels show QRPD estimation results of FDD*SCF for low and high nightlight intensity.

6.2 Labor input

As is shown in Tongwei et al. (Reference Tongwei, Luo, Choy, Li and He2020) and Meiyan (Reference Meiyan2010), labor input is one of the important factors in agricultural productivity. To investigate the heterogeneous effect of snowpack insulation based on labor input, we construct the annual average population density at the county-level administrative regions from 2001 to 2016Footnote ²⁰ and categorize them as low and high population density regions based on the median of population density. Figure 6 shows QRPD estimation results for FDD and FDD*SCF for low and high population density regions. Panel (a) shows that FDD is estimated to have a statistically significant negative effect on winter wheat yields for low population density regions but a statistically significant positive effect for high population density regions throughout the entire distribution of winter wheat yields. Panel (b) shows that FDD*SCF is estimated to have a significantly negative effect on winter wheat yields throughout the entire distribution with exceptions at the 90th and 95th quantiles for low population density regions. Comparing FDD and FDD*SCF estimates in figure 6 to those in figure 2, we find that the snowpack insulation effect persists at higher quantiles of the winter wheat yields distribution in low population density regions, indicating a heterogeneous effect of snowpack insulation exists. One of the possible explanations is that farmers are less likely to increase labor input in winter wheat cultivation in low population density regions, keeping the original natural effect of snowpack insulation.

Figure 6.

Heterogeneity analysis: population density. Notes: The dependent variable in all regressions is the mean of winter wheat yields. The top two panels show QRPD estimation results of FDD for low and high population density. The bottom two panels show QRPD estimation results of FDD*SCF for low and high population density.

6.3 Irrigation intensity

As shown in previous scientific literature (Sun et al., Reference Sun, Liu, Zhang, Shen and Zhang2006; Liu et al., Reference Liu, Wiberg, Zehnder and Yang2007; Zeng et al., Reference Zeng, Yao, Zhang, Yang, Bai, Zhang, Wang and Wang2021), irrigation plays an important role in winter wheat yields. We construct the annual average irrigation area at the county-level administrative regions from 2001 to 2016Footnote ²¹ and categorize regions as low and high irrigation regions based on the median of irrigation area. Figure 7 presents QRPD estimation results for FDD and FDD*SCF for low and high irrigation regions. Panel (a) shows that FDD is estimated to have a statistically significant negative effect on winter wheat yields for low irrigation regions but a statistically significant positive effect for high irrigation regions throughout the entire distribution of winter wheat yields. Panel (b) shows that FDD*SCF is estimated to have a significantly positive effect on winter wheat yields at higher quantiles of winter wheat yields distribution for both low (90th and 95th quantiles) and high (80th to 95th quantiles) irrigation regions. Based on FDD and FDD*SCF estimates in low irrigation regions which are consistent with estimates of FDD and FDD*SCF from figure 2, we find that the snowpack insulation effect persists at higher quantiles of the distribution for the low irrigation regions, suggesting a heterogeneous effect of snowpack insulation exists. One of the possible explanations is that high irrigation regions are more likely to be affected by artificial irrigation, distorting the original natural effect of snowpack insulation.

Figure 7.

Heterogeneity analysis: irrigation intensity. (a) QRPD estimation results: FDD, (b) QRPD estimation results: FDD*SCF. Notes: The dependent variable in all regressions is the mean of winter wheat yields. The top two panels show QRPD estimation results of FDD for low and high irrigation intensity. The bottom two panels show QRPD estimation results of FDD*SCF for low and high irrigation intensity.

6.4 Climate classifications

As shown in figure 1, winter wheat is cultivated in a large territory in China, covering different climatic regions from warm to snow as well as humid to arid. To investigate the heterogeneous effect of climate classifications, we split regions into three categories: arid, temperature and perhumid.Footnote ²² Figure 8 shows QRPD estimation results for FDD and FDD*SCF for the climate of arid, temperature and perhumid. As can be seen in panel (a), FDD is estimated to have a statistically significantly negative effect on winter wheat yields throughout the entire distribution for all three climate classifications. Panel (b) shows that FDD*SCF is estimated to have a predominantly statistically significant positive effect on winter wheat yields for the climate of arid and temperature but a predominantly negative effect for perhumid. Based on FDD and FDD*SCF estimates in figure 8, we find that the snowpack insulation effect persists for the climate of arid and temperature, indicating a heterogeneous effect of snowpack insulation for different climate classifications.

Figure 8.

Heterogeneity analysis: climate classifications. (a) QRPD estimation results: FDD, (b) QRPD estimation results: FDD*SCF. Notes: The dependent variable in all regressions is the mean of winter wheat yields. The top three panels show QRPD estimation results of FDD for arid, warm temp and snow climate classifications. The bottom two panels show QRPD estimation results of FDD*SCF for arid, warm temp and snow climate classifications.