2.1 Study Area and Data Collection

본 연구는 한국 전역의 하천을 대상으로 하며, 줄날도래속(Hydropsyche spp.) 3종(Hydropsyche kozhantschikovi, Hydropsyche valvata, Hydropsyche formosana)의 유충에 대한 출현유무와 환경변수(고도, 기온)를 바탕으로 서식처 적합도(Habitat Suitability)를 예측하고자 하였다. 한강수계관리위원회의 “기후변화가 수생태계에 미치는 영향과 대응전략(‘10–’12)”, 환경부의 “환경생태유량 산정기준 연구 및 시범산정(‘14)”, 환경부⋅국립환경과학원의 “수생태계 건강성 조사 및 평가(‘12–’13)”, “하천 수생태계 현황 조사 및 건강성 평가(‘14–’21)”, “수생태계 참조하천 선정 및 활용방안 마련 연구(‘15–’18)”, “4대강 보 개방에 따른 수생태계 변화 조사(‘18–’20)”, 국립환경과학원의 “생물측정망 모니터링 및 평가기법 개발연구(‘17–’18)”의 지점 또는 정점 조사 결과 중 2010–2021년 동안 수집된 총 3,333개 고유 지점의 결과를 기반으로 하였다(Fig. 1). 고유 지점은 한강, 낙동강, 금강, 섬진강, 영산강, 제주 수계를 포함한 전국에 분포하며, 0-888 m의 고도 내에서 다양한 서식환경이 포함되었다.

Fig. 1. Geographic distribution of 3,333 study sites surveyed across South Korea from 2010 to 2021.

2.2 Environmental Variables

2.3 Future Scenario Prediction

미래 분포 예측에는 기후변화 시나리오 중 Representative concentration pathway (RCP) 8.5를 적용하였다. RCP 시나리오는 IPCC 제5차 평가보고서(AR5)에서 사용된 기후변화 시나리오 체계로, 2100년까지의 방사 강제력(Radiative forcing) 수준에 따라 RCP 2.6, 4.5, 6.0, 8.5의 네 가지 시나리오로 구분된다. 이 중 RCP 8.5는 온실가스 감축이 이루어지지 않고 높은 배출이 지속되는 경로로, 기후변화가 극단적으로 진행될 가능성을 고려한 최악(Worst case) 시나리오에 해당한다^{(IPCC, 2013)}. 본 연구에서 온실가스 배출이 지속되어 기후변화가 극단적으로 진행되는 시나리오인 RCP 8.5를 적용한 이유는, 이러한 시나리오에서 생태적 취약성 평가 및 보전 전략 수립에 있어 잠재적 위험의 상한선을 평가할 수 있기 때문이다^{(Schwalm et al., 2020)}. 선행 연구에 따르면, 기후변화 시나리오에 따라 향후 수온은 전반적으로 상승할 것으로 예측된다. ^{Kim et al. (2022)}은 한강 유역을 대상으로 RCP 2.6–8.5 시나리오를 적용한 결과, 기온 상승폭이 최대 +3.51℃에 달하며, 이에 따라 수온은 기온 1℃당 평균 0.796℃씩 증가할 것으로 분석하였다.

^{Ahn and Han (2010)}은 RCP 8.5 시나리오를 바탕으로 전국 하천과 호소의 수온 변화를 예측한 결과 2100년 기준으로 하천 수온은 +2.45℃, 호소 수온은 +2.14℃ 상승할 것이라고 분석하였다. 한편, ^{Yu et al. (2021)}은 팔당호의 15년간 관측 데이터를 기반으로 RCP 8.5 시나리오를 적용하여 수온이 연평균 0.07℃씩 증가하고, 2070년에는 누적 상승 폭이 약 +1.8℃에 이를 것으로 전망하였다. 또한 ^{Lee et al. (2016)}은 금강 상류 용담댐 유역에 대해 HSPF 유역모형을 적용한 결과, 2080년대에 금강 기준 +0.5에서 +1.2℃, 진안 지역에서는 최대 +1.6℃까지 상승할 것으로 분석하였으며, 수온 증가는 주로 겨울과 봄철에 두드러질 것으로 예측하였다. 본 연구에서는 이러한 국내 예측값을 참고하여, 2100년을 기준으로 하천 수온이 +2.45 ℃ 상승하는 조건을 반영하였다. 산정된 수온은 HSI 분석, 랜덤 포레스트 모델의 입력변수로 활용되었으며, 동일한 구조를 갖는 예측 알고리즘을 통해 2100년 시점의 줄날도래속 3종의 분포를 예측하였다.

2.4 Habitat Suitability Modeling

줄날도래속 3종의 서식처 적합도 분석을 위해, 각종의 출현빈도 데이터를 기반으로 평형수온을 일정 간격으로 구간화하고, 구간별 상대 출현 빈도를 집계하여 PMF를 도출하였다. 이후 PMF에 7개 확률분포모형(Generalized exponential, Weibull, Inverse weibull, Lognormal, Generalized logistic power, Gamma, Beta)을 적용하여 확률밀도함수(Probability density function, PDF)를 근사하였다. 각 모형은 수평 이동 역치 값을 포함한 것이며, 표준화 평균제곱근오차(Normalized root mean square error, NRMSE)가 가장 낮은 값을 보인 모형을 최적으로 평가하였다^{(Kong et al., 2025)}. 선정된 분포모형의 확률밀도함수는 각종의 수온 반응 특성과 적정 서식 범위를 시각화하는 데 활용되었다. HSI 값은 PDF 값을 최빈값(Mode)에 해당하는 PDF 값으로 나누어 산정하였다. 즉 평형수온을 T, 최빈 수온을 T ^ , 각 종의 상대출현빈도에 대한 PDF를 f(T)라 할 때, HSI 함수인 f_HSI는 다음과 같다.

아울러, 현재 및 미래(+2.45°C 시나리오) 조건에서의 수온 구간별 가용 면적은 외부 공간자료 기반의 수온-면적 분포를 추정함으로써 산정되었고, 각종의 HSI 곡선과 중첩하여 공간적 서식 가능성과 겹침을 시각화하였다. 미래 조건에 대한 면적 함수 예측은 다음과 같은 방법론에 따라 수행되었다. 먼저, 현재 조건에서의 수온 구간별 면적 분포는 0.5°C 간격으로 구간화하였고, 평균 기온 상승량 +2.45°C를 각 수온 구간에 등량평행이동(Horizontal translation) 방식으로 적용함으로써, 단순화된 형태의 미래 구간별 면적 함수를 구축하였다. 이때 면적의 확률밀도 함수를 f_A라 하면, 현재 및 미래 조건의 종별 중첩 HSI 값(Combined HSI)은 다음과 같다.

2.5 Random Forest Modeling

줄날도래속 3종에 대한 서식처 적합도 예측을 위해 랜덤 포레스트(Random Forest, RF)를 적용하였다. 이 머신 러닝 기법은 분류(Classification) 및 회귀(Regression) 문제에 널리 활용되는 앙상블 학습 기법으로, ^{Breiman (2001)}에 의해 제안되었다. 이 알고리즘은 훈련 데이터로부터 다수의 의사결정나무(Decision Tree)를 생성하고, 이를 기반으로 분류의 경우 다수결 투표(Majority Voting), 회귀의 경우 평균 예측값(Averaged prediction)을 통해 최종 결과를 도출한다. 랜덤 포레스트는 각 나무마다 서로 다른 부트스트랩 샘플(Bootstrap sample)을 이용하고, 노드 분할 시 무작위로 선택된 변수 집합 중에서 최적의 분할 기준을 선택하여 모델 간 상관성을 낮추고 과적합(Overfitting)을 방지한다. 이러한 구조는 변수 간 상호작용이 복잡한 환경에서도 우수한 예측 성능을 나타내며, 국내 생태학 분야에서도 효과적인 서식처 적합성 예측 기법으로 활용된 바 있다^{(Kim et al., 2021)}.

본 연구에서는 종별로 개별 모델을 구성하였으며, 입력변수(Predictors)로는 고도와 평형수온 관계식^{(Kang, 2024)}을 사용하였다. 고도는 평형수온 산정에 일부 반영되지만, 유속, 산소 농도, 서식처 구조 등 수온 이외의 환경 요인에도 간접적으로 영향을 미칠 수 있으므로^{(Allan and Castillo, 2007)} 독립변수로 포함하였다. 이를 통해 고도의 생태적 효과와 수온 효과를 모두 반영할 수 있도록 하였다. 모델 학습은 파이썬(Python)의 주피터 노트북(Jupyter Notebook, 7.2.2) 환경에서 수행되었다. 학습 데이터셋은 2010–2021년 동안 수집된 줄날도래속 3종의 출현자료와 고도 자료를 통합하여 구성하였다. 종별로 출현 여부를 이진 반응 변수(1 또는 0)로 설정하였다. 전체 데이터셋을 사용하여 분류(Classification) 모델을 학습하였고, 하이퍼파라미터 설정 시, 기본적으로 사용되는 의사결정나무의 개수는 100개로 고정하여 각 의사결정나무의 최대 깊이는 제한하지 않은 상태에서 구성되었다.

모델의 일반화 성능을 향상하고 최적 구조를 도출하기 위하여, 하이퍼파라미터 튜닝을 5겹 교차검증(5-fold Cross-validation)을 통해 성능을 비교하였다. 그 결과, 모든 종에서 공통적인 하이퍼파라미터 조합이 최적값으로 도출되었다. 개별 의사결정나무의 최대 깊이를 20으로 제한하고 하나의 내부 노드가 분할되기 위해 필요한 최소 샘플 수를 5로 설정하였다. 또한 말단(Leaf) 노드에 최소 2개의 샘플이 존재해야 분할을 허용하도록 설정, 각 노드에서 분할 기준을 찾을 때 전체 입력변수 중 제곱근에 해당하는 수만큼의 변수들을 무작위로 선택하여 평가하였다. 또한, 모형의 예측 성능을 개선하고 과적합을 방지하기 위해 교차검증(Cross-validation)을 활용한 내부 검증을 수행하였다.

성능 평가는 Receiver operating characteristic – Area under the curve (ROC-AUC) 및 혼동 행렬(Confusion matrix)을 기준으로 수행하였으며, 모든 종에서 AUC가 0.90 이상으로, 양호한 수준의 예측력을 나타냈다. 학습이 완료된 랜덤 포레스트 모델은 한반도 격자 단위(Cell-by-cell)로 각 줄날도래속 종의 출현 확률을 산정하였고, 예측한 결과를 시각화하였다. 줄날도래속 3종에 대해 지도 학습(Supervised classification)을 위한 데이터셋을 구축하였으며, 고도 및 평형수온을 입력변수로 사용하였다. 이때 미출현 지점도 포함하여, 실제 자연 분포의 편향을 반영하고, 불균형 데이터 클래스 문제를 최소화하였다.

3.1 Current and Future Overlap between HSI Curves and Equilibrium Water Temperature Distributions

줄날도래속 3종의 PMF와 PDF를 기반으로 한 서식처 적합성 분석 결과, 열거한 종별 순서대로 Weibull, Generalized exponential, Inverse Weibull 모형이 가장 적합하였으며, 종별로 뚜렷한 열 생태 반응이 확인되었다(Fig. 2).

Fig. 2. Comparison of PMF–PDF fitted environmental response curves and predicted habitat suitability index (HSI) for three Hydropsyche species.

H. kozhantschikovi는 비대칭적으로 뚜렷한 정적편포(positively skewed distribution)를 보였으며 9–11℃ 구간에서 적합도가 가장 높았다. H. valvata는 서식 적합 온도 구간이 10–11℃였다. H. formosana는 상대적으로 낮은 출현도를 보였으며 서식 적합 온도 구간이 10℃ 내외로 온도 반응 구간도 좁았다.

3종(H. kozhantschikovi, H. valvata, H. formosana)의 HSI 곡선과 수온–면적 분포(PDF) 간의 중첩을 현재(파란색) 및 미래(붉은색) 조건에서 시각화한 결과(Fig. 3), 기후변화 시나리오에 따라 수온 분포가 고온 방향으로 이동함을 보여주었다. 종별 HSI와 면적 함수의 곱을 통한 중첩도는 3종 모두 현재 조건에서의 종별 중첩 HSI 값이 더 높았으며, 이는 수온 조건이 각 종의 적정 서식 온도 범위와 더 밀접하게 일치함을 시사한다. 반면, 미래(+2.45°C) 조건에서는 전체 면적 함수가 우측으로 이동하면서 각종의 최적 수온대에서 면적 비중이 감소하였고, 이에 따라 종별 중첩 HSI 값 또한 감소하였다(Fig. 3).

H. kozhantschikovi와 H. valvata는 미래 조건에서 온도 분포와의 겹침이 급격히 축소되었으며 적합성이 현저히 감소하였다. H. formosana는 현재 전체 면적 분포와의 겹침이 매우 제한적이었으며, 미래 조건에서도 적합 면적은 여전히 축소 경향을 나타냈다. H. kozhantschikovi의 중첩 HSI는 0.77 → 0.58, H. valvata는 0.77 → 0.56, H. formosana는 0.51 → 0.34로 감소하는 것으로 추정되었다.

Fig. 3. Changes in spatial overlap between habitat suitability index (HSI) and geographic distribution areas of three Hydropsyche species under current and the RCP 8.5 climate scenario (+2.45℃).

3.2 Current Distribution Modeling Based on Equilibrium Water Temperature

랜덤 포레스트 모델 성능은 ROC-AUC 기준으로 모두 0.90 이상으로 나타나 고도–수온 기반 이원적 입력만으로도 예측력이 양호한 것으로 평가되었다(Table 1). 이러한 결과는 고도와 수온이 줄날도래류의 서식처 적합도 결정에 중요한 변수임을 시사하며, 회귀 기반 평형수온 관계식이 공간적 수온 차이를 효과적으로 반영하고 있음을 뒷받침한다.

본 연구에서 종별 현재 분포 예측 결과는 다음과 같다. H. kozhantschikovi는 한반도 전역 중 북부와 중부의 중간 고도대에서 비교적 고르게 분포하였고, 동해안 북부와 서해안 북부의 일부 지역에서 예측 확률이 0.90 이상으로 나타났다. 따라서 특정 고도대에만 국한되지 않는 중온성(Mesothermal) 종의 분포 양상을 보이면서도 비교적 고고도에도 분포하는 경향을 보였다.

H. valvata는 전반적으로 낮은 서식 확률을 보였으나, 국지적으로 중부 내륙과 일부 남부 산지 지역에서만 0.40 내외의 제한적인 서식 가능성이 확인되었다. 따라서 중온성 종의 전형적인 분포 양상을 보이고 서식처 폭이 좁은 생태적 특성을 갖는 것으로 나타났다. H. formosana는 가장 낮은 서식 확률을 보여, 전국적으로 0.20 이하의 낮은 분포 가능성이 예측되었다. 일부 중남부 산지 지역에서만 희소하게 출현 가능성이 확인되었으며, 수온이 높은 지역에서 주로 서식하는 호온성(Thermophilic) 특성^{(Kim et al., 2024)}이 뚜렷하게 나타났다(Fig. 4).

Fig. 4. Spatial distribution of predicted habitat suitability (probability), equilibrium water temperature (Teq, ℃), and elevation (m) for three Hydropsyche species in South Korea under current and the RCP 8.5 climate scenario (+2.45℃)

고도 구간별 출현 확률 분석에서도 이러한 생태적 경향이 확인되었다. H. kozhantschikovi의 경우, 저지대와 중고도(0–800 m)에서 비교적 높은 출현 확률을 보였으며, 800 m 이상 고도에서는 급격히 감소하는 경향을 나타냈다. 이는 본 종이 저-중고도의 하천 환경에 적합하게 분포함을 의미한다. H. valvata는 전반적으로 낮은 출현 확률이었지만 200–400 m에서 가장 큰 출현 확률을 가지며 중고도 하천에서도 적응이 가능한 생태적 지위를 가진다. H. formosana는 전 고도 구간에서 낮은 출현 확률을 보여, 한국 내 분포가 제한적임을 보여준다.

특히 H. formosana는 0–400 m의 저지대에서만 다소 높은 출현 가능성을 보였으며, 이는 제한된 환경 조건에서만 생존할 수 있는 특성을 가짐을 의미한다(Fig. 5). 이러한 결과는 확률 분포 모델을 사용한 고도에 따른 줄날도래속 3종의 최적 적응 범위(H. kozhantschikovi, 242–574 m; H. valvata, 174–485 m; H. formosana, 137–272 m)와^{(Kim et al., 2024)} 비슷한 경향을 보였다. 따라서, 본 예측 모델은 수온과 고도만을 주요 변수로 사용하였음에도 종별 고도 구배 및 서식 적합성 차이를 명확하게 재현하였다는 점에서, 고도와 수온 기반의 분포 예측 실효성을 제시한다.

Fig. 5. Altitudinal distribution shifts of three Hydropsyche species under current and the RCP 8.5 climate scenario (+2.45℃).

3.3 Future Distribution Modeling under Climate Warming (+2.45℃) Based on Equilibrium Water Temperature

2100년을 기준으로 RCP 8.5 시나리오에 따른 기후변화 영향을 평가하기 위해, ^{Ahn and Han (2010)}이 제시한 전국 평균 수온 상승량(+2.45 °C)을 평형수온에 적용하여 미래 서식처 적합도를 예측하였다. 다른 환경변수는 현재 상태를 그대로 유지하였으며, 모델 구조와 학습 자료, 하이퍼파라미터는 현재 시나리오와 동일하게 설정하여 직접적인 비교가 가능하게 하였다. 또한 예측 확률에 대한 서식처 적합도 임계값은 실제 서식처 예측을 이진화하기 위해 0.5로 선정하였다^{(Hosseini et al., 2025)}.

H. kozhantschikovi는 평균 예측 확률이 전반적으로 감소(Δ = −0.1121)하였으며, 적합지 면적은 74,837.2 km²에서 57,842.7 km²로 16,994.5 km² 감소하였다(Table 2). 이진 변화를 보면 33,102.6 km²의 소실(1→0), 16,108.1 km²의 신규 출현(0→1)이 나타나 순감소를 보였다(Table 3). 미래에는 저고도(0–200 m) 및 일부 고고도(>1,600 m) 지역에서 적합성이 낮아져 전반적 축소를 나타냈다(Fig. 5). 위도별 분포에서는 36–38° N 부근의 고위도 구간에서 상대적으로 높은 적합성이 유지되어, 서식처가 북쪽으로 미세하게 치우치는 경향을 확인할 수 있었다(Fig. 6).

H. valvata는 평균 예측 확률 변화는 미미했으나(Δ = −0.0066), 적합지 면적은 21,757.5 km²에서 8,876.5 km²로 크게 줄어들었다(Δ = −12,881 km²)(Table 2). 이진 분석에서는 소실(18,560.6 km²)이 신규 출현(5,679.5 km²)보다 훨씬 많았다(Table 3). 고도대별로는 현재보다 800–1,200 m 구간에서 확률이 낮아져 적합지가 급격히 축소되는 경향이 확인되었다(Fig. 5). 위도별 예측 확률은 일부 구간에서 소폭 증가했으나 대체로 0.50 임계값을 넘지 못해 적합지 면적이 줄어든 것으로 나타났다(Fig. 6).

Fig. 6. Latitudinal distribution shifts of three Hydropsyche species under current and the RCP 8.5 climate scenario (+2.45℃).

H. formosana는 평균 예측 확률은 약간 증가(Δ = +0.0196)했으나, 적합지 면적은 5,252.9 km²에서 1,169.2 km²로 4,083.7 km² 감소하였다(Table 2). 이때 소실(5,193.7 km²)이 신규 출현(1,110 km²)보다 많았다(Table 3). 고도대별 분석에서도 대부분의 구간에서 적합 확률이 현저히 낮아져, 특히 저고도(0–600 m) 구간에서 거의 소실되는 결과를 보였다(Fig. 5). 위도별 변화는 일부 구간에서 확률이 소폭 상승했지만 역시 임계값을 넘지 못해 적합지가 급격히 축소되었다(Fig. 6). H. kozhantschikovi는 전반적으로 가장 큰 폭의 감소와 불안정한 분포를 보였으며, 면적 기준으로 약 45.6% 지역에서 0.20 이하의 감소가 나타났다. H. valvata와 H. formosana는 평균 확률 변화는 크지 않았으나 현재 분포가 임계값 0.5 부근에 집중되어 있어, 소규모 확률 변화가 적합지 소실로 이어졌다(Table 4).

예측 모델에 따른 기후변화(+2.45 ℃)는 수온에 민감한 줄날도래속의 저위도 및 저고도의 적합 서식처를 압축시키는 효과가 나타났다(Fig. 5; Fig. 6). H. kozhantschikovi, H. valvata, H. formosana는 모두 서식 면적이 감소하였으며, 특히 H. kozhantschikovi는 수온 상승에 따른 고고도로의 이동이 확인되었다. 서식 면적은 크게 감소되었으나, 기후변화 전후 적합지 면적은 가장 넓게 나타나, 고도보다는 수온이 더 밀접한 분포 제한 요인임을 나타냈다. 서식범위는 가장 넓게 유지하여 단기적으로는 취약종으로 보기 어려우나, 지속적인 온난화가 이어질 경우 한계수온을 넘어 급격한 서식처 붕괴가 발생할 가능성이 있으며, 잠재적 취약종으로 장기 모니터링과 조기경보 체계 구축이 필요하다.

H. valvata는 고도와 위도에 대한 분포 변화가 미미했으며, 평균 예측 확률이 적게 감소하였고 적합지 면적이 큰 감소를 하는 종으로 나타났다. 이는 최적 적응 범위에 대해 수온과 고도의 환경 요인에 크게 민감한 종임을 나타내며 기후변화 취약종으로 판단된다. H. valvata를 보전하기 위해서는 최적 서식 조건을 갖추는 것 중요하며, 고도와 위도에 대한 분포 변화가 미미하다는 점에서 현재의 서식처를 보전하는 것 또한 중요하다.

H. formosana는 위도에 대한 이동이 적었고, 고도에 대한 이동에서 고고도로의 이동이 예측되었으며, 평균 예측 확률이 증가했음에도 적합지 면적은 감소하여 H. formosana의 분포 제한 요인 중에서 고도보다는 수온의 영향이 크다는 것을 나타낸다. 또한 평균 출현 확률이 증가했지만 적합지 면적이 감소하여 종 보존 측면에서는 기후변화 취약종으로 판단될 수 있다. H. formosana를 보전하기 위해서는 제한적인 분포 특성상 서식처 연결성을 유지하고, 고립된 개체군을 보호하기 위한 서식처 복원, 서식처 네트워크 조성 등의 국지적 서식처 관리가 필요하다.

결론적으로 중-고위도 지역의 잠재적 피난처(refugia) 확보와 고지대 서식처 보호가 필요하며, 이에 따라 향후 연구에서 수온 변화와 더불어 수문 및 수질 변화 시나리오를 통합적으로 고려하여 종 간 상호작용, 기후-하천 수문 연계 모델링 등 생태적 복잡성을 반영한 통합적 예측 체계로의 발전이 요구된다.

3.4 Study Limitations and Future Improvements

본 연구는 남한 면적을 범위로 고도와 평형수온 관계식에 따른 평형수온을 입력변수로 사용한 모델에 기반하고 있다. 고도는 유속, 산소농도, 서식처 구조 등 환경 요인을 간접적으로 반영할 수 있다^{(Allan and Castillo, 2007)}. 하지만 다양한 환경 요인들을 모두 반영할 수 없다는 한계가 있다. 평형수온은 고도와 위도 및 해안선으로부터의 직선거리를 설명변수로 하여 수온을 추정할 수 있다. 하지만 모든 지역적 특수성을 반영하기 어렵다는 한계가 있다. 또한 미래 분포 예측에서 우리나라 해안 부근이 인근 내륙보다 높은 예측 확률이 나타났다. 이러한 결과는 고도와 평형수온 관계식 기반 평형수온을 주요 입력변수로 두 가지 입력변수만을 사용한 모델의 한계점을 나타내며, 생물 자료의 시공간적 불균형, 지점 간 조사 강도의 차이가 예측에 영향을 미칠 수 있다는 가능성을 보여준다. 따라서 본 결과는 최소한의 기후 영향만을 반영한 전망이며, 실제 미래 분포 변화는 서식처 파편화, 오염, 인간 활동 등 다른 요소들과의 상호작용에 따라 결과가 달라질 가능성이 있다.