Home

맥스웰 방정식은 시공간적으로 균질한 매질 안에서 전자기장이 어떻게 전파되어 나가는지를 알려주는 이론이다. 산란체 내부와 외부(예, 구름입자와 공기) 또는 산란체 내부에서의 불균질(예, 빙정 안에 갇혀 있는 공기방울)이 발생하는 경계면에서는 구성관계식의 변화와 불연속이 나타난다. 이 경우 경계면에서 서로 다른 매질의 장(예, 전기장)을 연결하는 추가적인 조건이 필요하다. 이 추가적인 경계 조건과 해를 구하는 방법에 따라서 단일산란 특성 직접 계산법의 분류가 이루어진다.

2.2.1.1 분석적 직접 계산법

직접 계산법은 크게 분석적(analytical)인 방법과 수치적(numerical)인 방법으로 나누며 수치적인 방법에는 많은 세부 방법이 존재한다. 맥스웰 방정식에 대한 분석적 직접 계산법은 변수 분리법(separation of variables method)을 사용하여 시간 조화 전기장(time-harmonic electric field)에 대해 벡터 Helmholtz 방정식을 푸는 문제로 귀결된다. 가장 대표적인 방법이 로렌츠-미(Lorenz-Mie) 이론이다. Figure 6은 로렌츠-미 코드를 사용하여 계산한 다양한 크기모수(size parameter, χ=2πr/λ)를 갖는 액체상 구형 입자의 단일산란 특성을 보인다. 여기서 r은 산란체의 반지름이다. 로렌츠-미 이론으로 대표되는 변수 분리법을 사용하는 분석적 직접 계산법으로 구한 맥스웰 방정식의 해는 구면(spherical), 원기둥(circular cylindrical), 회전 타원(spheroidal) 좌표계에 대해서만 존재한다(Liou and Yang, 2016). 따라서 변수 분리법은 구, 다층 구(multilayered sphere), 방사적으로 불균질한 구(radially inhomogeneous sphere), 원기둥, 회전 타원체 등의 비교적 단순한 모양의 산란체에 대해서만 적용 가능하다. 이 분석적 방법을 사용해서 구한 해는 단순한 모양의 산란체에 대해서도 이미 충분히 복잡하여 다음에 소개할 수치적 방법에 비해 크게 장점이 없다.

Fig. 6. 
Calculated (a) extinction efficiencies (Q_ext), (b) asymmetry parameters (g), (c) scattering phase functions (P₁₁), and (d) log(P₁₁) of spherical liquid particles using the Lorenz-Mie code (Bohren and Huffman, 1983) at λ = 0.55 μm. The refractive index of liquid water at this wavelength is 1.333 + i1.96 × 10^-9. The approximated scattering ranges of the Rayleigh, Mie, and geometric optics (GO) are indicated in panels (a) and (b).

2.2.1.2 수치적 직접 계산법

다양한 수치적 직접 계산법이 존재하지만 이 방법들은 크게 미분 방정식법(differential equation method)과 적분 방정식법(integral equation method)의 두 분류로 구분할 수 있다. 미분 방정식법은 벡터 파동 방정식의 해를 주파수 또는 시간 영역에서 구하여 산란장을 계산하는 반면, 적분 방정식법은 맥스웰 방정식의 면적 적분 또는 체적 적분을 통해 해를 구한다.

2.2.1.2.1 미분 방정식법

미분 방정식법에는 변수 분리법, 유한차분 시간영역법[finite difference time domain (FDTD) method](Yee, 1966), 시간영역 유사 스펙트럴법[pseudo-spectral time domain (PSTD) method] (Liu, 1998), 유한 요소법[finite-element method (FEM)] (White, 1985), PMM (point-matching method) (Oguchi, 1973) 등이 있다. 이중 FDTD와 PSTD 방법을 비구형 구름입자 단일산란 특성 계산에 주로 사용한다.

2.2.1.2.1.1 유한차분 시간영역법(FDTD)

FDTD는 Yee (1966)가 제안한 미분 방정식법으로서 전자기 파동 방정식의 해를 진동수 영역에서 구하는 일반적인 수치적 방법과는 달리 시간 영역에서 맥스웰 방정식의 해를 구한다. 이때 일반적인 수치적 방법은 경곗값 문제(boundary-value problem)인 반면 FDTD는 초기값 문제(initial-value problem)다. 수학적으로 경곗값 문제가 초기값 문제보다 더 어렵다. FDTD는 전기장과 자기장의 시간과 공간 미분을 유한차분 방법을 사용하여 나타낸다. 이때 계산의 오차와 수치적 계산의 안정성을 유지하기 위해서 계산 영역을 적절한 시공간 크기로 나누는 작업이 필요하다. 산란체는 유한한 계산 영역 안에 놓이며 흡수 경계조건(absorbing boundary condition)을 적용하여 무한한 영역에 대해 계산을 수행한다. FDTD는 시간 영역에서 근접장에 대한 계산을 수행하기 때문에 근접장의 계산 결과를 주파수 영역에서의 원접장 산란장 결과로 변환해야 한다(Yang and Liou, 1996a). FDTD는 개념적으로 간단하고 산란체의 모양에 제약이 없다는 장점을 갖고 있다. 또한 불균질한 산란체에도 적용 가능하다. 하지만 산란체의 크기가 커짐에 따라 계산 시간이 많이 요구되며 계산 자원이 급속히 증가하는 단점이 있다. 또한 입사하는 빛에 대한 산란체의 방향의 변화가 있는 경우 새로운 계산을 수행해야 한다. 즉, 공간적으로 불규칙 방향(random orientation)을 갖는 빙정과 같은 산란체의 단일산란 특성 계산의 효율은 높지 않다.

2.2.1.2.1.2 시간영역 유사 스펙트럴법(PSTD)

PSTD는 Liu (1998)이 개발한 방법으로 FDTD에서 파생되었다. FDTD는 맥스웰 방정식의 공간 미분을 위해 유한차분 방법을 쓰는 반면 PSTD는 fast Fourier transform (FFT)을 사용한다. FDTD와 마찬가지로 산란체의 모양이나 균질성에 영향을 받지 않는 장점이 있으며 FDTD에 비해 계산 시간과 계산 자원 측면에서 더 효율적이기 때문에 비구형 빙정과 에어로졸 입자 단일산란 특성 계산에 많이 사용한다. 불규칙한 방향성을 갖는 산란체의 단일산란 특성 계산에 대해서는 FDTD와 마찬가지로 그 효율이 낮다. FDTD와 PSTD는 미국 Tex as A&M 대학의 Ping Yang 교수 연구팀에서 빙정과 에어로졸의 단일산란 특성 계산을 위해 활발히 사용되고 있다(Yang and Liou, 1996a; Yang et al., 2000, 2004; Chen et al., 2008; Liu et al., 2012; Panetta et al., 2013).

2.2.1.2.2 적분 방정식법

적분 방정식법에는 모멘트법[method of moment (MoM)]과 이산 쌍극자 근사[discrete dipole approximation (DDA)]로 대표되는 체적 적분 방정식법과 T-matrix로 대표되는 면적 적분법이 있다.

2.2.1.2.2.1 T-matrix

입사한 빛과 산란된 빛의 전기장은 적절한 형태의 벡터 파동 함수로 확장될 수 있다. 이때 산란된 파동의 확장계수(expansion coefficient)와 입사한 파동의 확장계수는 전이행렬(transition matrix)을 통해서 연관 지어지며 이 전이행렬이 바로 T-matrix이다. T-matrix는 입사한 빛의 파장에 대한 산란체의 모든 단일산란 정보를 담고 있으며 산란체의 크기, 모양, 굴절률(refractive index), 좌표계에 대한 상대적 방향에만 의존한다. 즉, T-matrix는 FDTD와 PSTD와는 다르게 입사하는 빛과 산란체의 상대적인 방향과는 무관하기 때문에 공간적으로 불규칙 방향을 갖는 산란체의 단일산란 특성 계산에 매우 효율적이다.

T-matrix는 Waterman (1971)에 의해 제안되었고 이후 Michael I. Mishchenko 박사가 개발한 코드(Mishchenko et al., 1996; Mishchenko and Travis, 1998)가 물리학, 천문학, 대기환경과학을 포함한 광범위한 연구분야에서 사용되고 있다(Mishchenko, 2020). Figures 7, 8은 Mishchenko T-matrix를 사용하여 계산한 회전 타원체(spheroid)의 산란 위상함수, 비대칭 인자, 소산계수를 보인다.

Fig. 7. 
Calculated scattering phase function (P₁₁), asymmetry parameter (g), and extinction efficiency (Q_ext) of spheroids with AR = 0.5 at λ = 0.55 μm. The P₁₁ calculated using T-matrix (blue) and discrete dipole approximation (DDA, red) are compared in panel (a) and (b), while those between DDA and geometric optics method (GOM, black) are shown in panel (c) and (d). Comparisons of g and Q_ext calculated using T-matrix, DDA, and GOM as a function of a volume-equivalent-sphere size parameter (χ_ve = 2πr_ve/λ) are shown in panel (e) and (f), respectively, where r_ve is a volume equivalent radius. The refractive index of ice at this wavelength is 1.311 + i2.289 × 10^-9.

Fig. 8. 
Same as Fig. 7, but AR = 2.0.

T-matrix는 분류상 NFM (null-field method) 또는 확장 경계 조건법[extended boundary condition method (EBCM)]에 속한다. 하지만 엄밀한 의미에서 T-matrix는 NFM이나 EBCM에 국한되지 않는다. NFM이나 EBCM은 T-matrix를 계산할 수 있는 여러 방법 중에 하나다. T-matrix는 변수 분리법을 사용해서도 계산 가능하며(Schulz et al., 1998) 이론적으로 모든 직접 계산법은 T-matrix를 계산할 수 있다(Markkanen and Yuffa, 2017).

일반적인 T-matrix의 단점은 회전대칭 모양(예, 원기둥, 회전타원체, Chebyshev particle 등)을 갖는 산란체의 단일산란 특성만을 계산할 수 있다는 것이다. 하지만 대부분의 빙정은 비구형이며 또한 회전 대칭체가 아니기 때문에 엄격한 기준에서는 T-matrix는 빙정 단일산란 특성 계산에 사용할 수 없다. 다른 방법들과 마찬가지로 T-matrix 또한 산란체의 크기 증가에 따른 계산 시간과 계산 자원의 급격한 증가가 필연적이며 산란체의 종횡비[aspect ratio (AR)]가 1에서 멀어짐에 따라 계산의 수렴 및 정확도에 문제가 발생한다(Mishchenko and Travis, 1998).

Mishchenko T-matrix로 대표되는 일반적인 방법 이외에도 중첩(superposition) T-matrix (Mackowski and Mishchenko, 1996, 2011; Mackowski, 2014; Markkanen and Yuffa, 2017), Tsym (Kahnert, 2013), invariant imbedding T-matrix method (II-TM) (Johnson, 1988; Bi et al., 2013a, 2013b; Hu et al., 2019; Shcherbakov, 2019)와 같은 변형된 T-matrix 방법이 존재한다. 중첩 T-matrix는 접하고 있는 구의 집합체의 단일산란 특성을 계산할 수 있다. 미국 Auburn 대학의 Daniel W. Mackowski 교수가 공개한 병렬화된 포트란 코드가 블랙카본 부착물(black carbon aggregates)의 단일산란 특성 계산에 널리 사용되고 있다(Liu and Mishchenko, 2007; Liu et al., 2008; Zhao and Ma, 2009; Kahnert and Devasthale, 2011; Wu et al., 2016; Liu et al., 2017, 2019). 더욱 발전된 형태의 중첩 T-matrix는 구가 아닌 불규칙한 모양을 갖는 입자들의 집합체에도 적용이 가능하다(Markkanen and Yuffa, 2017).

일반적인 T-matrix는 회전대칭 모양을 갖는 산란체에만 적용 가능하기 때문에 비구형 구름입자들에 적용할 수가 없다. Kahnert (2013)가 개발한 Tsym 코드는 일반적인 T-matrix 계산 방법을 변형하여 기하적 대칭을 갖는 산란체의 단일산란 특성 계산을 가능하게 한다. 예를 들어 빙정의 기본 모양인 육각형 기둥은 장축을 중심으로 60^o 회전을 한 이후에도 본래와 같은 모양을 갖는 기하적 대칭이다. 기존의 T-matrix 계산 방법과 비교하면 큰 발전이지만 Tsym은 매우 작은 크기의 산란체에만 적용 가능하며 계산 결과의 수렴을 수동으로 조절해야 하는 큰 단점이 있다.

최근 단일산란 특성 계산 연구 분야에서 가장 큰 발전을 꼽자면 II-TM 코드의 개발이라 할 수 있다(Bi et al., 2013a, 2013b). II-TM은 Johnson (1988)이 제안한 방법으로 산란체의 모양과 균질성에 관계없이 T-matrix를 계산할 수 있는 매우 큰 장점을 갖고 있다. 오랜 기간 학계에서 주목받지 못하던 II-TM을 Bi and Yang (2014)이 수치적 해를 개발하면서 그 우수성이 알려지게 됐다. II-TM은 산란체를 n개 층으로 이루어진 불균질한 구로 나누고 n번째 층의 T-matrix를 더 안쪽에 있는 n-1번째 층의 T-matrix로부터 구하는 계산과정을 반복한다. 산란체의 중심에 위치한 첫 번째 층의 T-matrix는 0이다. 최근의 연구들은 II-TM이 빙정 단일산란 계산에 높은 정확도와 유연성을 보임을 보고 하였다(Bi and Yang, 2014; Yang et al., 2019). 현시점에서 II-TM은 불규칙 방향을 갖는 비구형 빙정의 단일산란 특성 계산에 있어 가장 정확도가 높다. 하지만 다른 방법들과 마찬가지로 산란체의 크기가 커짐에 따라 계산 시간과 계산 자원이 급격히 증가하는 단점이 있다. II-TM 코드는 최근 다른 연구자들에 의해서도 개발되었으나(Hu et al., 2019; Shcherbakov, 2019) 현재 공개된 코드는 존재하지 않는다.

2.2.1.2.2.2 이산 쌍극자 근사(DDA)

체적 적분 방정식법의 한 종류인 DDA는 DeVoe (1964)에 의해 그 개념이 제안되었으며 Purcell and Pennypacker (1973)에 의해 발전되었다. DDA는 연결된 쌍극자 방법(coupled dipole method)으로도 불린다. DDA를 사용한 단일산란 계산에서는 산란체를 편광이 일어나는 N개의 쌍극자로 표현한다. 각 쌍극자는 주어진 위치에서 전기장에 반응하며 다른 쌍극자에도 영향을 준다. 하나의 쌍극자에 영향을 미치는 것은 입사하는 전기장과 다른 모든 쌍극자들이 만들어 낸 전기장이다. 즉, N개의 쌍극자에 영향을 미치는 N개의 전기장에 관한 N개의 선형 방정식이 생성된다. 이 N개의 선형 방정식은 직접적 또는 반복적 방법으로 풀 수 있다. 모든 DDA 코드는 반복적 방법[예, Fastest Fourier Transform in the West 3 (FFTW3)] (Frigo and Johnson, 2005)을 사용하고 있으며 DDA 계산에서 가장 시간이 많이 걸리는 부분이다. DDA의 계산 결과의 정확도에 있어 산란체의 모양을 표현하기 위해 사용한 쌍극자의 개수는 매우 중요하다. 쌍극자는 산란체에 비해 현격히 작아야 하며 그 조건은 d = λ/(10|m|)로 표현한다. 여기서 d와 m은 각각 쌍극자의 지름과 복소 굴절률이다. 쌍극자의 크기가 입사하는 빛의 파장의 최소 0.1배가 되어야 함을 의미하며 더 정확한 계산을 위해서는 쌍극자의 크기는 더 작아야 한다. 계산 가능한 굴절률의 범위는 |m - 1| < 2이다(Yurkin and Hoekstra, 2011). 0.55 μm 파장에서 빙정의 복소 굴절률은 1.311 + i2.289 × 10^-9이기 때문에 DDA를 빙정의 단일산란 특성 계산에 사용하기 적합하다.

DDA와 MoM을 포함하는 체적 적분 방정식법은 어떠한 모양의 산란체에도 적용할 수 있으며 산란체의 불균질성과 이방성(anisotropic)에도 영향을 받지 않는 큰 장점이 있다. 또한 T-matrix와는 다르게 계산의 안정성과 결과의 정확도가 산란체의 종횡비에 무관하다. 하지만 산란체의 크기가 커짐에 따라 계산 시간과 요구되는 계산 자원이 급격히 증가하며 입사하는 빛에 대한 산란체의 방향의 변화가 있는 경우 새로운 계산을 수행해야 하는 단점이 존재한다.

DDA는 다양한 학문분야에 사용되고 있으며 다수의 공개된 코드가 있다. Discrete Dipole SCATtering (DDSCAT) (Draine and Flatau, 1994)이 가장 역사가 길며 많이 사용되고 있다. 하지만 DDSCAT은 슈퍼컴퓨터와 같은 큰 계산 자원에서 병렬 계산을 수행 시 하나의 CPU에 산란체의 특정 방향만을 할당하는 제약이 있기 때문에 진정한 DDA의 장점을 살리고 있지 못하다. 반면 Amsterdam discrete dipole approximation (ADDA) (Yurkin and Hoekstra, 2007, 2011)은 전체 쌍극자의 개수를 나눠서 각 CPU에 배분하기 때문에 불규칙 방향을 갖는 비구형 빙정의 단일산란 특성 계산이 DDSCAT 보다 효율적이다. Figures 7, 8는 ADDA로 계산한 회전 타원체의 산란 위상함수, 위상행렬 요소, 비대칭 인자, 소산계수이고 Figs. 9, 10은 육각형 기둥과 판 모양의 빙정에 대한 계산 결과다. 그 외 OpenDDA (Donald et al., 2009), MATLAB package describing discrete dipole approximation (MPDDA) (Shabaninezhad et al., 2021) 등의 공개된 DDA 코드가 존재한다.

Fig. 9. 
Calculated non-zero phase matrix elements of hexagonal columns with aspect ratio (AR) of 0.5 at λ = 0.55 μm as functions of scattering angle and size parameter (χ_Dmax) calculated using a discrete dipole approximation. Non-zero phase matrix elements (a) P₁₁, (b) - P₁₂/P₁₁, (c) P₂₂/P₁₁, (d) P₃₃/P₁₁, (e) - P₄₃/P₁₁, and (f) P₄₄/P₁₁ are shown. Here, D_max used to define χ_Dmax is the longer length of L or W of a hexagonal column.

Fig. 10. 
Calculated scattering phase function [P₁₁, (a) and (b)], asymmetry parameter [g, (c) and (d)], and extinction efficiency [Q_ext, (e) and (f)] of hexagonal columnar ice crystals with AR = 0.5 (left panels) and 2.0 (right panels) using the DDA and GOM. In order to distinguish lines in panels (a) and (b) 1, 10, and 100 is multiplied to each of the original values from small to large sizes, respectively. The results of GOM and DDA are indicated with black dotted lines and color solid lines, respectively, in panels (a) and (b).

2.2.1에서 소개한 직접 계산법은 산란체의 크기(즉, 크기모수)가 증가함에 따라 계산 시간과 자원이 급격히 증가하기 때문에 크기가 큰 산란체의 단일산란 특성 계산에 효율적이지 않다. 이에 직접 계산법보다 정확도는 낮으나 계산이 상대적으로 빠르며 계산 자원이 적게 요구되는 Rayleigh approximation, RGA (Rayleigh-Gans approximation), ADA (anomalous diffraction approximation) (van de Hulst, 1981), 섭동론(perturbation theories), 기하 광학법[geometric optics method (GOM)] 등의 근사법이 개발되었다.

2.2.2.1 RGA와 ADA

RGA는 산란과 관련된 물리 과정을 단순화한 근사법이다(Bohren and Huffman, 1983). RGA는 산란체 내부에서 전자기장 사이의 고차 간섭을 무시하여 단일산란 계산을 간단하게 만든다. 산란체를 잘게 쪼갠체적 요소들은 입사하는 빛에 대해서만 반응하며 체적 요소에 의해 발생한 산란장은 다른 체적 요소에 영향을 미치지 않는다는 가정을 바탕으로 전체 산란장은 각 체적 요소들의 산란장을 중첩하여 계산한다. RGA의 적용 가능 범위 조건은 χ << 1과 |m - 1| << 1이다. 즉, 산란체는 광학적으로 부드러워야(soft) 하며 입사하는 빛의 파장보다 훨씬 작아야 한다(Mishchenko et al., 2002). 위의 조건들은 대기과학 분야에서 마이크로파와 부착과정으로 형성된 크기가 큰 빙정들의 단일산란 계산에 적용 가능하여 RGA는 레이다 기상학에서 많이 사용되고 있다(Hogan et al., 2012, 2017; Tyynelä et al., 2013; Hogan and Westbrook, 2014; Leinonen et al., 2017).

ADA는 van de Hulst (1981)이 개발한 방법으로 소산은 산란체 내부를 투과하는 빛의 흡수와 산란체를 투과하는 빛과 산란체 주변을 지나는 빛 사이의 간섭으로 이루어 진다는 가정을 하고 있다. ADA의 적용 가능 범위 조건은 χ >> 1과 |m - 1| << 1이다. ADA는 육각형 기둥과 같은 간단한 모양의 빙정의 단일산란 특성 계산에 사용되었으나 소산효율(extinction efficiency)과 소산단면(extinction cross-section) 등 일부 결과에 한정되었다(Chylek and Klett, 1991; Sun and Fu, 1999; Yang et al., 2004).

2.2.2.2 기하 광학법

기하 광학법은 ray optics 또는 ray tracing 방법으로도 알려진 근사법이다. 특히 가시광선 영역에서 비구형 빙정의 단일산란 특성 계산에 널리 사용되었다(Wendling et al., 1979; Cai and Liou, 1982; Muinonen et al., 1989; Takano and Liou, 1989; Macke, 1993; Iaquinta et al., 1995; Yang and Liou, 1995, 1996b, 1997, 1998; Macke et al., 1996a, b; Baran and Labonnote, 2007; Um and McFarquhar, 2007, 2009, 2011).

기하 광학법은 산란체의 크기가 입사하는 빛의 파장보다 매우 큰 경우(즉, χ >> 1)에 적용 가능하다. 이때 산란체에 입사하는 빛은 광선 다발(bundle of rays)로 간주한다. 기하 광학법에서 산란된 빛은 회절, 외부 반사, 내부 반사 이후 굴절 결과의 합으로 결정한다. 산란체 표면에서 반사 및 굴절된 각 광선의 경로와 에너지는 스넬(Snell)의 법칙과 프레스넬 계수(Fresnel coefficients)로 계산한다. 회절에 의한 효과는 산란체에 의한 회절이 산란체의 투영 단면(projected cross section)과 같은 모양과 크기를 갖는 구멍(aperture)에 의한 회절 효과와 동일 하다는 배비넷의 원리(Babinet’s principle)로 계산한다. 이때 프라운호퍼(Fraunhofer) 회절 계산과 원접장 가정을 적용한다(Jackson, 1998). 회절에 의한 빛의 에너지와 프레스넬 광선(반사 및 굴절된 빛)에 의한 에너지는 각각 입사한 빛의 방향으로 투영한 산란체의 단면에 입사한 빛의 에너지와 같다고 가정한다. 따라서 기하 광학법에서는 소산효율(Q_ext)이 크기모수와 관계없이 2이다(Figs. 7f, 8f, 10e-f). 원접장에서 산란된 빛의 강도는 주어진 입체각으로 산란되는 각 광선들의 강도의 합으로 결정한다. 산란된 광선들 간의 위상 속도 차이로 발생하는 간섭 효과는 산란체가 불규칙 방향을 갖는 경우 무시할 수 있다고 가정한다.

기하 광학법은 산란체의 모양에 관계없이 적용 가능하며 크기모수가 매우 큰 경우에도 단일산란 특성 계산이 빠르고 간단하다는 장점을 가지고 있다. 하지만 기하 광학법은 근사법이기 때문에 복수의 단점이 존재한다. 기하 광학법을 사용하여 평행 평면들을 포함하고 있는 산란체(예, 육각형 기둥과 판, 총알장미 등)의 산란 위상함수를 계산하는 경우 산란된 빛이 전방산란각에 집중되는 delta transmission 또는 delta-function transmission (Takano and Liou, 1989)으로 알려진 현상이 일어난다. 이는 실제 자연 현상이 아닌 기하 광학법이 물리적 광학 효과를 포함하고 있지 않기 때문에 발생한 인위적인 현상이다. 같은 현상이 외부 반사한 빛에서도 일어난다(Mishchenko and Macke, 1998). 또한 일반적인 기하 광학법은 광선들의 위상 차이에 의한 간섭 효과를 고려하고 있지 않으며 이는 후방 산란된 빛의 강도에 영향을 미칠 수 있다(Mishchenko et al., 2002). 크기모수가 작은 산란체에 대한 기하 광학법의 적용 가능 범위는 정확히 알려져 있지 않기 때문에 이는 직접 계산법과의 비교를 통하여 신중히 결정해야 한다. Um and McFarquhar (2015)는 육각형 기둥과 판 모양의 빙정에 대한 기하 광학법 적용 가능 범위는 산란체의 종횡비에 따라 다르다는 것을 보였다. 일반적인 기하 광학법의 약점을 극복하기 위해 면적 적분 또는 체적 적분 방정식을 사용하여 근접장 계산 결과를 원접장에 사상(mapping)하는 physical geometric optics (Kirchhoff) method (PGOM) (Muinonen, 1989; Bi et al., 2011; Sun et al., 2017)와 PGOM을 단순화한 형태인 improved geometric optics method (IGOM) (Yang and Liou, 1996b)가 개발되었다.

공개된 기하 광학법 코드로는 독일 Leipzig 대학의 Andreas Macke 교수의 ray tracing 코드(Macke 1993; Macke et al., 1996a, 1996b))가 있다. Andreas Macke 교수의 코드는 여러 코드의 집합체로서 원기둥, 구, 회전 타원체, 육각형 기둥과 판, 총알, 총알장미 모양 등의 산란체에 대해 계산할 수 있으며 광범위한 분야에서 사용되고 있다. Figures 7과 8은 Andreas Macke 교수의 코드로 계산한 회전 타원체의 산란 위상함수, 위상행렬 요소, 비대칭 인자, 소산계수이고, Figs. 10과 11은 육각형 기둥과 판, 총알장미 모양의 빙정에 대한 계산 결과다.

Fig. 11. 
Calculated non-zero phase matrix elements of bullet rosette (BR, red), hexagonal column (Col, blue), and hexagonal plate (Pla, black) at λ = 0.55. Non-zero phase matrix elements (a) P₁₁, (b) - P₁₂/P₁₁, (c) P₂₂/P₁₁, (d) P₃₃/P₁₁, (e) - P₄₃/P₁₁, and (f) P₄₄/P₁₁ are shown. A hexagonal column with L = 200 μm and AR = 2.0, a hexagonal plate with W= 200 μm and AR=0.5, and a bullet rosette consisting of six bullets with L = 100 μm and AR = 2.0 are used. The corresponding calculated asymmetry parameters (g) are also embedded in panel (a).

또 다른 기하 광학법을 사용하는 공개된 코드로는 준구형 모양을 갖는 작은 산란체(Gaussian random sphere)의 단일산란 특성을 계산할 수 있는 핀란드 헬싱키 대학의 Karri Muinonen 교수의 SIRIS 코드(Muinonen et al., 1996)가 있다. Gaussian random sphere 모형은 작은 빙정의 단일산란 특성 계산에 사용되었다(Nousiainen and McFarquhar, 2004; Nousiainen et al., 2011; Um and McFarquhar, 2011).

Figures 7과 8은 불규칙 방향을 갖는 회전 타원체의 단일산란 특성을 T-matrix, DDA, GOM으로 계산한 결과이다. 수치적 직접 계산법인 T-matrix와 DDA로 계산한 산란된 빛의 강도를 산란각의 함수로 표현하는 산란 위상함수(P₁₁)의 계산 결과는 종횡비에 관계없이 두 방법이 잘 일치하고 있음을 보인다(Figs. 7a-b, 8a-b). 하지만 근사법인 GOM과 직접 계산법인 DDA의 계산 결과는 뚜렷한 차이를 보인다(Figs. 7c-d, 8c-d). 이 차이는 산란체의 크기가 증가함에 따라 줄어든다(Figs. 7d, 8d). 근사법인 GOM의 계산 결과의 정확도는 산란체의 크기가 증가함에 따라 높아지는 것을 비대칭 인자(Figs. 7e, 8e)와 소산효율(Figs. 7f, 8f) 계산 결과에서 알 수 있다.

Figure 10은 빙정의 기본 모양인 육각형 기둥(종횡비 2.0)과 판(종횡비 0.5)의 산란 위상함수, 비대칭 인자, 소산효율을 DDA와 GOM으로 계산한 결과다. 이 계산에서 빙정은 불규칙 방향을 갖는다. 빙정의 크기가 커짐에 따라 DDA 산란 위상함수 계산 결과는 육각형 빙정의 고유 특성인 22^o 와 46^o 무리와 130^o 와 165^o 사이에 만들어지는 얼음활(ice bow) 현상을 뚜렷이 나타낸다(Figs. 10a, b). 비대칭 인자(Figs. 10c, d)와 소산효율(Figs. 10e, f) 계산 결과는 근사법인 GOM의 적용 가능 범위를 보여주며, 이는 빙정의 종횡비 0.5와 2.0에 따라 다르다. Um and McFarquhar (2015)는 일반적인 GOM의 단일산란 특성 계산의 정확도는 빙정의 크기가 커짐에 따라 증가하며 빙정의 종횡비가 1에 가까울수록 더 작은 크기모수에서 DDA 계산 결과와의 일치 정도가 더 높은 것을 보였다. 또한 GOM으로 계산한 비대칭 인자와 소산효율의 정확도는 체적등가구 크기모수(volume-equivalent-sphere size parameter, χ_ve) 90에 대해서 각각 약 1.2%와 7.0%이며, 이는 체적등가구 크기모수 100에 대해서는 각각 약 0.8%와 3.3%임을 보였다. 위의 결과는 GOM과 같은 근사법의 적용범위 하한의 중요성을 강조한다.

Figures 6a와 6b는 레일리(Rayleigh) 산란, 미 산란, 기하광학의 세가지 산란 영역을 근사적으로 크기모수의 함수로 구분한다. 레일리 산란 영역은 산란효율이 크기모수의 증가에 따라 대략 선형적으로 증가하는 영역이다. 기상 레이다 연구에서 사용하는 대기 수상체(hydrometeor)와 마이크로파의 관계가 이 영역에 속하며 비교적 적은 계산 자원과 DDA와 같은 직접법 또는 RGA와 같은 근사법을 사용하여 단일산란 특성을 계산할 수 있다. 미 산란 영역은 산란효율이 크기 모수의 증가에 따라 진동하는 특성이 있다. 이는 회절한 빛과 투과된(굴절된 빛과 반사된 빛) 빛 사이의 간섭의 결과이다. 미 산란 영역에서 보이는 특징적인 산란효율의 진동은 크기모수가 커짐에 따라 진폭이 감소하며 가시광선 영역에서 2로 수렴하여 기하광학 영역으로 접어들게 된다. 가시광선 영역에서 빙정의 단일산란 계산의 가장 큰 문제는 미 산란 영역과 기하광학 영역이 중첩되는 부분에서 발생한다.

Figure 12는 구름입자 단일산란 계산에서의 어려움을 설명한다. 구형의 액체상 구름입자는 크기모수와 관계없이 로렌츠-미 이론의 적용이 가능하다. 로렌츠-미 이론은 계산 자원이 많이 요구되지 않으며 정확도 또한 매우 높은 직접적인 방법이다. 비구형 빙정의 경우 수치적으로 정확한 직접적 방법(예, DDA, FDTD, PSTD, T-matrix, II-TM 등)과 근사법(예, GOM, IGOM, PGOM 등)의 적용 가능 범위가 겹치는 회색영역(gray zone)이 발생한다. 이 영역에서 직접적 방법은 여전히 정확한 결과를 제공하지만 필요한 계산 자원과 시간이 기하급수적으로 증가하여 단일 비구형 빙정의 단일산란 특성 계산을 하기 위해서는 반드시 슈퍼컴퓨터가 필요하다. 반면 GOM과 같은 근사법은 매우 빠른 계산을 수행할 수 있으나 회색영역에서는 그 정확도가 낮다. 근사법의 적용가능 범위의 하한은 반드시 직접적 방법을 사용하여 계산한 정확한 결과와 비교해야 한다(Mishchenko et al., 2002). 또한 이 적용가능 범위의 하한은 비구형 빙정의 모양에 따라 다르다(Um and McFarquhar, 2015).

Fig. 12. 
Schematic diagram of numerical methods to calculate the single-scattering properties of spherical (i.e., Mie) and nonspherical cloud particles. There is a gray zone where the accuracy of a geometric optics method (red) is low and the computational cost for exact methods (blue) is very expensive for the calculations of single-scattering properties of nonspherical ice crystals. The Lorenz-Mie theory can be applied to spherical particles regardless of a size parameter without losing the accuracy.

현 시점에서 공개된 코드를 사용하여 회색영역을 포함한 모든 크기모수의 빙정의 단일산란 계산은 DDA와 GOM의 조합으로 가능하다. 하지만 회색영역에서의 계산은 주의를 기울여야 하며 많은 계산 자원과 시간을 필요로 한다. 공개하지 않은 코드를 사용하여 계산을 수행한다면 II-TM과 PGOM의 조합으로 더 효율적으로 수행 가능하다. 이미 Tex as A&M 대학의 Ping Yang 교수 그룹에서 이 계산 결과를 산란 도서관 자료로 만들어 놓았다(Yang et al., 2005, 2013). 이 도서관 자료는 대출만이 가능하며 다른 이들에 의한 수정 또는 추가는 불가능하다. 빙정의 미세물리 측면에서 이 도서관 자료를 평가하자면 장서가 매우 부족하다. 이는 다음에 설명할 빙정 미세물리와 매우 밀접하게 연관되어 있다. 하지만 위의 도서관 자료가 단파와 장파복사 파장대를 모두 아우르는 현존하는 거의 유일한 도서관 자료임에는 틀림이 없다.

식(10)~(12)에서 필요한 구름입자 미세물리 정보는 크기분포와 모양분포다. 산란단면 계산을 위한 구름입자들의 투영 단면적이 필요하지만 이는 모양분포 정보에서 추출 가능하다. 이런 구름입자 미세물리 정보는 주로 항공기를 사용한 관측을 통해서 획득한다.

Figure 1은 2.3 μm 해상도의 Cloud Particle Imager (CPI) (Lawson et al., 2019)로 관측한 구름 안에 존재하는 비구형 빙정과 구형 액체상 구름입자다. 액체상구름의 경우 식(10)~(12)에서 모양 빈을 고려할 필요가 없지만 얼음상 구름의 경우 빙정의 모양이 다양하며 빙정의 단일산란 특성이 모양에 따라 크게 달라지기 때문에 이를 반드시 고려해야 한다. Figure 1에 분류된 구름입자들의 모양을 바탕으로 Fig. 3과 같은 모양분포를 생산하고(Um and McFarquhar, 2009) 이를 Fig. 4의 관측한 크기분포와 결합하여 식(10)~(12)에 적용한다.

실험실에서 얼음상 구름의 환경 조건과 유사한 환경에서 성장한 빙정(Pfalzgraff et al., 2010; Magee et al., 2014)과 풍선을 사용하여 채집한 빙정(Magee et al., 2021)을 전자 현미경으로 관찰한 결과는 빙정의 표면이 다양한 형태의 마이크로미터 이하(sub-micron) 규모의 다양한 표면 거칠기(surface roughness)를 갖고 있음을 보였다. 또한 마이크로미터 이상(super-micron) 규모의 이상적인 모습으로부터 벗어나는 불규칙적인 내부와 외부 구조[예, 공동(hollowness)] 및 내부 첨가물(예, 공기방울, 황산암모늄, 검댕)의 발견은 실제 빙정과 단일산란 계산에 사용하는 이상화된 빙정 모형에 심각한 괴리가 있음을 보였다. Figure 13은 environmental scanning electron microscope (ESEM)을 사용하여 빙정의 다양한 성장 단계(예, 확산 성장과 승화)를 나노 미터 단위 이하의 고해상도로 관측한 결과다(Magee et al., 2014). 빙정은 일부 매끄러운 표면(Figs. 13a, d)으로 이루어져 있지만 더 많은 부분, 특히 승화에 의한 소멸 단계(Figs. 13b, e)에서 거친 표면과 불규칙한 외형을 보인다.

Fig. 13. 
Micrographs of ice crystals acquired by high resolution environmental scanning electron microscope (ESEM). Panels (a) and (d) show growing ice crystals, panels (b) and (e) show sublimating ice crystals, and panels (c) and (f) show ice crystals grown in an external diffusion chamber and transferred into the ESEM for imaging at equilibrium conditions.

Figure 13에 보인 관측된 빙정의 형태학적 복잡성(morphological complexities)은 단일산란 계산에 사용하는 이상화된 빙정 모형과는 매우 다르며, 이는 관측한 실제 구름의 복사적 특성과 이론적 계산의 큰 차이를 야기하는 원인으로 지목되었다(Sassen et al., 1994; Korolev et al., 1999). 이 차이를 극복하기 위해 마이크로미터 이하 규모의 다양한 빙정의 표면 거칠기를 이상화된 빙정 모형 표면에 직접적(Yang and Liou, 1998; Shcherbakov et al., 2006; Liu et al., 2013; Neshyba et al., 2013; Collier et al., 2016; Zhang et al., 2016)으로 또는 기하 광학법에 간접적(Macke et al., 1996a; Hess et al., 1998; Um and McFarquhar, 2007, 2009; van Diedenhoven et al., 2014a, 2014b; Geogdzhayev and van Diedenhoven, 2016)으로 표현하여 빙정의 단일산란 특성을 계산하였다. Macke et al. (1996a)는 기하 광학법을 사용하여 단일산란 특성 계산 시 반사 또는 굴절된 광선의 경로를 미리 정해놓은 범위 내에서 통계적인 왜곡(distortion)을 적용하는 방법을 개발하였다. Figure 14는 이 광선 경로의 통계적 왜곡을 적용하여 계산한 육각형 기둥, 판, 총알장미 모양의 빙정의 산란 위상함수와 비대칭 인자다. 광선 왜곡모수(ray-distortion parameter)가 증가함에 따라 산란 위상함수는 22^o와 46^o 무리와 같은 특성이 사라지면서 밋밋한 형태로 변하며 비대칭 인자는 감소한다.

Fig. 14. 
Calculated scattering phase function (P₁₁) of (a) column, (b) plate, and (c) bullet rosette using a geometric optics method with varying distortion parameter (t) of 0.0 (black), 0.1 (green), 0.2 (blue), and 0.3 (red). The corresponding calculated asymmetry parameters (g) are also embedded in each panel.

많은 연구에서 Fig. 13과 같은 빙정의 형태학적 복잡성은 빙정의 단일산란 특성에 큰 영향을 미침을 이론적 계산을 통해 보였다(Macke et al., 1996b; C.-Labonnote et al., 2001; Xie et al., 2009; Bi and Yang, 2014; Panetta et al., 2016; Smith et al., 2016). 이 연구들의 공통된 결과는 실제 관측 값에 더 가까운 특징 없는 산란 위상함수와 낮은 비대칭 인자 값(e.g., Fig. 14)이다(Knap et al., 1999; Yang et al., 2008; Baum et al., 2010; van Diedenhoven et al., 2013; Cole et al., 2014). 또한 빙정의 형태학적 복잡성을 고려하면 위성 관측 결과와 더 큰 일치를 보였다(Shcherbakov et al., 2006; Baum et al., 2011; Cole et al., 2014). 하지만 그 일치 정도는 비교에 사용한 구름에 따라 매우 많이 달라진다(Hioki et al., 2016). 큰 규모의 수치 모형에서 빙정의 표면 거칠기를 고려한 구름의 복사적 특성 모수화를 개발한 연구는 많지 않지만(Edwards et al., 2007; Yi et al., 2013; van Diedenhoven et al., 2014a, 2016; Smith et al., 2016; Järvinen et al., 2018)이 연구들은 빙정의 형태학적 복잡성이 중요함을 보였다. Yi et al. (2013)은 빙정의 표면 거칠기의 변화가 단파장에서 전구 평균 1~2W m^-2의 구름복사강제력과 10W m^-2 이상의 지역적 효과를 야기함을 보였다. 또한 Järvinen et al. (2018)은 빙정의 표면 거칠기에 의한 단파장에서 전구 규모의 1.12W m^-2의 냉각 효과를 보였다.

Järvinen et al. (2018)은 실제 관측한 빙정의 평균적인 산란 위상함수와 가장 큰 일치도를 보이는 이상화된 빙정 모형은 severely roughened hexagonal aggregates (SRHA) (Yang and Liou, 1998)임을 보였다. SRHA 빙정 모형은 현재 MODIS Collection 6 얼음상 구름 산출 알고리즘에서 사용하고 있는 모형이다(Platnick et al., 2017). SRHA 빙정 모형은 여섯 개의 육각형 기둥의 부착으로 이루어져 있으며 고정된 매우 거친 표면을 갖고 있다. 이 모형을 사용한 계산 결과는 특징 없는 밋밋한 산란 위상함수와 약 0.75의 작은 비대칭 인자 값을 보인다.

표면의 거칠기를 고려한 SRHA 빙정 모형의 개발 등의 노력에도 불구하고 현재 빙정의 형태학적 복잡성은 더욱 개선돼야 한다. 이는 빙정의 형태학적 복잡성이 주변 환경(예, 온도, 습도, 압력 등)의 변화에 따라 그 변화의 폭이 매우 넓고(Pfalzgraff et al., 2010; Magee et al., 2014, 2021) 현재 사용하는 빙정의 형태학적 복잡성을 단일산란 특성 계산에 적용하는 방법은 온도나 습도와는 무관한 고정된 거칠기 값을 사용하기 때문이다. 최근 수행된 실험의 결과들은 현재 사용하고 있는 빙정 표면 거칠기 모형이 실제 관측한 빙정의 형태학적 복잡성을 적절히 반영하고 있지 못함을 보였다(Magee et al., 2014; Schnaiter et al., 2012, 2016; Ulanowski et al., 2014). 즉, 현재의 이상화된 빙정 모형은 구름 미세물리 과정을 반영하고 있지 않다. 이는 단일산란 특성 계산을 위한 이상화된 빙정 모형이 구름 미세물리 과정을 고려하지 않고 관측된 특징 없는 산란 위상함수와 낮은 비대칭 인자 값을 맞추는 것을 목표로만 개발되었기 때문이다. 구름 미세물리적 특성과 복사적 특성 간에 큰 불일치가 존재하는 것이다. 예를 들어 SRHA 빙정 모형은 여섯 개의 육각형 기둥의 부착으로 이루어진 고정된 형태이기 때문에 구름 미세물리 과정에서 중요한 크기분포와 연관하여 표현하는 얼음수함량(ice water content)과 종단속도(terminal velocity) 등을 표현함에 있어 큰 문제를 야기한다. 구름의 복사적 측면에서는 SRHA 빙정 모형은 우수하지만 구름의 미세물리 과정에서는 사용할 수 없는 모형이다. 엄밀한 의미에서 SRHA 빙정 모형은 3.1에서 설명한 구름 전체의 복사적 특성을 계산하는 세 단계 중에 첫 번째 단계를 생략하고 두 번째 단계를 수행하는 것이다. 단일 빙정 입자의 단일산란 특성을 계산하는 것이 아니라 구름의 평균적인 단일산란 특성에 가장 가까운 결과를 만들어 내는 것이기 때문이다. SRHA 빙정 모형의 단점을 극복하기 위해 Voronoi aggregates 모형(Ishmoto et al., 2012)과 two-habit 모형(Liu et al., 2014) 등이 개발됐다. 이 모형들은 실제 빙정의 모양과는 형태학적으로 다르나 항공기 관측을 통해 수집된 얼음상 구름의 얼음수함량을 유사하게 표현할 수 있기 때문에 SRHA 빙정 모형과 비교 시 구름 미세물리 측면에서 더 의미가 있다. 하지만 이 모형들도 실제 빙정과는 다른 형태를 갖고 있기 때문에 여전히 부착, 결착, 강수 입자로의 전환 등을 표현하는 구름 미세물리 과정의 모수화에는 사용할 수 없다.

현재의 단일산란 계산에 사용하는 이상화된 모형과 구름 미세물리 과정에는 큰 괴리가 존재한다(Baran, 2012; van Didenhoven et al., 2014b; Senf and Deneke, 2017). 이는 두 분야가 독립적으로 발전하였기 때문이다. 또한 빙정의 성장과정에 따른 형태학적 복잡성을 주변 환경변수 또는 구름의 성장 단계에 따라 정량화하지 않았기 때문이다(Pfalzgraff et al., 2010; Magee et al., 2014, 2021). 긴 역사의 단일산란 연구분야에 비하여 상대적으로 역사가 짧은 구름 물리학 분야의 발전 속도가 더디다. 특히 환경 조건에 따른 빙정의 미세한 형태학적 복잡성 변화를 정량화 할 수 있는 관측이 부족한 것이 가장 큰 문제다. 또한 빙정의 거시적인 미세물리 특성(예, 크기분포와 모양 분포)의 다양한 환경 조건에 따른 변화를 특징 할 수 있는 충분한 전구 규모의 관측 자료가 생산되어야 한다. 구름 미세물리와 복사적 측면에서 물리적으로 일관된 빙정 모형 개발을 위해서는 두 분야의 연결 고리가 되는 항공기 및 구름 챔버 등의 실험실의 관측 자료가 환경조건에의 함수로서 충분히 확보되어야 한다.