본 연구는 전국 19개 농업기후지대 중 7개 이상의 시군구를 갖고, 지리적으로 이격되어있는 지대를 선정하였다. 이를 통하여 총 세 개의 상이한 농업기후지대(중북부 내륙, 소백서부 내륙, 영남 내륙, Choi and Yun, 1989)에 걸친 27개 시군구가 연구대상 지역으로 설정되었으며, 18년(2000년 ~ 2017년)에 걸친 시군구별 논벼 수량성 및 기상자료를 수집하였다(Fig. 1).

Fig. 1. 
A map of Korean mainland sigungus, with the three surveyed agroclimactic zones (as defined by Choi and Yun. 1985) colored in green, brown and orange

이앙이후 논벼 전체 생육주기(5월 말 ~ 9월 말, Table 2 참조)에 걸쳐서 중북부내륙지대(이하 ‘북부지역’)는 평균 27.68℃와 17.35℃의 일중 최고･최저온도를 보인 반면, 소백서부내륙지대(이하 ‘중부지역’)는 28.43℃와 18.35℃의 평균 최고·최저 온도를, 영남 내륙지대(이하 ‘남부지역’)는 평균 28.95℃와 19.34℃의 최고·최저 온도를 보였다. 생육주기 중 연평균 누적강수량은 각각 1061 mm, 901 mm, 893 mm로 상이하였다. 본 연구에서 종속변수로 활용한 수량성의 18년 평균값은 북부지역 627 kg/10 a(최대 696 kg/10 a, 최소 594 kg/10 a), 중부지역 평균 688 kg/10 a(최대 710 kg/10 a, 최소 662 kg/10 a), 남부지역 평균 653 kg/10 a(최대 666 kg/10 a, 최소 636 kg/10 a)로 각 지역 내에서 비교적 균일하였다.

본 연구는 통계청의 『농작물생산조사』를 참조하여 18년(2000년 ~ 2017년)에 걸친 분석대상 27개 시군구의 논벼 수량성 데이터를 수집하였다(단위: kg/10 a). 해당 데이터는 시군별 표본농가에 대한 현장실측을 통하여 기록되며, 표본추출틀은 매년 갱신된다(Statistics Korea 2019).

본 연구는 기상청(KMA 2020)의 1 km 단위 격자 기상데이터(gridded weather data)를 활용하였다. 해당 데이터는 실측자료에 해양도, 고도, 경사방향 등의 미기상적인(micrometeorological) 요소를 반영 후 재생산된 정밀기상데이터로, 이는 농림축산식품부 제공 팜맵(AGIS 2020)과 결합되어 시군구별 논벼 생산지 기상자료를 형성하는 데에 활용되었다. 즉, GIS를 활용하여 각 격자에 포함된 논벼 재배면적을 전체 시군구 논벼 재배면적 대비 비율로 분할한 뒤, 이를 가중치로 활용하여 최저온도, 최고온도 및 강수량의 가중평균을 도출하였다. 이와 같이 시군구 논벼 재배지의 지리적 분포를 감안한 유효기온 및 강수량의 도출은 선행연구(Schlenker and Roberts 2009)의 핵심적인 특성 중 하나로, 단일 관측소에서의 관측치를 활용하는 경우와 달리 기상변수의 공간적 분포를 반영한다는 장점을 갖는다.

예컨대 어떤 시군구의 논벼 재배지역이 5개 격자에 걸쳐 분포하고, 팜맵 조회 결과 각 격자 내 논벼 면적이 0.2 km², 0.3 km², 0.5 km², 0.4 km². 그리고 0.6 km²였다면 해당 시군구의 총 논벼 재배면적은 2 km²였으므로, 각각의 격자 기상데이터에 0.1, 0.15, 0.25, 0.2, 그리고 0.3의 가중치를 곱하여 시군구 내 논벼재배지의 평균 기상 데이터 값을 구하는 것이다. 이와 같이 추출된 기온변수는 기상청의 실측자료와 비교하여 평균 0.996(최저 0.985, 최고 0.998)의 Pearson 상관계수를 보였다.

본 연구에서 활용한 기상변수는 작물생육기간 중 기상변수와 극한기상 변수로 분류된다. 전자에는 선행연구(Schlenker and Roberts 2009)를 참조하여 온도구간별 누적소요시간과 강수량이 포함되었으며, 후자에는 기상청이 기후변화 시나리오 하에서 제공 중인 극한기후지수 중 강수량과 관련된 3종(5일 최다강수량, 최대 무강수지속기간, 호우일수)이 포함되었다.

극한기상의 수량성에 대한 영향력을 평가하기 위하여 각 격자별로 집중호우(일 강수량 80 mm 이상) 발생 여부를 확인한 뒤, 이 역시 앞선 방식과 동일하게 가중평균하여 해당 시군구 논 재배면적 중 집중호우 피해가 발생한 면적비중을 계산하였다.²⁾ 이외의 5일간 최대 누적강우량의 최대치와 연속 무강수일 계산은 시군구 단위로 종합된 강수량 데이터를 활용하였다(Table 1).

본 연구는 논벼의 생육주기별로 변화하는 기온 민감성을 반응하기 위하여 논벼가 노지 조건에 노출되는 시기³⁾를 분얼활착기(이하 ‘1기’), 유수형성 이후 출수개화 전(이하 ‘2기’), 그리고 등숙기(이하 ‘3기’)로 구분하였다. 또한, 지역별(북부, 중부, 남부)로 별도의 분석을 시도하여 기온조건에 대한 지역별 차이 역시 확인할 수 있도록 하였다.

논벼의 생육주기는 품종의 조만성(早晩性), 이앙시기 및 작기의 기후특성에 따라 달라진다. 본 연구는 분석의 일관성을 위하여 국립종자원의 정부 보급종 중 각 지역에 적합한 것으로 판단된 품종(북부지역 ‘오대’, 중부지역 ‘삼광’, 남부지역 ‘새일미’, KSVS, 2020)을 선택, 각 품종의 품종보고서(NICS, 2004; 2006; 2012)를 참고하여 생육주기를 선정하였다(Table 2).

전술된 바와 같이, 본 연구는 생육주기 내 온도구간별 소요일수를 독립변수로 활용하여 논벼 수량성을 추정한다. 온도구간을 설정하기 위해서는 구간설정의 시작온도와 종결온도를 정하는 것이 필요한데, 선행연구에서는 결빙점 이하 온도를 배제하는 방식을 선택하였다(Schlenker and Roberts, 2009)⁴⁾. 본 연구에서는 각 지역의 온도분포에 맞춰진 온도구간을 선정하기 위해 세 지역별로 각 생육주기의 일 최저온도와 일 최고온도를 집계하였다. 이어서 시작온도와 종결온도 간의 구간이 전체 집계일의 99.9%를 포함할 수 있도록 구간을 선정하였다(Table 3).

즉, 북부지역에서는 생육주기별로 각각 10개, 8개, 10개의 3도 너비(width) 온도구간을 설정하고, 중부지역에서는 9개, 8개, 11개의 온도구간을, 남부지역에서는 10개, 8개, 11개의 온도구간을 선정하되, 각 구간은 지역적 온도 분포를 반영하도록 했다.

이외의 기상변수(강수량 및 극한기상 3종)는 생육주기별로 서로 다른 계수를 가질 수 있도록 설정하였다. 끝으로, 시간에 따른 논벼 재배기술의 변화를 포함하기 위하여 추세항과 추세제곱항을 추가하였다.

계량경제학적 모형으로는 선행연구를 참조하여 미관측 효과 패널모형(unobserved effects panel model)을 활용하였다. 해당 모형의 추정은 대상 지역의 시군구별로 기상데이터와 추세항으로는 확인되지 않는 시간불변 특성(예컨대 토성(土城))이 존재하여 논벼 수량성에 영향을 미친다는 가정 하에 실시된다. 이 미관측 특성이 독립변수(기상, 극한기상 및 추세항)와 상관성을 갖는지 여부에 따라 해당 모형은 고정효과 모형 또는 임의효과 모형의 형태로 구체화된다(Wooldridge, 2010). 시군구별 미관측 특성의 존부를 확인하기 위해 F-test를 우선 실시하였으며, 해당 특성과 독립변수간의 상관성을 검정하기 위해 Hausman 검정을 실시하였다. 이를 종합한 모형의 형태는 아래와 같다.

y_it는 i 시군구 t년도에서의 수량성을, c_i는 i 시군구의 불변특성으로 인한 효과⁵⁾를, h_s,j,i,t는 i 시군구 t년도 s 생육주기 j번째 온도구간에서의 누적소요일수를 의미한다. p_s,i,t, HR_s,i,t, cumprcp_s,i,t, dconsec_s,i,t는 각각 s기 i 시군구 t년도에서의 강수량, 집중호우 피해면적 비율, 최대 5일 누적강수량, 최대 무강수일을 의미한다. T_t과 T_t²는 각각 추세항과 추세제곱항을 뜻하고, ε_i,t는 오차항을 지칭하며, E(ϵ_i,t|x_i, c_i)가정한다. 해당 수식에서 각 생육주기별 온도구간의 수(10개, 8개, 10개)는 북부지역의 사례이며, 중부지역과 남부지역의 생육주기별 온도구간 수는 상이하다.

이와 같이 계산된 모형의 오차항은 시간적(serial) 및 횡단면적(cross-sectional) 상관성을 가질 수 있다. 다시 말해, 특정 시군구에서 t년도의 수량성 충격은 t′년도(t≠t′)의 수량성 충격과 영향을 가질 수 있으며(시간적 상관성), 특정 연도에서 i 시군구의 수량성 충격은 i′ 시군구(i≠i′)의 수량성 충격과 상관성을 가질 수 있다(횡단면적 상관성). 각 상관성의 검정에는 Breusch-Godfrey 검정(Breusch 1978; Godfrey, 1978)과 Pesaran의 CD검정(Pesaran, 2004)을 활용하였으며, 귀무가설 기각 기준은 p-value 0.05였다. 오차항의 상관성이 존재할 경우 엄밀 표준편차(robust standard deviation)의 추정이 필요한데, 이를 위해서는 White와 Driscoll and Kraay 기법을 활용하였다(White, 1984; Driscoll and Kraay, 1998).

F-test 결과 세 지역 모두에서 0에 가까운 p-value로 시군구별 미관측 특성의 부재가 기각되었다. 미관측 특성변수의 형태를 결정하기 위해 실시한 Hausman 검정에서는 세 지역 모두에서 1에 가까운 확률로 귀무가설이 유지되어, 시군구별 미관측 특성이 독립변수와 무관함을 상정하는 임의효과 모형이 선정되었다.

오차항의 자기상관성 검정에서 북부지역에서는 0.26의 p-value로 시간적 자기상관성의 부재 가설이 유지되었으며, 횡단면적 자기상관성 역시 0.27의 p-value로 기각되지 아니하였다. 중부지역에서는 각각 0.005와 0.07의 p-value로 시간적 자기상관성 부재만 기각되었다. 남부지역에서는 각각 0.039와 0.033의 p-value로 두 자기상관성 부재가 모두 기각되었다.

추정 결과는 Table 4와 같다. 각 온도구간별 소요시간의 형태로 표현된 기온변수는 세 지역 모두에서 강수와 극한기상에 비하여 큰 중요성을 갖는 변수로 선정되었다.

세 지역은 각 생육주기 기온조건의 역할에 대해 판이한 해석을 제시하였다. 북부지역에서는 기온변수가 뚜렷한 영향력을 갖는 시기가 부재하였던 반면, 중부지역에서는 3기가, 남부지역에서는 2기가 수량성에 긍정적으로 기여하는 시기인 것으로 드러났다.

이와 같은 수량성 반응에서의 차이를 감안하면, 저온피해･고온피해 취약성 여부는 각 지역 수준에서 판단하는 것이 적합하다. 다시 말해, 세 지역이 각 온도구간에 대해 갖는 계수의 범위가 상이하므로, 온도장해 발생은 장해온도 구간에서의 계수 추정치가 동일 생육주기 내 다른 온도구간의 그것과 비교하였을 때 차이를 보이는지의 여부를 통하여 확인되어야 한다.

예컨대, 남부지역의 경우(Fig. 2)에서는 각 생육주기의 온도구간이 다른 온도구간에 비하여 어떠한 계수를 갖는가를 통하여 저온피해와 고온피해의 존부가 판단되어야 하는데, Fig. 2를 통하여 남부지역에서는 1기, 2기, 3기 모두 최말온도구간에서 고온피해가 관측되고, 저온피해는 3기에서만 관측되었음을 확인할 수 있다. 같은 맥락에서, 북부지역에서는 1기 고온피해, 3기 저온피해가 뚜렷하며, 중부지역에서는 1기 저온피해, 2기 저온피해가 드러난다.

Fig. 2. 
An illustration of the yield-temperature relationship, in the case of Southern region

각각 생육주기를 살펴보면, 북부지역과 남부지역은 1기 최말온도구간(bin 10, 각각 32℃ ~ 35℃, 35℃ ~ 38℃)에서 고온피해의 규모가 증가하는 것으로 나타났는데, 이는 33℃ 이상의 온도에서 논벼의 분얼발생이 억제된다(Yoshida 1973)는 선행연구와 일치한다. 또한, 냉해피해가 중부지역에서 발견된 구간(8℃ ~ 11℃) 역시 분얼활착기 논벼의 냉해피해 임계온도로 알려진 13℃(RDA 2014) 이하로, 농학분야 연구와 상당부분 일치한다.

2기에서는 냉해피해온도(15℃ ~ 17℃, RDA 2014) 이하 온도구간에서 중부지역에서 냉해피해가 드러났다. 이는 중부지역 이는 삼광벼 재배지역에서의 유수형성기 ~ 출수기 간 냉해피해가 임실률에 미치는 영향에 대한 농학분야 선행연구(Kim et al., 2009)와 일치한다.

3기에서의 고온피해온도는 남부에서 드러났는데, 이는 등숙기 하의 고온조건이 수량성에 부정적인 영향을 가진다는 선행연구(Lee et al., 2013)와 일치한다. 또한 냉해피해가 북부와 남부에서 관측된 온도구간(6℃ ~ 9℃, 5℃ ~ 8℃) 역시 선행연구의 등숙기 냉해피해 임계온도인 14℃(RDA, 2014) 이하로, 농학분야 연구와 일관적이다.

강수는 모든 지역과 생육주기에서 음의 계수를 가졌으나, 그 크기가 매우 작아(전체 논벼 생육주기 연평균 강수량 평균 960 mm) 수량성에 결정적인 영향을 미치지 못하는 것으로 보인다.

극한기상에 있어서, 집중호우 발생면적 변수(HR)는 유의성을 보인 경우 생육주기에서 부적합한 계수(즉 양의 계수)를 보였으나, 집중호우 강도(cumprcp) 변수는 유의성을 보인 지역 및 주기에서 음의 계수를 보였다. 한발해는 유의성을 보이지 아니하였다. 전반적으로 극한기상항의 제거는 북·중·남부 지역에서 각각 0.005, 0.037, 0.044의 Adj. R² 하락을 유발하였다.

추세항 및 추세제곱항은 중·남부지역에서 유의성을 보였는데, 전자는 기준연도(2000년)로부터 8.5년 이후 극소점을 달성하고 이후 상승하는 추세를 선정한 반면, 후자는 기준년으로부터 13년 이후 극대점을 달성하고 이후 하락하는 추세를 선정하였다.