2.1 연구 대상 지역 및 수리지질 특성

연구의 공간적 범위는 한반도 남동단에 위치한 부산광역시(동경 128°45’∼129°18’, 북위 34°52’∼35°23’)로 설정하였다. 이 지역은 복잡한 리아스식 해안선이 발달하여 내륙 지하수계가 해수면 조석 변동에 즉각적으로 반응하는 독특한 수리지질학적 환경을 갖추고 있으며, 이는 해수 침투에 대한 높은 취약성으로 이어진다^{(Lee and Kim, 2015)}. 기후적으로는 해양성과 몬순 기후가 공존하는데, 연평균 강수량(1,466.2 mm)의 절반 이상이 여름철에 집중되는 강한 계절적 편중성은 지하수위 변동성을 증폭시키는 핵심 요인으로 작용한다. 지질학적 측면에서도 내륙의 균열 암반층과 해안가의 다공성 충적층이 혼재되어 있어 강우 및 조석 등 외부 스트레스에 대한 수리적 반응이 구역마다 상이하게 나타나는 비선형적 거동 특성을 보인다. 특히 해안가에 발달한 다공성 충적층은 주로 투수성이 우수한 조립질 모래 및 자갈로 구성되어 있어 높은 유효공극률과 수리전도도를 지닌다. 이러한 미고결 퇴적층 구조는 조석 변동에 따른 해수의 수평적 유동과 수리적 분산을 물리적으로 촉진한다. 결과적으로 이는 유입된 해수가 내륙측 기반암의 단열대를 따라 내륙 깊숙이 침투할 수 있는 우선적 유동 경로를 제공하여, 해수 침투에 매우 취약한 수리지질학적 환경을 조성한다^{(Cleary et al., 2025)}. 본 연구에서 지하수위 예측 및 해석을 위해 선정한 주요 관측 지점의 지리적 분포는 다음과 같다(Fig. 1).

2.2 연구 방법론 및 모델링 프레임워크

본 연구에서는 지하수위 예측의 신뢰성을 확보하고 모델의 예측 정밀도를 극대화하기 위해, (1) 데이터 전처리(data preprocessing), (2) 변수 생성 및 선택을 위한 특성 공학(feature engineering), (3) 머신러닝 모델링을 통한 지하수위 예측(groundwater level prediction)으로 구성된 3단계 연구 프레임워크를 수립하였다. 각 단계는 수집된 시계열 데이터의 이상치 정제부터 예측 성능에 기여하는 최적 변수 추출, 그리고 하이퍼파라미터 최적화를 통한 모델 구축 및 성능 평가에 이르는 일련의 과정을 체계적으로 포함하고 있으며, 이에 관한 전체적인 연구 절차 및 흐름은 다음과 같다(Fig. 2).

Fig. 2. Conceptual framework of the proposed methodology for groundwater level prediction.

2.3 트리 기반 모델의 특징

지하수위의 복잡한 비선형 변동성을 예측하기 위해 단일 결정 트리(Decision Tree)와 5종의 앙상블 기법(XGBoost, Random Forest, GBT, LightGBM, AdaBoost)을 적용하였다. 이들 모델은 변수 간의 복잡한 상호작용 학습에 탁월하며, 설명 가능한 인공지능 기법인 SHAP와 결합하여 블랙박스 모델의 한계를 극복하고 예측 정확도와 해석 가능성을 동시에 확보할 수 있다^{(Lundberg and Lee, 2017)}.

2.4 모델 검증

트리 기반 머신러닝 모델의 지하수위 예측 성능을 객관적으로 평가하기 위해 평균 제곱근 오차(root mean square error, RMSE), 평균 절대 오차(mean absolute error, MAE), 내쉬-서클리프 효율지수(Nash-Sutcliffe efficiency, NSE), 그리고 결정계수(coefficient of determination, R²)의 네 가지 통계 지표를 활용하였다.

RMSE는 예측값과 관측값 차이의 제곱 평균에 대한 제곱근으로, 큰 오차에 가중치를 부여하므로 피크값 등의 이상치에 민감하게 반응하는 특성이 있으며, 값이 0에 가까울수록 예측 오차가 작고 모델이 관측값을 정밀하게 재현함을 의미한다. MAE는 예측값과 관측값 간 절대 오차의 평균으로 모든 오차를 동일한 가중치로 취급하며, 예측값이 관측값으로부터 평균적으로 얼마나 벗어나는지를 직관적으로 보여준다. NSE는 수문 및 수리 모델링, 수자원 시계열 예측 분야에서 가장 널리 사용되는 대표적인 모의 효율 평가 지표로, 모델의 잔차 분산을 관측값 분산과 비교하여 산출한다. 이 값이 1에 가까울수록 관측 평균을 기준으로 한 예측보다 뛰어난 성능을 나타내며, 0은 관측 평균과 동등한 수준, 음수는 평균 예측보다 성능이 낮음을 의미한다. 마지막으로 R²는 선형 회귀 및 머신러닝 모델의 성능 평가에서 광범위하게 사용되는 지표로, 모델의 독립 변수에 의해 설명되는 종속 변수의 총 분산 설명력을 정량화하며 1에 가까울수록 설명력이 높음을 뜻한다. 각 평가 지표를 산출하기 위한 수학적 계산식은 다음과 같다:

여기서 $n$은 데이터 포인트의 총 개수, $p_i$는 예측값, $o_i$는 관측값을 나타내며, $\overline{p}$와 $\overline{o}$는 각각 예측값과 관측값의 평균을 의미한다.

2.5 테일러 다이어그램

하이퍼파라미터 최적화가 완료된 6종의 트리 기반 모델을 대상으로 예측 성능을 종합적으로 평가하기 위해 테일러 다이어그램(Taylor diagram)을 활용하였다. 테일러 다이어그램은 상관계수(R), 표준편차(STD), 중심화된 평균 제곱근 오차 거리(RMSD)를 단일 도표에 시각화하여 모델과 관측값 간의 유사성을 종합적으로 평가하는 기법이다. 이 도표는 세 가지 주요 통계량이 코사인 법칙에 의해 기하학적으로 상호 연관되어 있다는 점을 기반으로 구성되며, 각 요소의 통계적 의미는 다음과 같다:

ㆍ 상관계수(correlation coefficient, R): 다이어그램의 방위각(azimuthal angle)으로 표시되며, 모델의 예측값이 관측값의 변동 패턴을 얼마나 유사하게 모사하는지를 나타낸다. 상관계수가 1에 가까울수록 모델과 관측값의 시계열적 패턴이 일치함을 의미한다.

ㆍ 표준편차(standard deviation, STD): 원점으로부터의 방사형 거리(radial distance)로 표현되며, 모델이 관측 데이터의 변동 강도(variability)를 얼마나 적절히 재현하는지를 나타낸다. 모델의 표준편차가 관측값의 표준편차 지점과 일치할수록 데이터의 변동 폭을 정확하게 모의한 것으로 간주한다.

ㆍ 중심화된 평균 제곱근 오차 거리(centered root-mean-square difference, RMSD): x축의 기준점(observation point)으로부터의 직선거리로 나타나며, 모델과 관측값 간의 정량적인 패턴 차이를 의미한다. 기준점으로부터의 거리가 짧아 RMSD가 낮게 산출될수록 모델의 예측 결과가 관측값에 더 근접함을 나타낸다.

2.6 SHAP 분석

머신러닝 모델의 고질적인 문제인 블랙박스 특성을 극복하고 예측 결과에 대한 해석 가능성을 확보하기 위해 SHAP 분석을 수행하였다. SHAP는 게임 이론(game theory)을 기반으로 각 입력 변수가 개별 예측 결과에 기여하는 정도를 샤플리 값(Shapley value)으로 정량화하여 산출하는 방법이다. 본 분석에서는 트리 기반 모델의 연산 구조에 최적화된 TreeSHAP 알고리즘을 적용하여 변수별 기여도를 평가하였다^{(Lundberg et al., 2020)}.

학습이 완료된 모델의 테스트 데이터셋을 대상으로 SHAP 값을 계산하였으며, 각 변수의 평균 절대 SHAP 값(mean absolute SHAP value)을 도출하여 전역적 중요도(global importance)를 산출하였다. 이러한 분석 과정은 모델 예측에 지배적으로 작용하는 수문 인자를 규명함으로써 예측 결과의 수문학적 타당성을 검증하고 결과 해석의 신뢰성을 제고하는 데 기여한다. 도출된 SHAP 분석 결과는 변수별 영향력을 시각화하여 결과 및 고찰 단계에서 심층적으로 논의하였다.

3.1 하이브리드 특성 선택 결과 및 해석

해안 지하수위 예측에 핵심적인 변수를 식별하기 위해, 피어슨 상관계수 기반의 필터 방식과 재귀적 특성 제거 기반의 래퍼 방식을 결합한 하이브리드 특성 선택 기법을 적용하였다.

먼저 1단계인 피어슨 상관계수 분석으로 도출한 다변량 상관 행렬을 바탕으로, 기상 및 수문 변수 간의 상호 의존성을 파악하였다(Fig. 3). 분석 결과, 일 강수량과 이동평균 강수량 변수들(MA_d10, MA_d20, MA_d30) 간에 강한 다중공선성이 확인되었으며, 두 가지 표준강수지수(SPI 6, SPI 12) 역시 통계적으로 유사한 수문학적 정보를 내포하고 있는 것으로 나타났다. 이에 따라 중복 변수 그룹 내에서 종속 변수인 지하수위와 가장 높은 상관관계를 보인 30일 누적 강수 지표(MA_d30)와 12개월 표준강수지수(SPI 12)를 최종 유효 변수로 채택하여 다중공선성 위험을 배제하였다.

단기 강수량보다 30일 누적 강수 지표 및 12개월 표준강수지수의 상관성이 더 높게 나타난 것은 연구 대상지의 수리지질학적 특성과 고도의 도시화가 복합적으로 작용한 결과이다. 본래 기반암의 균열 대수층은 강수에 대해 신속한 수리적 반응을 보이는 것이 일반적이지만, 대상 지역은 광범위한 불투수성 포장으로 인해 강수의 직접적인 수직 함양이 차단되어 지표수의 대수층 도달 시간(지체 시간)이 크게 연장된 것으로 분석된다^{(Shuster et al., 2005)}. 실제로 기후에너지환경부 환경공간정보서비스 자료에 따르면, 부산광역시의 불투수면 비율은 29.56%(불투수면 면적: 232.23 km², 행정구역 면적 785.50 km²)로 전국 특별⋅광역시 및 자치도 중에서 서울특별시(53.16%) 다음으로 높다. 이러한 환경은 선행 강우의 영향이 누적된 장기 수문 신호가 지하수위 변동을 지배하게 만드는 주요 원인으로 작용한 것으로 판단된다.

한편, 모든 변수 중 지연 지하수위(GWL(t-1)가 현재 지하수위와 가장 완벽한 수준의 양의 상관관계(r = 1.00)를 보였는데, 이어지는 래퍼 과정에서도 6개 모델 모두에서 지연 지하수위가 최우선 예측 변수로 공통 도출되었다(Table 3). 이는 수문 시계열 데이터가 지니는 강한 자기상관성과 수리적 관성에 따른 시간적 지속성을 모델이 타당하게 학습했음을 뒷받침한다.

또한, 전기전도도 역시 그 뒤를 이어 높은 상관관계(r = 0.84)를 나타냈으며, 대다수 모델에서 두 번째로 중요한 특성으로 산출된 것은 연구 대상지의 복잡한 리아스식 해안선 환경에서 염수 침투가 미치는 지배적인 영향력을 반영한다^{(Shim et al., 2002; Werner et al., 2013)}. 부산 해안 대수층은 복잡한 해안선과 다양한 투수성 매질이 혼재되어 있어 염수가 유입되기 쉬운 우선적 경로(preferential pathway)가 형성되기 쉬운 구조이며^{(Lee and Kim, 2015)}, 본 분석 결과는 이러한 지역 고유의 수리지질학적 특성을 명확하게 대변한다. 한편, AdaBoost 모델의 경우 알고리즘 내부 메커니즘에 의해 지연 지하수위와 전기전도도를 제외한 나머지 변수들의 기여도가 0으로 산출됨에 따라, 최종적으로 단 2개의 핵심 변수만이 예측 인자로 채택되는 특징을 보였다.

Fig. 3. Pearson correlation matrix of groundwater levels and various hydrometeorological variables.

3.2 모델 성능 비교 및 하이퍼파라미터 튜닝 결과

베이지안 최적화 적용 전후의 6종 트리 기반 모델에 대한 예측 성능(RMSE, MAE, NSE) 변화를 정량적으로 비교 분석하였다(Table 4). 튜닝 전 초기 상태에서는 부스팅 기반 모델들이 전반적으로 우수한 예측력을 보인 반면, Random Forest와 Decision Tree는 상대적으로 저조한 성능을 나타냈다. 구체적으로 AdaBoost는 RMSE 0.198, NSE 0.892로 초기 성능이 가장 우수하였으나, Random Forest와 Decision Tree는 NSE가 각각 0.564, 0.787에 그쳐 앙상블 기법 간의 뚜렷한 예측 성능 편차를 보였다.

하이퍼파라미터 튜닝 이후 대다수 모델의 성능이 향상되었으나, 각 알고리즘의 구조적 특성에 따라 개선 폭에는 유의미한 차이가 확인되었다. 특히 부스팅 기반 모델군(XGBoost, GBT, LightGBM, AdaBoost)에서 현저한 성능 향상이 두드러졌다. 가장 뛰어난 결과를 도출한 AdaBoost는 튜닝 후 RMSE 0.040, MAE 0.032, NSE 0.996을 기록하며 최적화 전 대비 NSE가 0.104 상승하는 높은 효율을 입증하였다. LightGBM 역시 NSE 0.992의 매우 정밀한 예측 성능을 달성하였으며, XGBoost와 GBT 또한 각각 0.985, 0.984의 우수한 NSE 수치를 기록하여 해안 지하수위 예측에 있어 부스팅 기법의 강세를 뒷받침하였다.

반면, Random Forest는 튜닝을 통해 NSE가 0.332 급증한 0.896에 도달하며 가장 큰 폭의 성능 개선을 이루었으나, 절대적인 예측력 측면에서는 여전히 부스팅 모델군에 미치지 못하였다. 단일 트리 모델인 Decision Tree는 튜닝 전후로 지표의 변화가 전혀 발생하지 않아, 해당 알고리즘 구조에서는 하이퍼파라미터 최적화를 통한 모델 성능 개선에 명확한 한계가 있음을 확인하였다.

이러한 최적화 효율의 차이는 각 알고리즘이 지닌 고유한 구조적 특성에 기인한다. 단일 모델인 Decision Tree는 단순한 이진 분할(binary split) 구조에 의존하므로, 복잡한 하이퍼파라미터 조합에 대한 성능 개선 반응성이 제한적이다. 이는 알고리즘 특성상 매 노드에서 국소적 최적해(local optimum)만을 탐색하며, 축에 평행한 계단형 결정 경계만을 생성할 수 있는 구조적 한계를 지니기 때문이다^{(Hastie et al., 2009)}. 즉, 베이지안 최적화를 통해 파라미터 탐색 공간을 확장하더라도 단일 트리 자체의 제한된 표현력으로 인해 특정 임계점 이후에는 수리적 성능 포화 상태에 도달하게 되며, 이로 인해 튜닝 전후의 유의미한 지표 변화가 나타나지 않을 수 있다. 결과적으로 단일 결정 트리는 파라미터 제어만으로는 독립적인 트리 구조 자체의 분할 한계를 극복하기 어렵다는 취약성을 지니며, 본 연구의 실험 결과(Table 4)는 이러한 이론적 한계를 잘 실증하고 있다. 또한, Random Forest는 개별 트리를 독립적으로 구축하는 병렬 앙상블 방식을 사용하기 때문에, 튜닝을 통해 일부 개별 트리의 성능이 개선되더라도 최종 앙상블 평균화 과정에서 그 효과가 상쇄되어 전체적인 성능 향상폭이 제한될 수 있다^{(Yang and Shami, 2020)}.

반면, 부스팅 계열 모델들은 각 학습 단계에서 이전 트리의 오차를 보정해 나가는 순차적 메커니즘을 취하므로, 하이퍼파라미터 최적화를 통해 이러한 잔차 보정 능력을 극대화할 수 있다. 특히 LightGBM의 리프 중심 분할(leaf-wise split) 전략과 GBT 등 기존 부스팅 기법들의 수평 중심 분할(level-wise split) 전략 등 알고리즘별로 특화된 트리 성장 방식이 베이지안 최적화의 탐색 과정과 성공적으로 맞물리면서 현저한 성능 향상을 이끈 것으로 판단된다.

실제 지하수위 관측값과 하이퍼파라미터 튜닝 전후의 모델 예측값을 비교한 시계열 거동을 살펴보면, 최적화 과정을 거친 이후 모든 앙상블 모델의 예측 궤적이 관측 수위의 변동 패턴을 훨씬 더 정밀하게 재현함을 확인할 수 있다(Fig. 4). 이러한 최적화에 따른 예측 모델의 성능 개선 효과는 관측값과 예측값 간의 산점도에 나타난 군집의 밀집도 및 결정계수의 뚜렷한 향상 결과를 통해서도 교차 검증된다(Fig. 5).

Fig. 4. Time-series comparison of observed and predicted groundwater levels (before and after optimization) using different ensemble learning models.

Fig. 5. Scatter plots of observed versus predicted groundwater levels for the six machine learning models: (a) XGBoost, (b) Random Forest, (c) Gradient Boosted Trees, (d) Decision Tree, (e) LightGBM, and (f) AdaBoost.

Fig. 6. Taylor diagram displaying the statistical comparison of predictive performance among the six machine learning models in terms of correlation coefficient, standard deviation, and root mean square error.

3.3 테일러 다이어그램을 통한 종합 성능 평가

테일러 다이어그램을 활용하여 관측값 대비 각 예측 모델의 상관계수(R), 표준편차(STD), 그리고 중심화된 평균 제곱근 오차 거리(RMSD)를 시각적으로 종합 비교하였다(Fig. 6). 하이퍼파라미터 최적화를 거친 앙상블 모델들은 전반적으로 관측 기준점 근처에 밀집된 군집을 형성하였다. 특히 AdaBoost와 LightGBM은 관측점에 가장 근접하게 위치하여, 0.99 이상의 높은 상관계수와 극히 낮은 RMSD를 동시에 달성하는 정밀한 예측 성능을 보였다.

부스팅 계열인 XGBoost와 GBT 역시 우수한 상관계수를 기록하였으나 앞선 두 모델 대비 상대적으로 높은 RMSD를 나타내어, 극단적인 수위 변동 폭을 완벽하게 모사하는 데에는 미세한 구조적 한계가 있음을 시사한다. Random Forest는 양호한 상관관계를 보였으나 RMSD 지표가 두드러지게 증가하였는데, 이는 개별 트리의 결과를 평균화하는 과정에서 발생하는 예측값의 평활화(smoothing) 경향이 반영된 결과이다. 한편, 튜닝을 통한 성능 개선이 없었던 단일 Decision Tree 모델은 도표상 관측점으로부터 가장 멀리 떨어져 있어 예측 효율이 가장 저조한 것으로 확인되었다. 결과적으로 본 시각적 분석은 트리 기반 앙상블 기법 중에서도 순차적 오차 보정 메커니즘을 채택한 부스팅 구조가 복잡한 해안 지하수위 거동 예측에 가장 적합함을 뒷받침한다.

그러나 이러한 부스팅 모델들의 우수한 통계적 지표는 모델이 지하수위 시계열 특유의 강한 시간적 지속성을 효과적으로 포착했음을 의미할 뿐, 외생 기상 변수를 활용한 실질적인 예측 능력이 온전히 확보되었다고 단정하기는 어렵다. 특히 모델들이 관측 데이터의 변동 강도를 정밀하게 재현한 것은 대수층 시스템 자체의 높은 자기상관성이 모델에 투영된 결과로 해석하는 것이 타당하다.

3.4 SHAP 분석을 통한 변수 중요도 및 기여도 해석

SHAP 분석을 수행하여 하이퍼파라미터 최적화 전후 트리 기반 모델들의 변수별 전역적 기여도와 개별 예측치에 대한 영향력을 정량적으로 평가하였다(Fig. 7-8). SHAP 요약 도표(summary beeswarm plot)는 각 데이터 포인트의 SHAP 값을 수평축 상의 밀집 분포로 시각화하며, 기준점(0)을 중심으로 양(+)의 영역은 지하수위 예측값의 상승을, 음(-)의 영역은 하강을 유발함을 의미한다. 도표 상의 개별 점들은 하나의 데이터 인스턴스를 나타내며, 각 변수의 특성값(feature value) 크기에 따라 적색(높은 값)에서 청색(낮은 값)으로 구분된다.

튜닝 전 초기 모델의 분석 결과, 모든 모델에서 지연 지하수위가 가장 높은 중요도를 나타냈으며, 전기전도도가 그 뒤를 이었다(Fig. 7). 전기전도도가 기상 변수들에 비해 높은 기여도를 유지하는 것은 연구 대상 지역의 복잡한 수리지질학적 환경에 기인한다. 암반과 충적층이 얽힌 부산 해안 대수층은 해수 유입에 민감하게 반응하며, 이로 인한 대수층 내부의 물리⋅화학적 상태 변화는 단기적인 기상 노이즈가 배제된 채 시계열적 연속성을 띠고 발현된다. 결과적으로 전기전도도는 해수 침투를 포함한 대수층의 거시적 상태를 대변하는 지시자로서, 지하수위 변동을 모의하는 데 있어 단기 기상 요인보다 결정적인 정보로 작용하는 것으로 해석된다.

알고리즘별 변수 중요도의 서열은 전반적으로 유사한 경향을 나타냈으나, 세부 기여 방식에는 일부 차이가 존재하였다. 특히 Random Forest와 GBT의 경우, 지연 지하수위와 전기전도도에 과도하게 편중되지 않고 다양한 입력 변수들을 복합적으로 활용하여 예측을 수행하는 다변량 의존 특성이 확인되었다.

하이퍼파라미터 최적화 이후에는 6개 모델 모두에서 지연 지하수위를 제외한 나머지 입력 변수들의 SHAP 분산 폭이 0 부근으로 급격히 수렴하며 기여도가 상쇄되는 뚜렷한 패턴 변화가 확인되었다(Fig. 8). 이러한 변수 기여도의 변화는 하이퍼파라미터 튜닝에 따른 모델 성능 지표(NSE)의 향상 추세와 밀접하게 연계된다(Table 4). 부스팅 계열 모델들은 최적화 이후 일제히 모의 효율이 상승하였으며, 튜닝 전 기저 성능이 가장 낮았던 Random Forest 역시 가장 높은 폭의 지표 향상을 달성하였다. 반면 단일 Decision Tree 모델은 튜닝 전후 성능 변화가 전혀 나타나지 않았다.

이러한 결과는 하이퍼파라미터 탐색 알고리즘이 예측 오차를 최소화하는 과정에서, 비선형성이 강한 기상 인자의 복합적 연산보다는 수문 시계열에 내재된 강력한 자기상관성에 가중치를 집중시키는 방향으로 모델 구조를 재조정했음을 시사한다. 즉, 순차적 오차 보정 메커니즘을 갖춘 부스팅 모델들은 최적화 과정을 통해 잔여 오차를 정밀하게 조정하는 과정에서 과거 수위 데이터에 대한 예측 의존성을 극대화함으로써 통계적으로 최적화된 모의 능력을 확보한 것으로 판단된다.

Fig. 7. SHAP summary beeswarm plots for the six machine learning models before hyperparameter tuning.

Fig. 8. SHAP summary beeswarm plots for the six machine learning models after hyperparameter tuning.

3.5 모델의 수리학적 특성 고찰 및 향후 연구 방향

트리 기반 모델과 SHAP 분석을 통합 적용한 본 연구 프레임워크는 해안 대수층의 지하수위 예측 메커니즘을 성공적으로 규명하고 모델의 데이터 분석적 특성을 도출하였다. 본 모델들이 달성한 높은 예측 정밀도는 전일 지하수위 데이터에 내재된 강한 시계열적 선행 정보를 알고리즘이 주도적으로 학습한 결과로 분석된다. 이는 현재 시점의 관측 데이터에 기반하여 근미래의 수위 거동을 추적하는 실시간 대수층 관리 및 단기 모니터링 도구로서 본 프레임워크의 실무적 유용성이 높음을 실증한다. 반면, SHAP 분석에서 강수량이나 기온 등의 외생 기상 변수들이 상대적으로 낮은 전역적 중요도를 나타낸 것은 트리 기반 알고리즘의 통계적 분산 설명 메커니즘에 기인한 현상이다. 즉, 시간적 연속성과 수리적 관성이 매우 강한 지하수위 시계열의 특성으로 인해 단기적인 기상 변수의 기여도가 통계적 마스킹 현상(masking effect)에 의해 가려진 것일 뿐, 기상 인자 본연의 수문학적 가치나 인과성이 배제된 것으로 해석해서는 안 된다.

이러한 수리학적 특성과 알고리즘적 구조를 보완하여 모델을 더욱 고도화하기 위해서는 향후 다각적인 외생 모델 개발이 요구된다^{(Wunsch et al., 2018)}. 기상 변수가 해안 지하수위에 미치는 순수한 영향력과 수리적 반응 메커니즘을 명확히 격리⋅검증하기 위해, 입력 변수에서 과거 지하수위 데이터를 배제한 모델 설계가 필수적이다. 오직 외생 기상⋅해양 인자만을 입력원으로 사용하는 외생 모델을 구축함으로써 단기적 지속성 모사를 탈피하고, 실질적인 기상 예보 자료와 연계하여 미래의 가뭄이나 침수 재해를 선행 예측할 수 있는 사전 방재 시스템으로 확장되어야 한다. 이와 더불어 계절별 특성 및 강우 강도별 이벤트에 대한 세분화 분석을 병행함으로써, 해안 대수층 시스템이 지닌 시공간적 이질성을 정밀하게 모사하려는 연구가 요구된다.

아울러, 복잡한 리아스식 해안 지형이 지닌 특수한 수리 구조적 취약성을 일반화하기 위해 향후 다양한 지질적 환경을 대변하는 해안 관측소 데이터를 확보하여 모델의 공간적 확장성과 일반화 가능성을 넓힐 필요가 있다. 특히 대수층 내 염수 침투를 지시하는 전기전도도의 경우, 해수 침투 현상 외에도 조위 변동이나 계절적 염도 구배 등 대수층 내부의 복합적인 환경 변화를 집약적으로 반영하므로 이에 대한 수리학적 세분화 분석이 유용하게 작용할 것이다. 최종적으로는 장단기 메모리(LSTM) 등 시계열 예측에 특화된 순환 신경망 계열의 딥러닝 모델이나 MODFLOW와 같은 전통적 물리 기반 모델과의 체계적인 비교 연구를 수행함으로써, 데이터 기반 머신러닝 모델이 도출한 예측 결과의 수문학적 타당성과 정밀도를 지속적으로 검증하고 보완해 나가야 할 것이다.