1. 서 론

일반적으로 밀도 결합 흐름 문제(density-coupled flow problem) 모의에 사용할 수 있는 수치모의 코드 중 가장 흔히 사용되는 것에는 유한차분법에 기반을 둔 SEAWAT (Guo and Langevin, 2002)과 유한요소법 기반의 FEFLOW (Diersch, 2002)와 SUTRA (Voss and Provost, 2002) 등이 있다. SEAWAT의 소스코드는 USGS의 지하수 모형인 MODFLOW의 수정 버전을 사용하여 염도가 존재하는 지하수 환경에 상응하는 담수 수두(equivalent freshwater head)를 주 종속변수로 모의한다. 그리고 일반적으로 많이 쓰이는 용질이동 코드 MT3DMS를 사용하여 염분 및 기타 용질의 이동 프로세스를 시뮬레이션한다. 유한 차분 방식은 계산이 비교적 간단하며 그 유도 과정이 비교적 직관적이라는 장점이 존재한다.

밀도 결합 지하수 흐름을 모의하는 수치해석 기법을 개발하고 검증하는데 유용한 수치예제들 중 널리 쓰이는 문제들을 Benchmark Problem이라 한다. Henry Problem은 반-해석적 해(semi-analytical solution)가 존재하는 유일한 밀도결합 수치예제로서(Henry, 1964) 오랜 기간동안 연구자들에 의해 검증되고 수정되어 왔다. 이 접근법은 해양 경계선을 향해 바다 쪽으로 이동하는 담수 흐름이 일정하다는 조건 하에서 담염수 경계면의 위치와 모양을 정의한다. Henry의 최초 연구 이후, Henry가 제시한 문제 해결을 위해 Pinder and Cooper (1970), Lee and Cheng (1974), Segol et al. (1975), Frind (1982), Huyakorn et al. (1987), Voss and Souza (1987) 가 몇 가지 다른 해석해와 수치해를 개발했다. 이들 연구는 모델의 매개변수들을 검증하여 그 중 일부를 Henry의 해에 맞춰 수정하는 방식으로 진행되었다. Croucher and O’sullivan (1995)은 이 모든 연구들을 검토한 뒤, Henry가 원래 제시한 해석해에 작은 오류가 있었을 수 있다고 결론 내렸다. Segol (1994)은 Henry의 해와 관련된 오류들을 재검토한 뒤 몇 가지 수정 사항을 제시했다. Simpson and Clement (2004)는 담수의 유량을 수정한 Modified Henry Problem을 선보였다. 또한 Goswami and Clement (2007)는 실내실험을 수행하여 새로운 벤치마크 데이터세트를 개발하였다.

본 연구의 목표는 포화 지하수계의 2차원 밀도 결합 흐름과 관련된 문제를 푸는 효율적인 유한차분 알고리즘을 개발하는 것으로서, Microsoft Excel에 내장된 Visual Basic을 플랫폼으로 활용하였다. 현재 시중에 나와있는 상용 프로그램들이 복잡하고 다양한 문제를 해결하기 위해 정교하고 정확한 코드를 개발하고 있지만, 수치모의를 교육용으로 접할 필요가 있거나 복잡하지 않은 실내실험의 결과를 검증을 하기 위해 간단한 모의 수준에서 수치모델을 구현할 수 있는 모의 기법개발의 필요성도 상시 존재한다. 이에 따라 본 연구는 널리 활용되는 플랫폼을 이용한 단순하면서 강력한 코드의 개발을 목표로 하였다. 이미 Torlapati and Clement (2013)은 1차원 반응성 이송 모델을 Microsoft EXCEL Visual Basic 플랫폼을 기반으로 성공적으로 구현한 바 있다. 본 연구에서 VBA (Visual Basic for application) 기반 2차원 밀도결합 지하수 유동 및 거동 모의 해석 기법을 VDFT (Visual Basic based Density- coupled Flow and Transport)로 명명하였으며 몇 가지 벤치마크 문제들을 이용해서 검증하였다.

이 기법은 압력 기반의 유한차분법을 이용해서 흐름의 지배방정식을 구성한다. 압력 기반 접근법은 다른 수두 기반 방식들에 비해 비교적 간단하고 직접적이다. 용질 이동의 해를 구하기 위해, 이송(advection) 프로세스는 양해(explicit)법인 3rd order TVD (Total Variation Diminishing)법으로 계산하여 담염수경계면 해석의 효율성을 꾀하였고, 분산(dispersion)/유입과 유출(source and sink) 조건은 음해유한차분법(Implicit Finite Difference Method)으로 시도하여 수치적인 안정성을 도모하였다. 본 연구의 차분 과정과 모의결과의 시각적 표기 방식은 Chang (2012)에서 제시한 방법 및 결과를 차용하여 수정하였다. MODFLOW의 연계 프로그램으로, 유한차분법에 기반하여 밀도 기반 지하수 흐름을 모의할 수 있는 검증된 코드인 SEAWAT (Guo and Langevin, 2002)으로 생성한 벤치마크 데이터세트들을 비교하였다.

2. 연구방법

2.1 지배방정식

2차원의 균질한 등방성 다공매질을 통과하는 밀도 결합 포화 흐름의 지배방정식은 다음과 같이 표현할 수 있다. 

(1)

여기서, Sp [M^-1LT²]는 방정식이 압력으로 표시되었을

때의 압력기준 비저류계수(specific storage), P [ML^-1T^-2]는 압력, [ML^-3]는 유체 밀도, 는 함수비, K [LT^-1]는 수리전도도, t [T]는 시간, x와 z는 데카르트 좌표로서 x는 수평방향 y는 수직방향이다. E는 용질의 밀도를 농도로 나눈 값으로서 밀도 경사라고 명명한다.

평형방정식 Eq. (1)의 우변은 미소체적내의 유체 질량 변화를 나타내고 있다. 우변의 첫 번째 항은 다공성물질의 압축성, 유체의 압축성 등에 의해 발생하는 지하수 저류 효과로 인한 유체 질량 축적률을 나타내고, 오른쪽 두 번째 항은 용질 농도 변화로 인한 유체 질량 축적률을 나타낸다.

2차원 용질이동은 이송-분산 방정식에 기반하는데, 지배 방정식의 비선형성을 극복하기 위해 연산자 분할 전략을 도입하면 이송방정식 Eq. (2a)의 해를 푼 이후에 분산방정식은 선형으로 차분화가 가능해진다.

 (2)

                      (2a)

              (2b)

여기서, 와 [L/T]는 지하수의 평균 선형 유속, C [ML^-3]는 용질 농도, 와 [L²/T]는 수리 분산지수를 나타낸다. 우선 연산자 분할을 이용해서 방정식 Eq. (2a)에 TVD 패키지를 적용하여 이송에 의한 농도변화를 update한다. 그리고 이송()에서 구한 농도를 사용해서 음해 유한차분법(Implicit FDM)으로 방정식의 분산 부분을 풀면 방정식 Eq. (2b)가 나온다. 각 격자에서의 평균 선형 유속을 계산하기 위해 단위면적당 유량(specific discharge)을 공극률로 나눈다. 공극률로 나누기 전 값인 다르시 유속(Darcy flux)의 일반식은 다음과 같이 쓸 수 있다:

 (3a)

 (3b)

여기서 , 는 단위 면적당 유량, P는 압력[ML^-1T^-2], 는 유체 밀도[ML^-3], 는 중력 가속도[LT^-2]를 나타낸다. 방정식 Eqs. (1) and (2)은 방정식 Eq. (4)을 통해 연결된다.

 (4)

국제 단위(SI unit)로 표기하였을 때, 일반적인 해수 농도는 약 35kg/m³이고, 담수 밀도는 1,000kg/m³, 해수 밀도는

1,025kg/m³이며, 이에 따라 계산되는 밀도 경사는 0.714이다

(). 실내실험 등의 경우

해수 농도 및 밀도 경사가 일부 달라질 수 있기 때문에 항상 연구 방법에서 매개변수 값을 명시해야 한다.

2.2 2차원 유한 차분 코드 구성

유한차분 모델인 SEAWAT과 VDFT기법의 가장 큰 다른 점은 차분화(Discretization)하는 방식이다. SEAWAT에서는 MODFLOW 프로그램의 수두 표현 방식을 차용하기 위해 밀도 모델의 압력 수두가 담수수두(equivalent freshwater head)로 변환된다. 이 과정을 거치면 흐름 모델은 유한차분 방식 특유의 간결성 및 직관성을 보기 어렵게 된다. 본 연구의 VDFT수치기법은 압력 수두를 주 종속변수로 고수하면서 좀 더 간결한 차분화 방식을 유도하였고 그 결과, 방정식 Eq. (1)의 편미분 방정식으로 표현했던 지하수 흐름 지배방정식의 유한차분 표기 방식은 다음과 같다.

 (5a)

여기서, 위첨자 n은 n번째 이산 시간 수준, tⁿ는 해를 알 때의 시간 수준, 위첨자 m은 Picard 반복 수준(Picard iteration step), 는 격자 ij에서의 압력 수두, 아래첨자 i+1/2, i-1/2는 인접한 두 셀 사이의 값을 나타내며, 두 셀 사이의 밀도는 평균값으로 표시한다. 격자 순서는, 본 연구에서 사용된 모든 식에 대해 VBA 언어에서 매긴 번호 순서를 따라서, 격자 ij를 기준으로 j는 오른쪽이 양의 방향이고 i는 아래쪽이 양의 방향이 된다. 현재 반복(iteration) 수준과 이전의 반복 수준은 각각 m과 m-1로 나타낸다.

(5b)

이후 유한차분식은 다음과 같이 선형방정식으로 정리된다:

(6)

여기서, 

방정식 Eq. (6)은 모든 내부 격자들에 적용하고 경계조건에서는 이 방정식을 수정한다. 그 결과, 다음과 같은 선형방정식 체계가 나온다:

 (7)

여기서, A는 정방 밴드 행렬(square band matrix)로 계수 행렬 (Coefficient Matrix), 는 현재 시간 수준 n과 반복 수준 m에서의 미지수 압력수두 열벡터 또는 해벡터s(Solution vector), b는 방정식 Eq. (6)의 RSH항을 표현한 열벡터이다, 이 선형대수방정식들의 집합은 가우스소거 도구(Gauss–elimination solver)의 밴드 버전으로 풀 수 있다.

방정식 Eq. (2)의 편미분 방정식으로 표현했던 용질 이송 및 분산 방정식의 유한차분 표기 방식은 다음과 같다. 이송에 TVD 기법을 적용하는데 사용한 선형 방정식에서 방정식 Eq. (7a)는 이송 방정식의 공간항, 방정식 Eq. (7b)는 시간항을 표현한다.

 (7a)

 (7b)

TVD 기법에서는 전체적인 플럭스 제한 절차(universal flux limiting procedure)를 채택하여, 농도가 급격히 변하는 지점(sharp concentration front)과 관련하여 일어날 수 있는 수치상의 진동(numerical oscillation)을 최소화한다. 해수침투의 경우 담염수 경계면이 농도가 급격히 변하는 지점에 해당한다. 이 기법은 수치가 과도하게 분산되지 않게 하면서 질량을 보존하며, 기본적으로 진동에서 자유롭다. 분산은 이송에 비해 시간간격 에 크게 구애받지 않기 때문에 안정적으로 음해 유한차분법을 아래와 같이 도입할 수 있다. 방정식 Eq. (8a)은 공간항을 기준으로, 그리고 방정식 Eq. (8b)은 시간항을 기준으로 표현한 것이다.

 (8a)

 (8b)

공간 항에서 사용된 분산 텐서(tensor)의 평균값이 반복의 이전 단계인 m-1에서 업데이트된 값이란 점에 주목했을 때, 이 값을 명시적으로 계산해서, 농도가 방정식에서 선형 행렬 계산에 사용할 유일한 독립변수가 되게 한다. 유한차분방정식 Eqs. (8a) and (8b)의 일반형은 다음과 같다:

(9)

여기서,

방정식 Eq. (9)는 모든 내부 격자들에 적용된다. 경계조건에 해당하는 격자에서는 적절한 경계조건이 반영되도록 이 방정식을 수정한다. 그 결과로 도출되는 선형 방정식 체계는 밴드식 가우스소거 도구로 푼다.

경계 조건의 경우, 다음 절에서 소개하는 예제 1과 2 모델의 담수경계에는 수두의 값이 지정되는 Dirichlet 경계조건을 사용하며 예제3은 담수 유량을 설정하는 Neumann 경계 조건을 사용한다. 또한, 모델 바닥은 모든 검증문제에서 흐름이 없는 경계로 설정한다. 해수침투는 경계의 농도를 바꾼다. 자세한 변경내용은 각 검증문제별 모형 설정에서 표시된다.

시간 간격을 설정하는데 있어 시간 간격 를 제한하려면 쿠란 조건(Courant condition)이 필요하다. 시간 간격을 충분히 작게 유지하기 위해 모든 모의에 쿠란 조건을 적용한다. 쿠란 조건으로 시간 단계 길이를 계산하여, 용질이동 방정식 중 이송 방정식의 안정성 제약과 정확성 요건을 충족시킨다. 본 연구의 모의에서는 쿠란 값(Courant number)을 0.5로 설정하였다.

3. 연구결과

3.1 수치예제 1: Sinking plume model 모의 계산 결과

첫 번째 검증문제는 고농도 용질이동 사례이다. Goswami et al. (2012)는 SEAWAT/MT3DMS에서 제공하는 수치기법 3개를 이용해서 2차원 흐름 탱크의 실험해와 수치해를 제공했다. 이 실험들은 하강하는 지하수 플룸과 상승하는 지하수 플룸을 나타내도록 설계되었다. 이 문제는 오염물질 이동, 해안 대수층의 쓰나미 문제 같은 현실 환경문제에도 도움이 될 수 있다. 영역조건(domain condition)과 경계조건은 피압대수층을 묘사한다. 크기가 고정된 담수 압력 수두가 주 영역의 왼쪽 경계와 오른쪽 경계에 위치한다(Fig. 1 참조). 2차원 모델링 영역은 길이 22 5mm, 높이 180 mm이다. 종분산 값은 0.1mm로, 횡분산 값은 0.01mm로 설정한다. 3개의 스트레스 기간을 채택하여, 다음 3개의 실험단계를 모의했다: (1) 주입 전 정류상태 조건, (2) 담수/염수 슬러그 주입, (3) 슬러그가 탱크를 통과하며 이동. 시간 단계 크기는 처음에는 최소값 0.1초 (second)와 최대값 1초로 정했다. Courant 값은 0.5로 정했다. 이 모델에서 사용한 값들은 Table 1에 나와 있다.

세 번째 상태(이동 단계)의 모의 결과를 SEAWAT 결과와 비교한 내용은 Fig. 2를 통해 확인할 수 있으며 SEAWAT 결과와 비교하였을 때 비교적 유사함을 알 수 있다. 밀도 의존 하강 연구는 전형적인 해수 또는 염분 농도와는 구성이 다른 유체에도 적용될 수 있다. 예를 들어 Kanel et al. (2008) 은 영가철 나노입자(INP)의 확산과 이동을 시뮬레이션한 후, 포화된 정류상태 흐름조건에서 폴리아크릴산(poly acrylic acid)(S-INP)을 이용해서 열가철 나노입자를 다공성매체 속에 안정화시켰다.

3.2 수치예제 2: 해수침투 실내실험 모델 모의 계산 결과

이전까지 수치모델 결과만으로 벤치마킹 하던 방식이었던데 반해 Goswami and Clement (2007)는 수치예제 문제에 사용할 수 있는 실내실험 데이터세트를 제안했다는 점에서 예제 2는 예제 1및 예제 3과의 차이를 보인다. Fig. 3은 Goswami and Clement (2007)의 실내실험으로 사용한 수조 및 실험 설정들을 개념적으로 구현한 도식이다. 모델 모델 영역은 피압대수층을 묘사한며, 값이 고정된 압력 수두가 주 영역의 왼쪽 경계와 오른쪽 경계에 위치한다(Fig. 3참조). 1.0cm 크기의 격자를 전체적으로 적용해서 53.0cm × 26.0cm의 차원모델을 구획했다(x축으로 54칸, z축으로 26칸). 염수 수두는 왼쪽 경계 25.5cm 높이로 유지하고, 담수는 오른쪽 정수두 경계(constant-head boundary)에서 모델 영역으로 유입된다. 비저류계수 값(S_s)은 모든 격자에서 10^-6cm^-1로 정했다. 모형 바닥에는 무흐름경계로 설정하였다. 종분산지수와 횡분산지수는 각각 0.05cm와 0.005cm를 사용했다. 이 모델에서 사용한 모든 매개변수 값들은 Table 2에 나와 있다.

3개 정류상태 조건을 선택해서, 본 연구에서 개발한 수기기법으로 모의한 데이터를 SEAWAT 모의 결과와 비교했다. Fig. 4의 VDFT SS1, VDFT SS2, VDFT SS3는 각각 담수 경계가 26.7, 26.2, 26.55cm으로 설정되었을 때의 수치모의로 얻은 정상상태의 염수-담수 경계면 위치이다. 본 연구에서의 경계면 위치는 염도(salinity)가 50%인 위치로 정의하며, 각 격자의 농도를 이용해 내삽법(interpolation)으로 계산한다. 정상상태의 모의 결과인 세 개의 경계면의 위치는 SEAWAT 모의 결과와 잘 부합하며, 해수침투조건과 후퇴조건에 대한 부정류상태 모의값 역시 SEAWAT 결과와 잘 부합한다(Fig. 5 참조).

3.3 수치예제 3: Flux-type Henry Problem 모의 계산 결과

초기의 Henry Problem (Henry, 1964)와 Modified Henry Problem (Simpson and Clement, 2004)에 기반한 개념모형들을 모의했다. Simpson and Clement (2004)은 물리적으로 실내실험 규모의 매개변수값을 수치 예제인 Henry Problem에 부여하였다. 밀도에 의한 흐름 및 이동 모의 기법은 밀도항에 의해 결합(coupling)되어 반복 연산이 실행되는데, 유체 내 밀도의 영향을 좀 더 명확히 알 수 있도록 Henry Problem의 모든 매개변수값을 유지한채 담수 유입 유량을 절반으로 줄인 Modified Henry Problem을 제시하였다. 밀도 항에 의한 두 방정식의 결합 및 반복 연산 작업이 정확하지 않아도 담수 유입량이 큰 경우는 어느정도 모의 결과가 검증된 수치예제와 유사하게 보일 수 있다고 판단했기 때문이다. 본 연구의 세번째 예제는 Simpson and Clement (2004)이 제시한 모델 설정을 따른다. Fig. 6에서 보이는 바와 같이 200cm × 100cm크기의 모델을 5.0cm 크기의 격자로 구획하여 2 차원 모형을 구성했다. 우측 경계에 해수가 위치해있으며 수두값은 100cm이고 왼쪽 경계에는 유량이 고정된 담수유입이 있다. 비저류계수 값(S_s)은 모든 격자에서 10^-6cm^-1로 설정했다. 모든층을 피압대수층으로 설정하였고 모델 바닥면은 무흐름경계로 설정하였다. 종분산지수와 횡분산지수는 모두 0으로 가정한다.

Henry Problem (a)와 Modified Henry Problem (b)에 본 연구의 수치기법을 적용한 결과와 SEAWAT 결과를 비교하여 Fig. 7에 나타내었다. 예제 1과 예제 2와는 달리 VDFT가 SEAWAT 의 모의 결과와 일부 영역에서 차이가 나는 것을 볼 수 있다. 두 기법 간 차이는 Modified Henry Problem에서 명확히 드러난다. Simpson and Clement (2004)가 담수의 유량을 반으로 줄이며서 밀도 영향에 대한 코드의 민감성을 더욱 크게 보여줄 수 있다고 주장한대로, VDFT 코드와 SEAWAT 모의 결과의 차이가 Modified Henry Problem에서 상대적으로 크게 나타난다 이는 Henry Problem에서 정의한 무차원 변수 a 값의 비교로도 알 수 있는데

Henry Problem의 무차원 변수값 a는 0.263, Modified Henry

Problem의 무차원 변수값 a= 0.1315 이다. 의

정의에 따라 계산한 예제 2의 무차원 변수값 a은 0.910이 나와 상대적으로 예제 2의 담수 유입량이 큰 것으로 볼 수 있다. 이 계산과정에서 Henry Problem의 무차원변수에 적용한 매개변수값은 Table 3에 나와 있다.

4. 결 론

본 연구는 2차원 밀도 결합 유한차분 코드인 VDFT를 EXCEL 기반 Visual Basic 플랫폼을 이용하여 개발한 것이다. 교육용 및 실내 실험 검증을 위해 사용하고자 비교적 간단하고 강력한 차분 방식을 도입하였으며 그 결과 실내 실험 결과를 바탕으로 한 예제 문제를 벤치마크 문제로 비교 검증하였을 때는 매우 유사한 결과를 성공적으로 도출해내었다. VDFT 코드가 Henry Problem 타입의 모델에서는 실제 실내 실험보다 담수의 영향력을 약화시킨 경계 조건의 구성으로 인해서 기존 코드와 일부 차이가 발생하였기 때문에, 본 연구에서 개발한 모의 기법은 유체 내 밀도가 지하수 흐름에 미치는 영향을 정교하게 예측하고 정확성을 높이는 결과를 얻기 위해 향후 연구가 지속되어야 한다. 현재 개발 수준에서도 실험 결과의 비교를 위한 활용 가능성을 보여주고 있는데 본 연구에서 개발한 수치기법은 해수침투 연구 뿐만 아니라 밀도 차이를 두는 유체 간의 상호작용 문제인 Hele-Shaw 타입의 문제에서도 사용 가능할 것으로 보인다. 향후 본 연구결과는 VDFT 코드 개발을 위한 더 많은 예제 발굴에 활용될 수 있을 것으로 판단된다.