2.1 활성슬러지 종속영양미생물 성장-사멸 동역학 배치 실험

하수에 포함된 Readily Biodegradable COD (RBCOD)와 하수처리 시스템 활성슬러지의 호기성 미생물 분율을 정량 분석을 위하여 미생물 성장-사멸 동역학 실험(Microbial Growth Kinetic Test)을 수행하였다(^{Choi et al., 2017}; ^{Lee et al., 2006}; ^{Wentzel et al., 1995}). Modified Ludzack-Ettinger (MLE) 타입 하수처리 시스템에서 추출한 활성슬러지와 하수를 혼합하여 3 L 회분식(Batch) 반응기에 투여 후 실험을 진행하였다(Fig. 1a). 반응기 용존산소 측정을 위한 용존산소 프로브(YSI Inc., 미국)와 산소 공급을 제어할 수 있는 PLC 시스템을 활용하여, 회분식 실험이 진행되는 동안 산소호흡율(Oxygen Uptake Rate, OUR)를 지속적으로 측정하였다. 용존산소 농도가 2 mgO/L 이하로 하락하면 산소공급이 자동으로 시작되며, 4 mgO/L 이상에 도달하면 산소공급이 차단되고, 산소공급이 중단된 기간 동안 측정된 용존산소 농도 저감율을 바탕으로 산소호흡율을 산정하였다. 또한, 회분식 실험이 진행되는 동안, 일정 시간 간격으로 여과 샘플을 확보하여 질산염 및 아질산염 분석을 수행하였고 질산화에 따른 산소호흡율을 산정하였다. 본 연구는 유기물 분해 종속영양미생물(Ordinary Heterotrophic Organisms)에 한정하여 미생물 성장-사멸 동역학 모델을 수립하였으므로, 질산화에 따른 산소호흡율(OUR_N)을 총 산소호흡율(OUR_T)에서 감하여 계산한 종속영양미생물에 의한 산소호흡율(OUR_H)을 모델 매개변수 보정 작업의 실험데이터로 활용하였다(Fig. 1). 회분식 실험 시작과 종료 시 샘플에 대한 COD 농도를 측정하여 COD 물질수지 수립에 활용하였다. 미생물 성장 동역학 실험의 산소호흡율 자료에 대하여 해석적 계산법(Analytical Solution)을 적용하면, 하수의 RBCOD 분율과 활성슬러지에 포함된 호기성 미생물 분율을 산출할 수 있다(^{Lee et al., 2014}; ^{Wentzel et al., 1995}). 하지만, 본 연구는 해석적 계산법을 적용하지 않고 연립 미분방정식 형태의 수학적 모델을 베이지안 보정 기법과 연계하여 하수 RBCOD 분율과 활성슬러지 호기성 미생물 분율을 산정하였다.

Fig. 1. (a) Schematic diagram of a batch test for microbial growth kinetics, and (b) Experimental results obtained from the batch test (i.e., changes of nitrate and nitrite concentration and oxygen uptake rate). OURTindicates total oxygen uptake rate. OURNand OURHindicate oxygen uptake rates by nitrification and heterotrophic microorganisms, respectively.

2.2 응집 동역학 실험

응집 동역학의 최적 모델 계수 산정 도구로써 베이지안 보정 기법의 활용성을 평가하기 위하여 응집 동역학 테스트를 수행하였다. 생체고분자물질(Extracellular Polymeric Substances)의 일종인 Xanthan Gum (Sigma-Aldrich, USA)과 점토질 입자를 포함한 현탁액을 사용하여 2 L 표준 쟈테스터(Phipps and Bird, USA)에서 응집 동역학 테스트를 수행하였다(^{Kim and Lee, 2020}). Xanthan Gum은 다당류 중합체로써 입자 응집을 증진시키는 것으로 알려져있다(^{Villacorte et al., 2015}). 또한, 0.1~4.0 μm 크기 분포의 Kaolinite (Sigma-Aldrich, USA)를 점토질 입자로써 실험에 사용하였다. 응집 동역학 실험 시작 후 2분 동안, 반응조 현탁액을 200 rpm으로 급속 혼합하였으며, 이후 응집 동역학이 평형 상태에 도달할 때까지 현탁액의 난류 강도를 제어하기 위해 혼합속도를 100 rpm으로 고정하여 교반하였다. 응집 동역학 실험을 진행하는 동안, 반응조와 광산란식 입도분석기(LISTT-200X, Sequoia Inc., USA)를 6mm 튜브와 정량펌프로 연결하여 현탁액의 입도 분포를 실시간 관측하였다(Fig. 2a). LISST-200X 입도 분석기는 Mie 산란 이론에 기반하며 레이저가 수중 입자에 부딪힐 때 산란되는 빛을 분석하여 입도분포를 실시간 분석하는 장비로써, 본 연구의 연속류 흐름(Flow-through) 장치와 연계하여 매 30초 간격으로 반응조 내 입도분포를 관측하였다. 현탁액의 입자 농도가 높을 경우, 광산란식 입도분석기는 레이저 경로상 중첩에 의한 오류가 발생하므로, 본 연구의 점토질 현탁액의 농도는 0.12 g/L으로 제한하였다. 관측된 입도 분포는 응집 동역학 모델의 주요 종속변수인 응집체 평균 직경(D_F)과 부피 분율(V_F)을 환산되어, 응집 동역학 모델과 베이지안 보정 기법 적용성 평가에 활용되었다(Fig. 2b).

Fig. 2. (a) Schematic diagram of the batch test for flocculation kinetics and real-time measurement of floc size distribution, (b) Experimental results obtained from the batch test (i.e., mean diameter and volumetric fraction of flocs).

2.3 활성슬러지 종속영양미생물 성장-사멸 동역학 모델

본 연구의 생물학적 공정(반응) 동역학에 대한 베이지안 보정 기법 적용성 평가를 위하여, ASM1 (Activated Sludge Model No.1) 모델에 기반한 호기성 종속영양미생물 성장 동역학 모델을 활용하였다(^{Henze et al., 1987}; ^{Jeppsson, 1996}). ASM1 호기성 미생물 성장 동역학 모델은 Monod 타입 모델을 기초로 최대 성장속도 상수와 반포화 상수를 포함한 성장 동역학과 1차 반응식 형태의 미생물 사멸(자가분해) 동역학으로 구성되며(Equation 1), 미생물 기질(먹이)은 빠르게 분해되는 Readily Biodegradable COD (RBCOD)와 느리게 분해되는 Slowly Biodegradable COD (SBCOD)로 구성된다(Equation 2~3). 특히, SBCOD는 미생물 가수분해(Hydrolysis) 작용에 의해 일정 분량은 RBCOD로 변환되며 나머지는 난분해성물질로 변환되는 것으로 가정된다. 미생물 성장 동역학 모델은 각 호기성 미생물, RBCOD, SBCOD에 대한 연립 미분방정식 형태로 구성되며, 산소소모율(Oxygen Uptake Rate; OUR)은 시스템의 상태 변수로써 산출될 수 있다(Equation 4). 호기성 미생물 성장 동역학의 연립미분방정식은 Matlab (The MathWorks Inc., USA) 내장함수를 활용하여 적분되었으며, 베이지안 보정 기법과 연계하여 다수 모델 매개변수 조합에 대하여 연립미분방정식을 반복적으로 연산하였다. 호기성 미생물 성장 동역학의 연립미분방정식 계산을 위한 Matlab 코드는 Appendix 1~3에 수록하였다.

여기서, X_H = 호기성 미생물 농도 (mgCOD/L), S_S = Readily Biodegradable COD (RBCOD) 농도 (mgCOD/L), X_S = Slowly Biodegradable COD (SBCOD) 농도, Y_H = 기질의 생체 전환율 (-), μ_H = 최대 성장속도 상수 (/d), K_S = RBCOD 분해 반포화 상수 (mgCOD/L), k_H = 가수분해 속도 상수 (/d), K_X = SBCOD 분해 반포화 상수 (mgCOD/L), f_p = SBCOD 가수분해에 따른 난분해성물질 전환율 (-), b_H = 호기성 미생물 자가분해 상수 (/d), O = 용존산소 농도 (mgO/L), 그리고 OUR = 산소소모율 (mgO/L/d). 본 연구의 베이지안 보정 기법 적용 시, 배치 실험에서 획득된 종속영양미생물의 산소호흡율(OUR_H) 데이터를 실험-모델 접합을 위한 상태 변수로 활용하였으며, 미생물 성장 동역학 모델 매개변수 4개(μ_H, K_S, k_H, K_X) 그리고 각 종속변수 초기값 3개(X_H,0, S_S,0, X_S,0)를 데이터-모델 접합(Fitting) 계수로 활용하였다. 미생물 사멸 동역학 계수 b_H 및 화학양론 계수 Y_H, f_p는 불확실성이 크지 않으므로 ASM1 모델에서 제시된 대푯값(0.62, 0.67, 0.08)을 사용하였다(^{Henze et al., 1987}; ^{Lee et al., 2014}) (Table 1; Appendix 2).

2.4 응집 동역학 모델

본 연구의 응집 공정(반응) 동역학에 대한 베이지안 보정 기법 적용성 평가를 위하여, 이군집(Two-Class) 응집 동역학 모델을 적용하였다(^{Lee et al., 2011}). 물환경에서 생체고분자물질(EPS) 매개되는 응집은 응집핵(Flocculi)과 응집체(Floc) 이군집 입도분포를 형성하는 것으로 알려져 있다(^{Kim and Lee, 2020}; ^{Lee et al., 2012)}. 따라서, 본 연구는 입경이 고정된(~5 μm)을 가진 응집핵(Flocculi)과 변화(성장)하는 입경을 가진(10~500 μm) 응집체(Floc)로 구성된 이군집입자개체군수지식(Two Class Population Balance Equation; TCPBE)을 적용하였다. TCPBE는 응집핵 농도(dN_P/dt), 응집체 농도(dN_F/dt) 및 플록에 결합된 응집핵의 총 농도(dN_T/dt)의 시간 변화율을 예측하는 연립 미분방정식 형태로 구성된다(Equation 5~7). TCPBE는 응집체에 결합된 응집핵의 수(NC)를 응집체 크기 지표로 사용하였으며, 충돌 빈도 인자(β)와 응집체 파괴 동역학 상수(a_F)에 대한 커널은 수체 전단력(G) 및 응집체 크기와 상관한 공식을 사용하였다(Equation 8~9)(^{Burd and Jackson, 2002}; ^{Lee et al., 2011}). TCPBE 연립미분방정식은 Matlab (The MathWorks Inc., USA) 내장함수를 활용하여 적분되었으며, 베이지안 보정 기법과 연계하여 다수 모델 매개변수 조합에 대하여 연립미분방정식을 반복적으로 연산하였다. TCPBE 연립미분방정식 계산을 위한 Matlab 코드는 Appendix 4~6에 수록하였다.

여기서, D_P와 D_F는 각각 응집핵, 응집체의 평균 직경을 나타낸다. f는 응집체 파괴에 의해 생성된 응집핵의 분율이다. α_PP, α_PF, α_FF는 각각 응집핵-응집핵, 응집핵-응집체, 응집체-응집체 충돌 시 최종 부착되는 충돌 효율을 나타내며, 일반적으로 모델 보정 매개변수로 사용된다. 충돌 빈도 인자(β)는 수체 난류에 따른 입자 충돌의 확률을 나타내며, 아래첨자 i와 j는 응집체(P) 또는 응집체(F)를 나타낸다. E_b는 응집체 파괴 동역학의 효율, μ는 유체의 절대 점도, G는 전단 속도, F_y는 응집체 전단 강도, p와 q는 실험적으로 결정된 경험적 매개변수이다. 본 연구의 베이지안 보정 기법 적용 시, 응집 동역학 실험에서 획득된 응집체 평균입경(D_F) 및 부피분율(V_F)을 실험-모델 접합을 위한 상태 변수로 활용하였으며, 충돌 효율 인자(α_PP, α_PF, α_FF), 응집체 프랙탈 상수 (n_f), 응집체 파괴 동역학의 경험적 매개변수(E_b, p, q), 응집체 파괴에 의해 생성된 응집핵 분율(f)을 데이터-모델 접합 매개변수로 활용하였다. 이와 더불어, 응집핵의 초기 부피분율(V_P,0)을 실험-모델 접합 매개변수로 활용하여 TCPBE 연립미분방정식의 초기조건(N_P,0,, N_F,0,, N_T,0)을 산출하였다(Table 2; Appendix 5). 불확실성이 비교적 작은 응집 동역학 매개변수(Fy)는 기존 문헌의 대푯값(1.0×10^-10)을 적용하였으며(^{Winterwerp and van Kesteren, 2004}; ^{Lee et al., 2011}), 반응조 유체 흐름 및 입자 특성 매개변수들은 표준 및 실측값(μ= 1.0×10^-6, G = 95, D_P = 5.0×10^-6, ρ_P = 2500, D_F,0 = 15.0×10^-6, C_mass = 0.12)을 활용하였다(Table 2).

2.5 베이지안 보정 기법

베이지안 보정 기법은 베이지안 추론(Bayesian Inference)을 기반하는 확률론적 데이터-모델 보정 방법으로써, 관측 데이터-모델 결과 간 최적 보정(Fitting) 매개변수를 확률론적으로 추론하는 방법이다. 베이지안 보정 기법은 알려진 데이터를 기반으로 미지의 사후 분포를 추론하는 방법론이다(Equation 10). 베이지안 보정 기법은 비선형 또는 선형 관계를 설명할 수 있으며, 모델의 불확실성 정량 평가할 수 있는 장점을 가진다(^{Beckers et al. 2020}).

본 연구에서는 베이지안 추론을 기반으로 하는 DREAM (Differential Evolution Adaptive Metropolis) 소프트웨어 패키지(https://faculty.sites.uci.edu/jasper/software/)를 사용하여 미생물 성장-사멸 동역학 및 응집 동역학 모델 매개변수를 추정하였다. 베이지안 보정의 결과는 실험 및 관측 데이터에 대하여 모델 접합(Fitting) 매개변수의 불확실성을 반영하여 매개변수의 사후 확률 분포를 나타낼 수 있다. 매개변수의 사후 분포는 사전 분포와 실험/관측 데이터(e.g., 미생물 성장-사멸 동역학 실험의 산소호흡률, 응집 동역학 실험의 응집체 크기 및 분율)를 고려한 매개변수의 확률(빈도) 함수로부터 계산된다. 본 연구는 베이지안 보정 기법의 적용에 초점을 두었으며, 자세한 수학적 알고리즘은 DREAM 개발자가 출간한 연구논문에 수록되어 있다(^{Vrugt, 2016}).

3.1 최적 접합 매개변수 추적 알고리즘

본 연구의 베이지안 보정 기법은 DREAM 소프트웨어(^{Vrugt, 2016})를 사용하였으며, 특별히 DREAM_(ZS) 알고리즘을 사용하여 최적 접합 매개변수 추적하였다. Fig. 3(a)와 같이, DREAM_(ZS) 알고리즘은 병렬연산(Parallel Computation)을 활용하여 1회 계산 배치에 독립적인 3개 매개변수 조합에 대하여 시뮬레이션을 동시에 수행하여 연산 시간을 단축하였다. 각 회차 계산 배치의 연산 결과를 바탕으로 최적 접합에 더 근접한 매개변수를 도출하여 후속 연산에 활용하므로, 시뮬레이션 결과가 점진적으로 최적 접합에 도달하게 된다. 본 연구의 미생물 성장-사멸 동역학 X_H,0 매개변수 경우, 대략 1,000회 계산 배치 수행 후 최적값인 49.2 mg/L에 도달하는 것으로 나타났다(Fig. 3a). 또한, 계산 배치가 진행되면서 더 이상 발산하지 않고 안정화되는 여부(Convergence)를 판단하기 위하여 R-diaganostic 인자를 바탕으로 결정된다(^{Gelman and Rubin, 1992}). 본 연구 미생물 성장-사멸 동역학 테스트의 경우, 12,000회 계산 배치 후 안정화에 도달하는 것으로 나타났다(Fig. 3b).

Fig. 3. (a) Trace plot of sampled values of a fitting parameter (XH,0of the microbial growth-decay kinetics), and (b) Evolution of theRconvergence diagnostic ofGelman and Rubin (1992).

3.2 미생물 성장-사멸 동역학

미생물 성장-사멸 동역학 배치 실험을 수행하여 획득한 산소호흡율(OUR) 곡선은 Fig. 4 같이 나타나며, 실험초기 호기성 미생물의 RBCOD 분해에 의해 OUR의 급격한 증가가 나타나며 RBCOD 분해가 종료되는 시점에 OUR이 급격히 하락하며, 이후 SBCOD 및 미생물 자가분해(사멸)에 따라 OUR이 일정하게 유지되는 경향을 가진다. 미생물 성장-사멸 동역학 실험의 이와 같은 OUR 변화 경향을 활용하여 하수에 포함된 RBCOD 농도(S_S,0)를 계산할 수 있으며, 또한 활성슬러지의 호기성 미생물 농도(X_H,0) 계산에 활용할 수 있다(^{Lee et al., 2014}; ^{Wentzel et al., 1995}).

Fig. 4. Time series of oxygen uptake rate measured (dotted), best-fit simulation (red line), and 95% simulation uncertainty intervals due to total uncertainty (gray band).

기존의 RBCOD 및 호기성 미생물 농도 계산은 연립미분방정식을 간략화, 실험 데이터에 선형화된 모델을 최적 접합(Fitting)하는 해석적 계산법(Analytical Methods)을 활용하였으나, 본 연구는 베이지안 보정 기법을 활용하여, 비선형 연립미분방정식 변경 없이, 사용자가 정의한 모델 매개변수의 사전 분포(Prior Distribution)에 대한 결과, 즉 사후 분포(Posterior Distribution)를 확률론적인 기법으로 산출하게 된다. Table 3과 같이, 사용자가 정의한 일정 범위를 가진 모델 매개변수의 사전 분포(Prior Distribution) 값(수천~수만 개)들에 대하여 연립미분방정식을 연산하고 데이터-모델 접합도를 평가함으로써 사후 분포(Prior Distribution) 값들을 얻게 된다. 본 연구에서는 번거로운 모델의 선형화, 단순화 과정 없이, 호기성 미생물 성장 동역학 모델의 주요 매개변수인 하수 RBCOD 농도(S_S,0), 활성슬러지 미생물 농도(X_H,0)를 평가할 수 있었다. 다만 베이지안 보정 기법은 수천~수만 개 연립미분방정식 풀이 연산을 해야 하지만, 현재의 가정용 PC 수준으로도 10분 내로 결과를 도출할 수 있다.

모델 매개변수의 사후 분포는 Fig. 5와 같이 일반적으로 정규분포의 형태를 가지지만, 일부 민감도가 낮은 매개변수는 명확한 정규분포를 띄지 않고 평탄한 형태를 가진다. 이와 같이 비교적 평탄한 사후 분포를 가지는 매개변수들은 모델 연산 결과에 큰 영향을 미치지 아니하므로 데이터-모델 접합 작업에서 제외할 수 있다. 또한, 모델 매개변수들 중 상관성이 매우 높은 경우가 발생하며, 예를 들어, μ_H~X_H,0 또는 k_h~S_S,0의 경우, 사용자의 판단에 따라 상응하는 두 개 매개변수 중 하나를 취사선택하여 데이터-모델 보정 작업의 복잡성을 줄일 수 있다. 이와 같이 베이지안 보정 기법은 최적 매개변수의 도출 기능뿐만 아니라, 모델 예측 결과에 대한 매개변수의 영향과 매개변수들 간 상관성에 대한 정보를 추가 제공할 수 있다. 최근 대세인 머신러닝, 인공지능 기법들과 비교하여서도, 베이지안 보정 기법은 매개변수 특성 평가에 관한 특별한 장점을 가지고 있다.

Fig. 5. Histogram of posterior distribution of microbial growth-decay kinetic parameters used in Bayesian calibration and correlation plots between different parameters

3.3 응집 동역학

본 연구에서는 응집핵-응집체 이군집 응집 동역학 모델(TCPBE)을 적용하였으므로, 응집체 크기(D_F)와 전체 응집핵-응집체 농도 대비 응집체 분율(V_F)을 실험 데이터-모델 보정 작업에 활용하였다. 응집 동역학 실험 동안, 응집을 통하여 응집핵의 상대적 분율을 감소되고, 응집체의 크기와 분율은 지속하여 커지며 더 이상의 변화가 없는 평형 상태에 도달한다(Fig. 6). TCPBE는 이와 같은 응집핵-응집체의 변화 경향을 모사할 수 있으며, 베이지안 보정 기법을 적용하여 실험 데이터-모델 예측치 간 최적 접합의 결과를 도출할 수 있었다(Table 4).

Fig. 6. Time series of (a) floc diameter and (b) volumetric fraction measured (dotted), best-fit simulation (red line), and 95% simulation uncertainty intervals due to total uncertainty (gray band).

앞서 소개한 미생물 성장-사멸 동역학의 경우와 유사하게, 사용자가 정의한 응집 동역학 모델 매개변수의 사전 분포(Prior Distribution)에 대한 결과(Posterior Distribution)를 확률론적인 기법으로 산출하게 된다. 데이터-모델 접합을 위한 다수 매개변수들에 대한 일괄적인 평가는 어려운 작업이며, 인위적인 시행착오(Trials and Erros)로 산출하기는 거의 불가능하다. 베이지안 보정 기법은 수천, 수만회의 반복 작업을 컴퓨터를 활용하여 수행하며, 그 결과를 통계적으로 해석함으로써 비교적 손쉽게 데이터-모델 보정 작업을 수행하게 된다.

응집 동역학 모델 매개변수의 사후 분포(Posterior Distribution)는 대체로 정규분포를 따르며, 특별히 응집체 충돌 효율 계수(α_FF)는 사용자가 정의한 사전 분포의 상단 경계에 최빈값이 형성되는 것으로 나타났다(Fig. 7). α_FF는 물리적으로 1을 넘을 수 없으므로, 매개변수의 사전분포 범위 설정에서 현실의 물리적인 범위 내로 한정하였기 때문에, 사후분포는 상단 경계 1을 넘지 않고 그 주변으로 집중되어 있었다. α_FF와 응집체 프랙탈 차원(형상) 계수(n_f)는 강한 상관성을 나타내며, 이들 매개변수가 공통으로 응집 동역학을 증대시키는 것으로 파악할 수 있었다. 또한, 응집체 파괴와 관련된 세 개의 매개변수들(E_b, p, q) 간 뚜렷한 음 또는 양의 상관성을 나타내고 있다. 이와 같이, 모델 매개변수들의 상관성을 파악함으로써, 모델의 구조적 특성 이해를 도모하고 사용자의 모델 적용의 방향에 따라 매개변수를 취사 선택할 여지를 제공할 수 있다.

Fig. 7. Histogram of posterior distribution of flocculation kinetic parameters used in Bayesian calibration and correlation plots between different parameters.