Feasibility Test for Hydraulic Conductivity Characterization of Small Basin-Scale Aquifers Based on Geostatistical Evolution Strategy Using Naturally Imposed Hydraulic Stress

자연 수리자극을 이용한 소유역 규모 대수층 수리전도도 특성화: 지구통계 진화전략 역산해석 기법의 적용 가능성 시험

Abstract

In this study, the applicability of the geostatistical evolution strategy as an inverse analysis method of estimating hydraulic properties of small-scale basin was tested. The geostatistical evolution strategy is a type of data assimilation method that can effectively estimate aquifer hydraulic conductivity by combining a global optimization model of the evolution strategy and a local optimization model of the ensemble Kalman filtering. In the applicability test, the geometry, hydraulic boundary conditions, and the distribution of groundwater monitoring wells of Hanlim-Eup were employed. On the other hand, a synthetic hydraulic conductivity distribution was generated and used as the reference property for ease of estimation quality assessment. In the estimations, two different cases were tested where, in Case I, both groundwater levels and hydraulic conductivity measurements were assumed to be available, and only the groundwater levels were available, in Case II. In both cases, the reference and estimated hydraulic conductivity fields were found to show reasonable similarity, even though the prior information for estimation was not accurate. The ability to estimate hydraulic conductivity without accurate prior information suggests that this method can be used effectively to estimate mathematical properties in real-world cases, many of which little prior information is available for the aquifer conditions.

1. 서론

지하수자원은 일반적인 지질자원과는 달리 일부 지역에 편재되지 않으며, 시-공간적으로 동적인 특성을 가진다. 현행 수내지 수십 개의 수량·수질 관측정에 기초하여 한 지역의 지하수를 관리하는 방법론은 통계적 개념에서 표본을 통한 대표성을 얻고, 이에 기초하여 그 지역의 지하 수자원을 관리하고자 하는 개념이다. 그러나 관측정의 수와 위치가 한정적인 현실적 상황에 비추어 볼 때, 관측정측정에 기반한 지하수 관리만으로는 수량·수질의 대표성을 얻기 어렵고 긴급한 수량·수질 사고에 체계적으로 대응하기 어렵다는 문제가 있다. 또한, 장기적인 수자원 관리계획의 유효성을 관측정 위주의 지점별 관측에 입각하여 평가할 경우, 광역적 효과를 가늠하는 데 있어서의 부정확성 및 불확실성이 수반될 수 있다. 이러한 관측정 위주의 지하수 관리의 제한성을 극복하기 위하여, 지하수의 유동과 지하수가 잠재적으로 포함할 수 있는 오염물질 확산의 공간적 범위를 포괄하는 기본 단위인 소유역을 지하수 관리의 기본단위로 하는 개념적 전환은 이미 국내·외적으로 이루어지고 있는 상황이다(Yeh et al., 2008; Graymore et al., 2008). 지하수 자원 관리의 기본 단위를 관측정-중심에서 소유역-중심으로 전환할 경우 관리에의 효율성과 지속개발 가능성을 도모할 수 있으며, 궁극적으로 지하수 관리의 건전성을 확보할 수 있다. 그러나 아직까지 소유역-중심의 지하수 관리를 위한 구체적인 방법론은 확립되어 있지 않은 상황이다.

소유역을 기본단위로 하는 지하수 수량·수질 관리를 위해서는 유역하부 대수층의 공간적 연결성을 정의하는 지질매질의 수리적 물성분포 특성화가 우선되어야 한다. 대수층 매질의 공간적 분포는 지하수 유동 및 오염물질 거동 특성 및 영향범위를 규정한다. 일반적으로 지하수 수량·수질에 관련되어 다양한 규모의 사업을 수행함에 있어, 대수층 매질의 공간적 분포는 가장 핵심적인 정보이다. 예를 들어, 신규 지하수 양수정의 설치, 지하수 상류 보전구역 설정, 잠재오염 시설의 입지 결정, 오염 정화를 위한 최적 위치, 신규 관측정의 최적 위치 결정 등은 대수층 매질의 분포에 기반하여 이루어지는 것이 일반적이다. 전통적으로 대수층 매질의 수리특성을 얻기 위한 방법으로 양수 및 순간 주입 시험을 포함하는 대수성 시험이 이용된다. 그러나 대수성 시험결과 얻어지는 결과는 대수성 시험의 세트업과 대수층 매질의 특성, 그리고 미지의 경계조건에 의해 결정되는 지지부피(support volume)를 대표하는 평균적인 값이며, 해당 시험의 결과가 시험정을 포함하는 대수층 전반의 고유값으로 볼 수 없다. 지질의 분포와 수리물성의 분포가 상응할 것이라는 가정 하에 지질매질을 편의적으로 분류하고 이에 근거하여 수리물성을 구역화하는 수리층서법(hydrostratigraphy)은 비록 국내에서 다수 이용되고 있는 상황이다. 그러나, 수리층서법은 특정 지질(e.g., 퇴적암)에서 만 성립되며 타 결정질암에 대해서는 그 근거가 희박하다. 국내에서는 적용된 사례가 없지만, 대상지역의 대수층 매질 분포를 정밀하게 해석하기 위한 방법으로 수리적 토모그래피(hydraulic tomography) 기법이 적용될 수 있다(Kitanidis, 1995; Yeh et al., 1996). 그러나 토모그래피 기법은 많은 양수 또는 추적자 시험의 수행과 그 결과에 대한 관측이 시행되어야 하기 때문에 비용-효율적인 방법이라 하기 어렵다. 또한, 토모그래피 기법을 적용하여 얻어지는 특성화된 대수층 규모 역시 소유역 규모와 비교할 때 매우 협소하여 대수층 전반을 파악하기 위한 방법이라 하기 어렵다.

한편, 국내에서는 2000년대 이후 현재까지 전국적으로 높은 밀도의 지하수 관측망이 설치·운영되고 있으며 이를 통하여 장기간의 지하수 시계열자료가 얻어지고 있다. 국내 관측망은 관측 목적에 따라 각기 다른 정부부처에서 운영되고 있으며, 이에는 환경부의 국가지하수관측망과 지하수수질측정망, 농림축산식품부의 농촌지하수관리 관측망과 해수침투관측망, 그리고 지자체에서 운영하는 지역 지하수관측망 등이 포함된다. 또한, 지하수가 유일한 담수자원인 제주도의 경우 육지에 비하여 상대적으로 높은 밀도의 특정목적 관측정이 설치·운영되고 있다. 지하수 관측망에서는 수위 또는 지하수 수질항목에 대한 연속적인 관측이 이루어지고 있어, 주기성을 고려한 지하수 수량 및 수질의 베이스라인 자료로 활용될 수 있다. 그러나 추가적인 지하수 관측의 필요성 및 정당성을 확보하기 위해서는 관측되는 자료의 직접적이고 구체적인 활용성을 제시할 필요가 있다. 이에는 당면한 기후변화에 따른 지하수 영향 파악, 지하수 개발 유망지 선정, 지하수 오염 심화 및 광역화에 따른 오염확산 억제책 수립, 도시화 및 대규모 지표·지하 시설 건설에 따른 지하수 오염 취약성 평가, 범국가적 수자원 관리사업의 유효성 검증 등을 꼽을 수 있을 것이다. 그러나 지하수 관측 자료만으로는 전술한 문제들에 대한 직접적 방법론이 되기는 어려우며, 대수층 불균질성에 대한 규명을 포함한 다양한 해석 방법론이 확보되어야 한다.

본 연구에서는 소유역-중심의 지하수 자원 관리에 핵심정보로 활용될 수 있는 대수층 수리물성의 공간적 분포를 특성화하기 위하여 유역 내에서 이루어지고 있는 지하수 관측을 이용할 수 있는 방법론을 제시하고자 한다. 즉, 관 측정으로부터 측정된 지하수위와 지하수위 변동을 유도하는 자연적 수리자극인 강우-함양을 이용하여 대수층 매질 수리전도도의 공간적 분포를 특성화할 수 있다는 것을 가상의 사례를 통해 입증하는데 본 연구의 목표가 있다. 이에는 관측기간 내 인위적으로 가해진 수리자극이 존재하지 않는다는 것을 전제로 한다. 또한 연구지역 형태, 경계조건, 그리고 관측정 분포는 제주 한림유역을 기준으로 설정하였으나, 실제 입력된 지하수위 시계열 자료는 가상의 것을 이용하였다. 이는 실제 한림 유역 대수층의 수리전도도 분포 자료가 부재한 상황에서, 추정된 수리전도도 분포와 레퍼런스 수리전도도 분포의 유사성 비교를 용이하게 하기 위한 방편이다. 본 연구에서는 특성화 방법론으로 지구통계 진화전략 역산해석 기법(Park, 2020; 박은규, 2020)을 이용하였다. 현재까지 지구통계 진화전략 기법은매우 이상적인 기하(정사각형 혹은 직사각형)를 갖는 도메인에만 적용이 이루어졌으며 본 연구에서와 같이 실제지역의 기하와 경계조건을 이용한 사례는 없다. Park(2020)에 의해 최초 소개된 지구통계 진화전략 역산해석기법은 광역최적화 기법인 진화전략(Hansen and Ostermeier, 2001), 지역최적화 기법인 앙상블 칼만 필터(ensemble Kalman filter, Evensen, 2003), 그리고 가우시안 프로세스 회귀분석(Gaussian process regression)을 연계하고 연산을 효율화하는 다양한 알고리듬을 개발·추가하여 적은 연산자원으로 효율적으로 역산해석을 수행하여 대수층의 수리전도도 분포를 추정하는 기법이다. 본 논문에서는 연구내용 전달의 효율성 및 편의성을 도모하고 이론 소개에 할애되는 논문의 페이지를 최소화하기 위하여 상세 이론 및 수학적 전개 전부를 과감하게 생략하였다. 다만, 이용된 기법에 대하여 보다 상세한 이론적 배경에 관심 있는 독자는 Park(2020) 및 박은규(2020) 연구를 참조하길 추천한다.

2. 본론

2.1. 연구지역의 위치 및 개략적 수리 정보

지구통계 진화전략 기법의 실제 소규모 유역 적용 가능성을 시험하기 위하여 선정된 지역은 제주특별자치도 제주시 한림읍이다(Fig. 1). 비양도를 제외한 한림읍의 면적은 약 90 km² 이며 동-서로 약 15 km 그리고 남-북으로 14 km의 길이를 갖는다. 한림읍의 위치는 제주도 서북부이며 서쪽으로는 한경읍 그리고 동쪽으로는 애월읍과 접하고 북서측으로는 바다와 접해있다. 한림읍 내에는 제주특별자치도에 의해 관리되는 총 6개의 지하수위 관측소(JD한림1, JD한림2, JD협재, JM명월, JR금악1, 및 JR금악2)가 운영되고 있으며, 한림읍 경계 인근 2개의 지하수위 관측소(JI원종장1 및 JD판포2)가 분포한다(Fig. 1). 따라서, 본 연구에서는 총 8개(내부 6개 경계 2개)의 지하수위 관측정이 있는 것으로 가정하였다.

Fig. 1. Study area (Hanlim-eup, Jeju, South Korea) and location of groundwater monitoring wells.

한림읍 내 최상류에 위치하는 관측정은 JR금악2 관정으로 최근 10년내(2011.1.1.-2020.11.7.) 평균 수위는 약 70 m이다. 동일 기간 인근 지역에 위치하는 관측소(고산)의 연평균 강수량이 약 1235.7 mm인 점을 고려할 때 한림읍 하부 대수층의 평균적인 수리전도도는 매우 높을 것으로 예상된다. 정류 상황을 가정한 평균적인 수리전도도의 경우 정확한 지하수 함양율이 주어질 경우 간단한 해석학적 해를 통해 산정이 가능하지만 본 연구에서는 이 값을 평가하지 않았다. 대신, 본 연구에서는 가상으로 재현된 수리전도도 분포를 실제 수리전도도로 가정하고 지구통계 진화전략 기법의 적용을 통한 수리전도도 분포 추정을 실시하였다. 또한 제주특별자치도에 의해 구분된 한림읍을 포함하는 실제 유역의 형태와 한림읍의 형태는 일치하지 않으나 본 연구에서는 행정구역상의 한림읍 경계를 유역의 경계로 가정하였다. 따라서, 수치모델 구성 시 유역 경계를 따라 북동측 선분은 바다와 접하고 있으므로 0 m의 수두를 갖는 고정수두(specified head) 경계로 설정하였으며, 유역의 남서측 및 북동측 선분은 지하수 흐름방향과 평행하다는 가정 하에 0 m/day의 비유량을 갖는 고정유량(specified flux) 경계로 설정하였다. 지하수 함양 시나리오 역시 가상의 자료가 이용되었으며, 해당 지역과 유사한 연평균 강우량을 갖는 국내 관측소를 무작위로 선정하여 해당 일강우량 중 7-9월(90일) 자료를 이용하였으며 지하수 함양율은 30%로 가정하였다. 일 강우를 이용함에 따라, 지하수 유동 모델링에서의 시간-차분의 개수는 90개로 설정하였으며, 초기조건 구성을 위하여 연평균 강우량의 30%를 함양량으로 적용한 정류모사를 실시하고 그 이후 부정류 모사를 실시하였다.

2.2. 도메인 차분 및 가상 수리전도도 분포

수치모사 대상 모델 도메인은 16 km × 16 km의 면적을 가지며, 100 × 100 격자에 의해 차분되었다. 따라서, 총 셀의 개수는 10000이고 각 셀은 0.0256 km^₂ 의 면적을 갖는다. 총 셀의 개수 중 한림읍 영역에 해당하는 격자는 3434개이다. 본 연구에서는 지구통계 진화전략과 USGS의 지하수 유동 모사 코드인 MODFLOW-2005(Harbaugh, 2005)를 결합하였으며, 앞서의 정보들을 고려하여 모델링 경계조건(i.e., Basic Package의 IBOUND 값)을 구성하였다. 수직적으로 모델 도메인은 전 영역에서 동일하게 100m의 두께를 지닌 한 개의 층으로 이루어져 있다고 가정되었다. 따라서, 본 연구의 지하수 유동 모사는 2차원 모사이며 추정 대상은 엄밀하게는 투수량계수이다. 그러나, 만약 대수층 두께의 수직적인 분포자료가 가용할 경우, 이를 손쉽게 수치모델에 반영할 수 있으며 결과적으로 수리 전도도로 변환할 수 있다. 한림읍 지역 대수층은 전체에 걸쳐 자유면 대수층으로 가정되었다.

가정된 추정대상(이하 레퍼런스) 수리전도도 분포 제작을 위해서는 순차 가우시안 모사기법(sequential Gaussian simulation, SGSIM)이 이용되었으며, GSLIB(Deutsch and Journel, 1992) 소프트웨어가 이용되었다. 해당 모델은 속성값의 공간적인 분포가 가우시안 프로세스를 따른다는 가정하에 추계론적으로 순차적 값을 생성하며, 이에 크리깅(kriging) 기법이 이용된다. 공간적 상관거리의 경우 남-북 방향을 따라 8 km 그리고 동-서 방향을 따라 2 km로 설정하였으며, 비참조(unconditional) 방식의 재현을 실시하였다. 해당 소프트웨어에 의해 생성된 값들은 평균 0 그리고 분산 1의 정규분포에 가까우며, 일반적인 수리전도도 분포가 log-normal 분포를 따르는 특성을 고려하여 생성된 값에 지수함수를 취하여 수리전도도 값을 제작하였다. 여기서, 로그 변환된 수리전도도의 평균값은 2로 그리고 분산은 1.5로 조정하였으며, 이는 약 7.4 m/day 평균 수리전도도에 해당한다. 그러나 실제 모델링 영역(i.e.,애월읍 경계 내부)으로 한정하였을 경우 평균과 분산은 각각 1.89 및 1.89로 산정되었다. 모델링 영역 내 최소 및 최대 로그 수리전도도는 –2.2와 5.3이며, 이는 0.11 및 207.3 m/day 수리전도도에 해당한다. 비산출율과 비저유계수는 전체 대수층에 걸쳐 각각 0.002 및 2 × 10^-5 m^-1의 값을 가진다고 가정되었다.

2.3. 추정 입력인자

지구통계 진화전략 기법을 이용한 수리전도도 추정을 위해서는 핵심적 정보(core information)로 각 관측정에서 얻어진 지하수위 시계열 자료가 필요하다. 이를 위하여 위에서 언급된 정보를 포워드 모델인 MODFLOW-2005에 반영하여 각 관측정 위치에서의 값을 지하수위 시계열 정보로 이용하였다. 추가적인 정보(auxiliary information)로는 각 관측정에서 사전 조사를 통하여 얻어진 수리전도도 값을 이용할 수 있으며, 본 연구에서는 수리전도도 측정치가 있는 경우와 없는 경우 모두를 가정하여 추정을 실시하였다. Fig. 2는 지하수위 수문곡선을 보여주며, 각 수문곡선에 해당하는 관측정은 Fig. 1과 같다. 그림에서 살펴볼 수 있는 바와 같이, 해안 인근에 위치한 관측정의 경우 지하수위 및 수위의 변동폭은 수리적 경계조건에 의해 제한된다. 반면에 상류로 갈수록 지하수위 및 변동폭은 대체로 증가하는 패턴을 보여준다.

Fig. 2. Groundwater levels monitored in the monitoring wells in Fig. 1. The monitoring wells in the order of the upper left to the lower right are JDHanlim1, JDHyeopjae, JMMyeongwol, JDPanpo2, JRGeumak1, JRGeumak2, JIWonjongjang, and JDHanlim2, respectively

2.4. 추정 결과의 평가

지구통계 진화전략 기법에 의해 추정된 수리전도도의 공간적 분포를 평가하기 위해서는 레퍼런스 수리전도도 분포와의 비교가 필요하며, 크게 세 가지 방법이 이용되었다. 첫 번째는 레퍼런스 및 추정 수리전도도 분포의 격자간 평균 제곱근 오차(root mean squared error, RMSE)이며, 그 범위는 0에서 ∞이다. 두 번째는 레퍼런스 및 추정 수리전도도 분포의 피어슨 상관계수()이며, 그 범위는 0에서 1이다. 여기서, 첫 번째와 두 번째 방법은 레퍼런스와 추정 격자들 사이의 값들 간 유사성을평가하는 방법이다. 그러나 이 두 가지 방법으로는 수리전도도 구조의 공간적 유사성을 설명하기에는 어려움이 있다. 실제, 대수층 내에서 지하수의 유동 및 오염물질의 거동에 있어 연장된 형태의 고투수성 구조의 역할은 상당히 크기 때문에 전체 도메인에 걸친 구조적 유사성은 매우 중요하다고 할 수 있다. 따라서 본 연구에서는 세 번째 방법으로 구조적 유사성을 평가하는 지표로 이미지 분석 분야에서 다수 활용되고 있는 구조유사성 인덱스 척도(structural similarity index measure, SSIM)를 이용하였다. SSIM의 범위는 0에서 1이며, 1의 값에 가까울수록 구조적인 유사성이 높음을 의미한다. SSIM의 논리적 배경에 대하여 보다 관심있는 독자는 Wang et al. (2004) 논문의 참조를 추천한다. SSIM은 역산해석으로 추정된 대수층 매질의 정확성을 측정하는 척도로 Sánchez-León etal. (2020)에 의해 이용된 바 있다. 레퍼런스 지하수위 분포와 추정된 지하수위 분포를 비교하기 위해서도 마찬가지로 세 가지 척도(RMSE, p, 및 SSIM)가 이용되었다.

2.5. 컴퓨터 시스템

본 연구의 지구통계 진화전략을 이용한 소유역 규모 대수층 수리전도도 추정에서 이용된 전산시스템은 2.5 GHz의 28개 코어와 768 GB 메모리를 지닌 워크스테이션이다. 개발된 코드는 Mathworks사의 MATLAB 상에서 운용되며 동시에 다수의 MODFLOW-2005 포워드 모델 연산을 위하여 Parallel Toolbox가 이용되었다.

3. 결과 및 토의

3.1. Case I: 관측정 지하수위와 수리전도도 측정값이 추정에 반영된 경우

앞서 언급된 바와 같이, 본 연구에서는 각 관정에서 관측된 지하수위 시계열 자료와 수리전도도 측정값을 동시에 수리전도도 추정에 이용할 수 있다. 우선, 두 정보가 모두 이용된 추정 사례(Case I)에 대하여 언급하고 이후 다음 장에서 지하수위만 추정에 이용된 사례(Case II)를 언급하고자 한다.

일반적으로, 8개의 관측정만으로 약 90 km2 면적을 가진 한림읍 전체의 수리전도도 분포를 평가하는 것은 매우 어려운 일이며 제한조건의 부족으로 인하여 과소결정에 의한 비유일성 해가 도출될 수 있다. 이를 위해서는 다수의 가정이 이용되었으며, 일단 각 관정에서 측정된 지하수위와 수리전도도 측정값의 지지부피는 1개 셀에 해당하는 2.56 × 103 km3(160 m × 160 m × 100 m)이다. 많은 경우, 이와 같은 가정은 합당할 수 있지만, 관정의 위치가 셀의 크기 대비 상대적으로 소규모인 매우 높은 투수성 구조(예. 용암동굴)와 일치할 경우 위의 가정은 성립하지 않을 수 있다. 두 번째 주요 가정으로는, 한림읍 지역 전체를 통하여 공간적 강수량 분포 및 지하수 함양율이 동일하다는 가정이다. 한림읍과 같이 완만한 지형경사를 가지며 광역적인 도시화가 이루어지지 않은 지역에서는 이러한 가정이 성립될 수 있다. 그러나, 유역 내 지표구배가 큰 경우 강수량이 동일하지 않을 가능성이 있으며, 일부 광역적인 도시화가 진행된 경우 지표피복에 의해 동일한 함양율을 적용하기 어려울 수 있다. 세 번째 주요 가정으로는, 지하수위의 변동이 강수의 함양에 의해서만 발생한다는 것이다. 만약 인위적인 지하수위 변화가 개입될 경우 지하수위의 변화를 발생시킨 인위적 수리자극의 크기 및 스케쥴에 대한 정보가 필요하다. 본 연구에서는 강수가 집중되는 우기 동안에 한정하여 강수량을 적용하였으며, 이는 해당 기간 동안 지하수의 이용이 비교적 적을 것이라는 가정을 내포한다. 또한, 열거된 가정 외에 실용적 지하수 모사를 위해 적용되는 일반적인 가정들이 본 연구에서 마찬가지로 적용되었다.

지구통계 진화전략 기법을 이용한 추정에 있어, 총 진화 세대수는 50회 그리고 각 세대에서의 최대 포워드 모델링 수행 횟수는 10회로 제한하였다. 지하수위 및 로그 수리전도도 측정값의 오차는 각각 0.05 m 및 0.5이며, 상관거리는 남-북 및 동-서 방향으로 각각 0.24 및 0.08 km라 가정하였다. 이는 실제 레퍼런스 로그 수리전도도 분포 제작을 위하여 이용하였던 값들과 매우 다른 값이며, 예측 초기 공간적 상관성에 대한 정보가 부재하다는 실질적 상황을 가정한 것이다. 실제 90 km2 면적의 지역에서 8개 관측정의 측정만으로 공간적 상관거리를 얻는 것은 매우 어렵다. 총 8개 관측정에서 측정된 로그 수리전도도의 평균은 1.92로 이는 전체 도메인에 걸친 평균 로그수리전도도 1.89에 비하여 약간 큰 값이다. 8개 관측정로그 수리전도도 분산은 2.42로 산정되었으며, 이는 설정된 로그 수리전도도 분산 1.89에 비하여 28% 가량 과대평가 되었다.

Fig. 3은 레퍼런스 로그 수리전도도 및 지하수위(부정류모사 90일) 그리고 추정된 로그 수리전도도 및 지하수위를 보여준다. 육안으로 비교하였을 때 레퍼런스와 추정 이미지들은 상당한 유사성을 보여준다. 해당 추정을 위하여 총 소요된 연산시간은 33.4초이며, 이 중 22.1초가 지하수 유동을 예측하는 MODFLOW-2005 수행(총 225회)에 쓰였다. 레퍼런스와 추정 로그 수리전도도의 RMSE는 0.74였으며, 이는 설정된 측정 오차인 0.5에 비하여 약간 크다. 두 로그 수리전도도 격자간의  값은 0.76으로 격자간의 유사성이 매우 큼을 지시한다. 이미지의 유사성을 지시하는 SSIM 값의 경우 0.74의 값을 보여 구조적으로도 레퍼런스 및 추정 수리전도도 분포 사이의 유사성이 상당히 큼을 확인할 수 있다.

Fig. 3. Reference hydraulic conductivity and groundwater level (left hand side) compared to the corresponding hydraulic conductivity and groundwater level estimated by the geostatistical evolution strategy used in this study (right hand side). In the estimation, both the groundwater levels and the hydraulic conductivity measurements of all monitoring wells were used (Case I).

레퍼런스와 추정 지하수위의 RMSE는 4.73 m로 나타났으며, 해안에 가까울수록 오차가 크고 상류구배로 갈수록 오차가 작아지는 경향을 보였다(Fig. 4). 두 지하수위의 격자간 로우 값은 0.999로 격자간 유사성은 매우 크다. SSIM 값의 경우 1에 가까운 값(0.9997)의 값을 보여 레퍼런스 및 추정 지하수위 분포의 구조적 유사성이 거의 일치에 가까울 정도로 높은 것을 확인할 수 있다. Fig. 5 는 초기 및 50세대에 걸친 추정 이후 추정 불확실성의 변화를 보여준다. 초기 불확실성 분포는 가우시안 프로세스 회귀분석을 통해 얻어진 것이므로 크리깅 분산과 동일하다. 그러나 크리깅 분산은 직접정보인 관측정에서의 수리전도도 측정치만 이용할 수 있는 반면, 자료동화 기법을 이용할 경우 직접정보와 간접정보(i.e., 지하수위 시계 열)를 모두 이용하여 상대적 불확실성을 산정하므로 보다 의미 있다고 할 수 있다. Fig. 5의 값 분포는 상대적인 불확실성 정도를 보여주는 그림으로 절대적 불확실성을 의미하지는 않는다. 그림에서, 전반적으로 추정 불확실성은 감소하는 경향을 보여주지만, 일부 위치에서는 추정 이후에도 여전히 상대적으로 높은 불확실성이 유지되고 있음을 확인할 수 있다. 이와 같이 높은 불확실성이 유지되는 위치는 추가적인 관측정 설치가 필요한 위치라고 볼 수 있다. 만약 그림에서 상대 불확실성이 높은 위치들에 추가 관측정을 설치할 경우 전반적인 추정 불확실성은 크게 감소할 것이라 예상된다.

Fig. 4. Groundwater levels monitored in the monitoring wells and the corresponding groundwater level estimations by the geostatistical evolution strategy. The monitoring wells in the order of the upper left to the lower right are JDHanlim1, JDHyeopjae, JMMyeongwol, JDPanpo2, JRGeumak1, JRGeumak2, JIWonjongjang, and JDHanlim2, respectively

Fig. 5. Comparison of the initial and the final relative uncertainty assessed by the geostatistical evolution strategy. The final relative uncertainty distribution was obtained at the end of 50 generations of estimations.

전반적인 Case I 추정결과에 대한 관찰을 통하여, 비록 관측정의 수가 적지만 관측정마다 측정된 수리전도도가 있을 경우, 상당히 정확한 수리전도도 분포를 예측할 수 있음을 확인할 수 있다. 또한 추정에 이용된 공간적 상관거리가 레퍼런스의 상관거리와 매우 다름에도 불구하고 예측의 질이 매우 우수함을 확인하였다. 이는 전통적인 역산해석 기법들에서는 상관거리에 대한 사전 정보가 예측의 질에 민감한 영향(Mu et al., 2020)을 주는 반면, 광역최적화에 기반한 지구통계 진화전략의 경우 사전 정보들에 대한 의존성이 상대적으로 낮음을 의미한다. 실제 광범위한 조사가 이루어지지 않은 지역에 대하여 역산해석 기법을 적용하여 수리상수 분포를 예측할 경우, 사전정보는 부정확할 가능성이 매우 높다. 본 연구에서의 관찰은 지구통계 진화전략 기법이 사전 정보가 부재하거나 부정확한 상황에서 실용적으로 활용될 수 있을 수 있음을 지시한다.

3.2. Case II: 관측정 지하수위만 추정에 반영된 경우

앞서 Case I에서는 8개 관측정에서 지하수위 관측과 수리전도도 측정값이 모두 가용한 상태임을 가정하고 이를 추정에 반영하였다. 그러나 관측정에서의 지하수위 관측은 가용하지만 수리전도도 측정값이 부재한 경우가 보다 일반적이라고 할 수 있다. 따라서, Case II에서는 8개 관측정에서 90일간 관측된 지하수위만이 이용 가능한 정보라는 가정 하에 수리전도도 분포 추정을 실시하였다. 8개 관측정에서 측정된 수리전도도 값이 부재하다는 것을 제외하면, Case II의 모든 가정은 Case I과 동일하다. CaseII에서는 정보량의 축소를 고려하여 지역최적화(i.e., EnKF)에 대한 추정 가중치를 광역최적화와 동일하게 0.5로 부여하였다. 이전 Case I의 경우 광역최적화 및 지역 최적화에 주어진 가중치는 각각 0.9 및 0.1이었다. 이러한 가중치 조정은 일부 수리전도도 값이 반영되는 상황에서 추정의 자유도 감소 및 근소한 차이를 보이는 추정들만 나타남에 따라 발생할 수 있는 동종교배(inbreeding) 효과(Houtekamer and Mitchell, 1998) 및 추정결과의 발산을 배제하기 위함이었다. 추정에 이용할 수 있는 정보량의 감소는 추정 수렴속도의 감소를 가져오므로 Case II에서의 총 진화 세대의 수는 Case I에 비하여 2배 증가시킨 100을 이용하였다. 추정을 위한 사전정보로 평균 로그 수리전도도 값 0을 입력하여 추정을 실시하였으며, 이는 실제 레퍼런스에 이용된 평균 로그 수리전도도인 2에 비하여 상당히 낮은 값이다. 이와 같이 부정확한 사전정보를 입력한 이유는 지구통계 진화전략의 광역최적화 특성을 확인하기 위함이다.

Fig. 6은 앞서의 Fig. 3과 마찬가지로 레퍼런스와 추정로그 수리전도도 및 지하수위(부정류 모사 90일)를 보여준다. 육안으로 비교하였을 때 레퍼런스와 추정 이미지들은 좋은 유사성을 보여주지만 앞서 Case I의 추정에는 미치지 못함을 확일할 수 있다. 추정에 이용된 총 연산시간은 61.5초이며, 이 중 40.3초가 지하수 유동을 예측하는MODFLOW-2005 수행(총 314회)에 이용되었다. 레퍼런스와 Case II 추정 로그 수리전도도의 RMSE는 0.77이었으며, 이는 앞서 Case I의 추정 결과와 유사하다. 두 로그 수리전도도 격자간의 p 값은 0.68로 격자간의 유사성이 상당히 크지만 Case I의 결과에는 미치지 못한다. 추정 결과의 SSIM 값은 0.6의 값을 보여 레퍼런스에 대한 구조적 유사성이 큰 것으로 나타났지만 이 역시 Case I에 비하여 낮은 수치를 보여준다. 즉, Case II에서 추정의 질을 평가하는 모든 지표들은 Case I에 비해 상대적으로 낮은 값을 보였으며, 이는 입력 정보의 감소에 기인한 것으로 볼 수 있다.

Fig. 6. Reference hydraulic conductivity and groundwater level (left hand side) compared to the corresponding hydraulic conductivity and groundwater level estimated by the geostatistical evolution strategy used in this study (right hand side). In the estimation, only the groundwater levels were used whereas the hydraulic conductivity measurements of all monitoring wells were not used (Case II).

Case II 추정을 통하여, 고무적인 사항은 레퍼런스의 평균 로그 수리전도도와 비교했을 때 상당히 낮은 값인 0을 사전정보로 이용했음에도 불구하고 추정된 평균 로그 수리전도도는 1.89로 레퍼런스에 이용된 값 1.89와 동일였다. 그러나 추정된 로그 수리전도도의 분산은 0.68로 산정되어 레퍼런스의 분산 1.5를 저평가 하였다. 이와 같은분산의 감소는 포물선형 방정식을 이용하는 역산해석에서 일반적으로 관찰되는 현상이다. 앞서의 Case I의 경우 관측정에서의 수리전도도 측정값이 이용되었으므로 이러한 분산 저평가 경향이 크지 않았다고 볼 수 있으며, CaseII의 경우 수리전도도 측정값이 이용되지 않았으므로 저평가 경향이 보다 심화된 것으로 해석할 수 있다. 따라서, 측정값이 이용되지 않은 수리전도도 추정결과를 해석할 경우, 실제 수리전도도 분포의 분산이 추정된 것 보다 클 수 있다는 점을 유념하여야 한다.

레퍼런스와 추정된 지하수위 간의 RMSE는 4.21 m로 나타났으며, 이는 Case I에 비하여 약간 작은 값이다. 관측정이 해안에 가까울수록 오차가 크고 상류구배로 갈수록 오차가 작아지는 경향은 Case I과 유사하게 나타났다(Fig. 7). 두 지하수위의 격자간 p 및 SSIM 값은 각각 0.999 및 0.9996로 나타났으며, 이는 레퍼런스와 추정 지하수위의 형태가 거의 일치함을 지시한다. 또한 Fig. 7에서와 같이 레퍼런스와 추정 부정류 지하수위가 모든 관정에서 거의 일치하는 현상을 확인할 수 있다. 이와 같은부정류 지하수위의 일치는 지구통계 진화전략 기법에 의해 보정된 지하수 모델에 의한 미래 지하수위 변화 예측이 신뢰할 수 있음을 의미한다. 지하수 모델이 정류 수위에 기초하여 보정되었을 경우, 이는 과소결정 상태 하에 있다고 볼 수 있다. 따라서, 대부분의 과소결정된 지하수모델은 미래에 대한 예측 능력이 없다고 판단할 수 있으며 최소 동일 지점에서 서로 다른 시기에 2회 이상 측정한 수위를 이용하여 보정을 하여야 미래에 대한 예측능력이 있다고 할 수 있다.

Fig. 7. Groundwater levels monitored in the monitoring wells and the corresponding Case II groundwater level estimations . The monitoring wells in the order of the upper left to the lower right are JDHanlim1, JDHyeopjae, JMMyeongwol, JDPanpo2, JRGeumak1, JRGeumak2, JIWonjongjang, and JDHanlim2, respectively.

Case II 추정을 통하여 추정에 이용되는 정보의 양이 추정의 질에 결정적인 역할을 한다는 일반론적 사실을 관찰하였다. 또한, 지구통계 진화전략 기법의 경우 광역최적화 특성에 의해 사전정보에 크게 좌우되지 않음을 부정확한 평균 로그 수리전도도 입력을 통하여 확인하였다. 앞서 Case I에서 토의하였던 바와 유사하게, 실제 소유역 규모의 지역에 대하여 평균 수리전도도를 파악할 수 있을 만큼 많은 조사가 이루어지기는 매우 어렵다. 따라서, 소유역 규모 대수층에 대하여 역산해석을 통해 수리전도도 분포를 파악하는 경우 부정확한 평균값을 이용하는 것이 보다 일반적이라 할 수 있으며, 사전정보에 민감한 역산해석 모델을 이용하는 경우 추정결과의 신뢰성을 담보하기 매우 어렵다. 지구통계 진화전략은 광역최적화 기법에 근거하므로 사전정보에 대한 추정결과의 민감도가 상대적으로 낮다는 것을 Case II 예제를 통하여 확인하였다. 그러나 Case II에서와 같이 관측정에서의 측정된 수리전도도가 이용되지 않을 경우 포물선형 편미분 방정식을 이용하는 역산해석의 특성상 평활화 문제로 인한 추정 수리전도도의 분산 감소 문제는 심화될 수 있음이 확인되었다.

4. 결론

본 연구에서는 소유역 하부 대수층의 수리전도도 분포를 추정함에 있어 지구통계 진화전략 역산해석 기법의 적용 가능성이 시험되었다. 지구통계 진화전략 기법은 광역최적화 모델인 진화전략과 지역최적화 모델인 앙상블 칼만필터를 연계하여 수리 물성의 공간적 분포를 효율적으로 추정하는 자료동화 기법의 일종이다. 소유역 대수층 특성화 시험을 위하여 제주도 한림읍의 형태, 수리적 경계조건, 그리고 지하수위 관측정 분포가 이용되었다. 또한, 각 관측정에서의 수리전도도 측정값이 있는 경우(Case I)와 없는 경우(Case II)를 상정하여 예측을 수행하였다. 본연구에서는, 역산해석 기법을 통해 추정된 수리전도도 분포와 추정 대상 수리전도도 분포의 비교 용이성을 위하여 가상의 수리전도도 분포를 이용한 가상의 지하수위 시계 열자료를 적용하였다. Case I 추정 결과, 사전정보로 입력된 공간적 상관거리의 부정확성에도 불구하고 대수층 수리전도도의 공간적 분포는 구조적 측면과 격자들 간의 수리전도도 값 측면에서 매우 유사하게 나타났다. Case II 추정 결과, 각 관정에서의 수리전도도 측정값이 없는 것으로 가정하였음에도 불구하고 추정된 수리전도도 분포는 가정된 수리전도도 분포와 상당한 유사성을 보였다. 이와 같이 사전정보에 추정결과가 민감하지 않은 것은 지구통계 진화전략이 광역최적화 기법이기 때문인 것으로 해석되며, 따라서 실제 추정 대상 값들의 사전정보를 얻을 수 없는 실제적인 경우에 실용적으로 이용될 수 있을 것으로 판단된다. 또한, 본 연구는 지하수위 관측정을 통한 지속적인 시계열 관측의 중요성과 정당성을 제시하는 주요 근거가 될 수 있다.

본 연구는 가상의 자료를 이용한 제한성을 지니고 있다. 또한 지하수위를 변동시킬 수 있는 다양한 요인들이 모두 고려되지 않았다는 제약성을 가진다. 이에 따라 현재 실제 지하수위 자료를 이용한 유역 규모 수리특성화 연구가 진행되고 있으며, 그 결과는 관련 학술지에 투고될 예정이다.