地磁場(chǎng)發(fā)電機(jī)數(shù)值模擬綜述*

董 超，張 懷，石耀霖

(中國(guó)科學(xué)院大學(xué)地球科學(xué)學(xué)院 中國(guó)科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100049)

地球磁場(chǎng)由液態(tài)外核中的磁流體動(dòng)力學(xué)(magnetohydrodynamics, MHD)過(guò)程所產(chǎn)生和維持，這一機(jī)制被稱為地磁場(chǎng)發(fā)電機(jī)[1]。熱浮力和組份浮力為這一過(guò)程提供能量，驅(qū)動(dòng)導(dǎo)電流體形成對(duì)流。發(fā)電機(jī)數(shù)值模擬意在通過(guò)數(shù)值計(jì)算，得到自持式發(fā)電機(jī)過(guò)程[2]。從1950年代到1990年代，發(fā)電機(jī)理論已經(jīng)證實(shí)外核發(fā)電機(jī)過(guò)程的可行性和基本條件。然而，只有當(dāng)MHD數(shù)值計(jì)算足夠成熟后，地磁場(chǎng)發(fā)電機(jī)的數(shù)值結(jié)果才可以展示所觀測(cè)的地磁場(chǎng)詳細(xì)性質(zhì)[3]。

隨著MHD數(shù)值模擬技術(shù)的成熟與大規(guī)模高性能計(jì)算機(jī)的飛速發(fā)展，發(fā)電機(jī)模擬在1995年迎來(lái)了突破。Glatzmaier 和Roberts[4]的發(fā)電機(jī)模型GR95，不僅得到了三維自持發(fā)電機(jī)解，也展現(xiàn)了地磁場(chǎng)的一些基本形態(tài)特征和物理過(guò)程，例如偶極占優(yōu)、地磁西漂和磁極倒轉(zhuǎn)。此外，GR95模型還研究?jī)?nèi)外核黏性耦合以及電磁耦合給內(nèi)核差動(dòng)旋轉(zhuǎn)帶來(lái)的影響[5]，并被隨后的地震波研究所證實(shí)[6]。其后陸續(xù)有新的發(fā)電機(jī)模型被建立，如Kageyama等; Kuang 和 Bloxham; Christensen 等; Fearn; Dormy 等; Jones等; Busse; Zhang 和 Schubert; Kono 和 Roberts; Glatzmaier 等和Takahashi等等[7-10]。大多數(shù)模型可以近似吻合地磁場(chǎng)的空間譜、長(zhǎng)期變化以及地磁場(chǎng)倒轉(zhuǎn)等現(xiàn)象。近年來(lái)，科學(xué)家已經(jīng)能夠逐漸運(yùn)用發(fā)電機(jī)數(shù)值模擬結(jié)果解釋地球、行星[11]和恒星[12-13]磁場(chǎng)的一些特殊性質(zhì)與機(jī)制[14]。

然而，由于發(fā)電機(jī)模型參數(shù)還沒(méi)有達(dá)到完全類地行為所需的極端分辨率，所以這些數(shù)值模型還不能解釋關(guān)于地磁場(chǎng)發(fā)電機(jī)是如何運(yùn)作的很多基本問(wèn)題。特別是，期望外核的流場(chǎng)展現(xiàn)大范圍的長(zhǎng)度特征尺度，從黏性邊界層的厚度(～0.1 m)到地核半徑(～7×106m)，外加一個(gè)相當(dāng)完整的時(shí)間范圍，都是現(xiàn)在模型無(wú)法實(shí)現(xiàn)的[15]。研究人員通過(guò)采用比實(shí)際地核值大很多的黏性擴(kuò)散和熱擴(kuò)散值來(lái)解決這個(gè)問(wèn)題，例如假設(shè)超大的黏性值用以壓制小尺度的湍流運(yùn)動(dòng)。但是對(duì)小尺度湍流的抑制是否會(huì)嚴(yán)重影響發(fā)電機(jī)過(guò)程的物理本質(zhì)仍然是一個(gè)懸而未決的問(wèn)題[3]。

本文簡(jiǎn)述地磁場(chǎng)發(fā)電機(jī)基準(zhǔn)模型的方程組、無(wú)量綱方案、初始條件邊界條件、數(shù)值方法和標(biāo)度律等基本問(wèn)題，并且以MoSST模型為例展示數(shù)值模擬結(jié)果，在文章的最后展望地磁場(chǎng)發(fā)電機(jī)數(shù)值模擬的未來(lái)發(fā)展。

1 地磁場(chǎng)發(fā)電機(jī)模型基礎(chǔ)

一個(gè)成功的自持式地磁場(chǎng)發(fā)電機(jī)模型需要滿足一些基本的條件：1)一個(gè)合理的動(dòng)力機(jī)制來(lái)驅(qū)動(dòng)導(dǎo)電流體運(yùn)動(dòng)。大部分模型假設(shè)的驅(qū)動(dòng)力是熱動(dòng)力和/或組份對(duì)流，少部分模型假設(shè)的驅(qū)動(dòng)力是地球進(jìn)動(dòng)[16]；2)Cowling反發(fā)電機(jī)理論證明，軸對(duì)稱磁場(chǎng)不可能由軸對(duì)稱運(yùn)動(dòng)所維持[17]。因此地磁場(chǎng)發(fā)電機(jī)模型中的流場(chǎng)和磁場(chǎng)必須都是三維的，球形的幾何結(jié)構(gòu)也是再現(xiàn)地球全球偶極磁場(chǎng)占優(yōu)的根本要素；3)模型必須是旋轉(zhuǎn)的。因?yàn)榭剖狭κ菢?gòu)建流場(chǎng)結(jié)構(gòu)和產(chǎn)生全球偶極磁場(chǎng)占優(yōu)的關(guān)鍵因素[1]。

相比于以上因素，關(guān)于地核其他性質(zhì)的假設(shè)就沒(méi)有那么重要了。在大多數(shù)模型中，忽略密度的可壓縮性和地核的徑向熱結(jié)構(gòu)。取而代之的是使用布辛涅斯克近似，溫度和組份的改變只通過(guò)Navier-Stokes方程(N-S方程)的浮力項(xiàng)輸入，進(jìn)而改變密度[18]。從這個(gè)角度看，地磁場(chǎng)發(fā)電機(jī)的數(shù)值模型和氣體行星或恒星發(fā)電機(jī)數(shù)值模型有著較大的差異[19]。大多數(shù)的地磁場(chǎng)發(fā)電機(jī)模型的物質(zhì)屬性都是常數(shù)，外核是一個(gè)可旋轉(zhuǎn)的液態(tài)球殼。在球殼內(nèi)部，內(nèi)核的電導(dǎo)率和外核一樣(有些模型把內(nèi)核簡(jiǎn)化為絕緣體)，在球殼外部是絕緣的地幔。

根據(jù)以上這些方法可以構(gòu)建一種最簡(jiǎn)單的模型，把它定義為地磁場(chǎng)發(fā)電機(jī)數(shù)值模型的基準(zhǔn)模型[20]。在接下來(lái)的章節(jié)，我們都是以基準(zhǔn)模型為參考，討論各種發(fā)電機(jī)模型的不同修改。

1.1 地磁場(chǎng)發(fā)電機(jī)方程組

球坐標(biāo)下發(fā)電機(jī)模型的控制方程組由動(dòng)量方程(Navier-Stokes方程)、磁感應(yīng)方程、熱對(duì)流方程以及速度場(chǎng)和磁場(chǎng)的無(wú)散約束構(gòu)成。其中，N-S方程包含科氏力項(xiàng)和洛倫茲力項(xiàng)；由麥克斯韋方程組推出的磁感應(yīng)方程中忽略位移電流[21-22]：

(1)

一些模型忽略N-S方程中的慣性項(xiàng)[4-5]，或做簡(jiǎn)化處理只保留軸對(duì)稱分量[8,23-24]。這樣處理除簡(jiǎn)化的因素外，是由于對(duì)地球而言，慣性力相比于科氏力和洛倫茲力要小得多。然而沒(méi)有任何一個(gè)模型直接比較有或沒(méi)有慣性力項(xiàng)的情況，在地核條件下慣性力也許起的作用很小，這仍是一個(gè)未知問(wèn)題。

大多數(shù)發(fā)電機(jī)數(shù)值模型假設(shè)流體不可壓縮并采用Boussinesq近似，然而Glatzmaier和Roberts[23]及隨后的一些模型考慮了密度分層結(jié)構(gòu)和其他非Boussinesq近似的影響，如擴(kuò)散率對(duì)熱方程的貢獻(xiàn)。在這些模型中，使用滯彈性假設(shè)，不同于通過(guò)壓制聲波的完全可壓縮形式。Gastine等[25]對(duì)此進(jìn)行系統(tǒng)的研究，結(jié)果表明當(dāng)層間密度差超過(guò)4倍時(shí)，磁場(chǎng)將會(huì)變?yōu)榉桥紭O場(chǎng)。如果地球外核是密度差為1.2倍的弱密度分層結(jié)構(gòu)，則可能會(huì)對(duì)磁場(chǎng)形態(tài)和發(fā)電機(jī)系統(tǒng)有重要的影響[26-27]。

1.2 無(wú)量綱方案

大部分模型對(duì)方程組(1)進(jìn)行無(wú)量綱化處理，然而無(wú)量綱化方案有很多種，所以不同的地磁場(chǎng)發(fā)電機(jī)模型使用不同的無(wú)量綱參數(shù)[4-5,8,24,27]。對(duì)于對(duì)流驅(qū)動(dòng)型發(fā)電機(jī)模型來(lái)說(shuō)，基本有4個(gè)獨(dú)立的無(wú)量綱參數(shù)。在地磁場(chǎng)發(fā)電機(jī)基準(zhǔn)模型的定義中，以D=ro-ri為長(zhǎng)度特征尺度(ro為外核半徑，ri為內(nèi)核半徑，基準(zhǔn)模型中ri/ro=0.35)，黏性擴(kuò)散時(shí)間t=D2/ν為時(shí)間特征尺度，B=(ρμ0ηΩ)1/2為磁場(chǎng)特征尺度，以內(nèi)外核邊界的溫度差ΔT為溫度特征尺度，將方程組(1)無(wú)量綱化，得到地磁場(chǎng)發(fā)電機(jī)模型的無(wú)量綱化方程組：

(2)

式中：?jiǎn)挝幌蛄縵為旋轉(zhuǎn)軸的方向，假設(shè)重力與半徑為線性的關(guān)系g(r)=gor/ro，go為外核邊界(CMB)的重力值。式中的無(wú)量綱參數(shù)定義如表1所示。

表1 地磁場(chǎng)發(fā)電機(jī)數(shù)值基準(zhǔn)模型各無(wú)量綱參數(shù)定義Table 1 Definitions of dimensionless parameters in geodynamo benchmark

注：α為外核流體熱膨脹系數(shù)[3,15]。

表1中的外核估值根據(jù)參數(shù)定義計(jì)算得到，所用到的擴(kuò)散系數(shù)等物理參數(shù)的估值可參考Olsen等[28]、Christensen和Wicht[3]的研究結(jié)果。從表1可以看出，這些參數(shù)在數(shù)值模擬中的取值與實(shí)際值仍存在較大差距，目前的計(jì)算能力還遠(yuǎn)遠(yuǎn)無(wú)法達(dá)到參數(shù)的估計(jì)數(shù)量級(jí)。

無(wú)量綱方程組和無(wú)量綱參數(shù)的不同是由于所選取的特征尺度不同造成的。對(duì)于長(zhǎng)度特征尺度，部分模型選擇ro代替D。對(duì)于時(shí)間特征尺度，部分模型選擇磁擴(kuò)散時(shí)間或者熱擴(kuò)散時(shí)間來(lái)代替黏性擴(kuò)散時(shí)間。同樣對(duì)于磁場(chǎng)特征尺度和溫度場(chǎng)特征尺度的定義也有差異[3]。由此就會(huì)導(dǎo)致不同的模型定義不同的無(wú)量綱參數(shù)，從而造成數(shù)值計(jì)算的復(fù)雜及多樣性，但是通常來(lái)說(shuō)其他無(wú)量綱參數(shù)都可以由基準(zhǔn)模型中的4個(gè)無(wú)量綱參數(shù)表示[10]。

地磁場(chǎng)發(fā)電機(jī)模型的數(shù)值解也需要用無(wú)量綱診斷參數(shù)來(lái)表示。比較常見(jiàn)的幾個(gè)無(wú)量綱診斷參數(shù)為磁Reynolds數(shù)Rm，Reynolds數(shù)Re，Rossby數(shù)Ro，Elsasser數(shù)Λ，其定義如表2所示。

表2 地磁場(chǎng)發(fā)電機(jī)數(shù)值基準(zhǔn)模型各無(wú)量綱診斷參數(shù)定義Table 2 Definitions of dimensionless diagnostic parameters in geodynamo benchmark

注：urms為速度場(chǎng)的均方根值，Brms為磁場(chǎng)的均方根值，具體定義參考下一章節(jié)；Elsasser數(shù)等效于磁場(chǎng)特征尺度的平方[3]。

除用無(wú)量綱診斷參數(shù)表示數(shù)值解以外，還可以用診斷圖來(lái)判斷數(shù)值模擬的結(jié)果是否達(dá)到穩(wěn)定狀態(tài)，具體的診斷圖將會(huì)在下章討論。

1.3 初始條件

地磁場(chǎng)發(fā)電機(jī)數(shù)值模型的初始條件由Christensen等所定義[20]，速度場(chǎng)初始條件為

u=0.

(3)

溫度場(chǎng)初始條件為

3x4-x6)sin4θcos4φ.

其中x=2r-ri-ro；

(4)

磁場(chǎng)初始條件為：

1.4 邊界條件

邊界條件對(duì)于地磁場(chǎng)發(fā)電機(jī)數(shù)值模型非常重要，控制著地幔和內(nèi)核的差異轉(zhuǎn)動(dòng)。這種差異轉(zhuǎn)動(dòng)是由凈扭矩產(chǎn)生的，扭矩主要包括黏性扭矩和電磁扭矩。由于所關(guān)注的問(wèn)題有差異，地磁場(chǎng)發(fā)電機(jī)數(shù)值模型的邊界條件分為很多種類。同時(shí)，不同的邊界條件也會(huì)造成定義不同的無(wú)量綱參數(shù)。這一節(jié)系統(tǒng)討論不同的邊界條件對(duì)于發(fā)電機(jī)數(shù)值模擬的影響。

1.4.1 速度邊界條件

對(duì)于速度場(chǎng)來(lái)說(shuō)，最簡(jiǎn)單的發(fā)電機(jī)模型假設(shè)外核為一個(gè)剛性的、不可穿透且共旋的球殼。這意味著在旋轉(zhuǎn)坐標(biāo)下，速度場(chǎng)所有分量在邊界處均為0，即no-slip邊界條件：

[u]=0.

(8)

[ ]表示通過(guò)邊界的區(qū)別。Kuang和Bloxham[8]認(rèn)為如果使用no-slip邊界條件，模型中所使用的過(guò)大的Ekman數(shù)會(huì)導(dǎo)致過(guò)大的Ekman-layer影響，所以他們使用free-slip速度邊界條件從而消除Ekman-layer的影響：

1nu=1n×(σν·1n)=0.

(9)

式中：1n為邊界法向量，σν為黏性應(yīng)力張量。

現(xiàn)有數(shù)值模型不能解釋為什么內(nèi)核會(huì)和地幔共旋。部分模型的坐標(biāo)系下內(nèi)核和地幔的旋轉(zhuǎn)是有差異的[4-5,8]。而轉(zhuǎn)速的差異是由凈扭矩產(chǎn)生的，當(dāng)忽略重力扭矩時(shí)，Glatzmaier和Roberts[29]發(fā)現(xiàn)通過(guò)內(nèi)核邊界的耦合，內(nèi)核每年多旋轉(zhuǎn)幾度。隨后地震學(xué)證實(shí)存在這樣的差異旋轉(zhuǎn)，并證明差異旋轉(zhuǎn)速率很小[30-31]。當(dāng)考慮內(nèi)核和地幔的重力耦合后的發(fā)電機(jī)模型，也得出了較小的內(nèi)核旋轉(zhuǎn)速率[32]。

1.4.2 磁場(chǎng)邊界條件

對(duì)于磁場(chǎng)來(lái)說(shuō)，當(dāng)假設(shè)地幔和內(nèi)核均為絕緣體時(shí)，液核磁場(chǎng)在內(nèi)部和外部均為一個(gè)連續(xù)勢(shì)場(chǎng)?；谇蛑C函數(shù)的譜方法很容易滿足這個(gè)條件[33]，而局域法則不能滿足。在一些數(shù)值模型中，運(yùn)用水平磁場(chǎng)分量為0的邊界條件[7,34-36]。在邊界處只有非零徑向場(chǎng)分量，即稱為“準(zhǔn)真空條件”，是不合實(shí)際的。然而其得到的數(shù)值解，與更真實(shí)的絕緣外部區(qū)域條件所得的數(shù)值解性質(zhì)相似[34]。所以如果不關(guān)注外部場(chǎng)結(jié)構(gòu)的研究，則可使用準(zhǔn)真空條件簡(jiǎn)化數(shù)值方法[3]。

假設(shè)內(nèi)核為絕緣體顯然也是不合實(shí)際的，在通常的數(shù)值模型中，假設(shè)內(nèi)核與外核有相同的電導(dǎo)率。這種情況下，必須在內(nèi)核內(nèi)解磁感應(yīng)方程，其速度場(chǎng)只是單純地描述固態(tài)內(nèi)核坐標(biāo)系下的固體旋轉(zhuǎn)速率。在內(nèi)核邊界，應(yīng)用磁場(chǎng)和電場(chǎng)水平分量的連續(xù)條件。Hollerbach和Jones[37]提出的一個(gè)簡(jiǎn)化的平均場(chǎng)發(fā)電機(jī)模型結(jié)果表明，內(nèi)核使用有限電導(dǎo)邊界條件，對(duì)于使得發(fā)電機(jī)產(chǎn)生一個(gè)穩(wěn)定的偶極場(chǎng)起重要作用：

[B]=[1n×J]=[1n×E]=0.

(10)

E為無(wú)量綱電場(chǎng)。在一系列的MHD模型中，包括一些偶極倒轉(zhuǎn)的模型，Wicht[38]發(fā)現(xiàn)當(dāng)其他條件一致時(shí)，內(nèi)核絕緣和可導(dǎo)結(jié)果相差無(wú)幾。作為對(duì)比，Dharmaraj和Stanley[39]指出，內(nèi)核電導(dǎo)率可以強(qiáng)化模型的偶極占優(yōu)并且減少倒轉(zhuǎn)次數(shù)。他們認(rèn)為Wicht的研究結(jié)果的差異主要是其使用固定熱流邊界條件代替了固定溫度邊界條件。

另外一些模型[5,24]沒(méi)有假設(shè)地幔是可導(dǎo)的，而是用一個(gè)電導(dǎo)率適中的薄層作為地幔的底部。Kuang 和 Bloxham[24]的模型在CMB外面包含一個(gè)電導(dǎo)率為其1/10的D”層。其主要目的是使地幔和外核通過(guò)地磁扭矩達(dá)到核幔耦合。目前，電導(dǎo)層對(duì)于發(fā)電機(jī)過(guò)程的影響并沒(méi)有系統(tǒng)的研究，但是被認(rèn)為影響很小[3]。

1.4.3 熱邊界條件

對(duì)于地球外核來(lái)說(shuō)，驅(qū)動(dòng)其對(duì)流的是熱浮力和組份浮力兩部分，前者的來(lái)源是內(nèi)核固化和地核長(zhǎng)期冷卻所釋放的潛熱，后者則是由于輕元素的釋放從而加固了內(nèi)核[40-41]。實(shí)驗(yàn)證明地核內(nèi)可能存在放射性元素鉀，但是仍有爭(zhēng)議[42-43]。而對(duì)流的地?？刂浦睾说臒釗p耗，其有效地給地磁場(chǎng)發(fā)電機(jī)在CMB強(qiáng)加了一個(gè)固定熱流的邊界條件

(11)

熱流是隨時(shí)空變化的，地幔對(duì)流的時(shí)間尺度是幾十到上百萬(wàn)年，相對(duì)于發(fā)電機(jī)的時(shí)間尺度要更長(zhǎng)一些。因此在地磁場(chǎng)發(fā)電機(jī)模型中，除Driscoll和Olson的模型[44]外，其余的模型均忽略時(shí)間變化，而在一些模型中，考慮了空間變化[3]。

許多模型在內(nèi)核和外核邊界之間強(qiáng)加一個(gè)固定溫度邊界條件

θ=0.

(12)

這個(gè)簡(jiǎn)化條件沒(méi)有物理根據(jù)，但是卻證明熱邊界條件似乎對(duì)發(fā)電機(jī)過(guò)程沒(méi)有影響。然而Sakuraba和Roberts[45]、 Hori等[46-47]認(rèn)為，熱邊界條件對(duì)于流場(chǎng)尺度和磁場(chǎng)形態(tài)有很大影響。內(nèi)核邊界條件的物理形式不如外核邊界條件那么清晰。內(nèi)核任何一個(gè)可能的平移運(yùn)動(dòng)也許會(huì)導(dǎo)致鐵的半球熔化，從而形成穩(wěn)固的分層結(jié)構(gòu)[48]。這類復(fù)雜的研究工作目前還沒(méi)有開(kāi)展。

1.5 數(shù)值方法

目前發(fā)電機(jī)模型的數(shù)值方法有兩類：譜方法和局域法。大多數(shù)發(fā)電機(jī)模型采用的是譜方法。Bullard和Gellman[49]首次應(yīng)用譜方法解發(fā)電機(jī)方程組的微分方程。最初使用的Galerkin方法[49-51]已經(jīng)被更有效的譜方法所替代。對(duì)于地磁場(chǎng)發(fā)電機(jī)來(lái)說(shuō)，首先由Glatzmaier和 Roberts[5]提出(其最早應(yīng)用于恒星發(fā)電機(jī)模型)，隨后被其他模型使用。在譜方法中，未知數(shù)在徑向上采用Chebyshev展開(kāi)，球面上用球諧函數(shù)進(jìn)行展開(kāi)。這使得譜方法求導(dǎo)簡(jiǎn)單，并且有較好的收斂性，但是其在并行計(jì)算中需要消耗大量的節(jié)點(diǎn)進(jìn)行計(jì)算和通信。一些發(fā)電機(jī)模型采用超擴(kuò)散從而保證計(jì)算的穩(wěn)定性[4-5,8,23-24,52-53]。本文將在下一個(gè)章節(jié)具體介紹超擴(kuò)散的引用。

近些年一些發(fā)電機(jī)數(shù)值模型發(fā)展使用了局域法：有限差分法、有限元法、有限體積法和譜元法。選擇局域法有以下幾個(gè)原因：1)可以對(duì)球殼進(jìn)行更靈活的分解；2)一些局域法更容易實(shí)現(xiàn)，而且所有局域法的并行化都比譜方法的并行化簡(jiǎn)單；3)局域法避免了在每個(gè)時(shí)間步的Legendre轉(zhuǎn)換，從而節(jié)省大量的計(jì)算資源[54]。比較不同架構(gòu)上不同代碼的性能是一項(xiàng)艱巨的任務(wù)，特別是考慮絕對(duì)運(yùn)行時(shí)間或計(jì)算性能時(shí)。盡管局域法在大型并行計(jì)算上的運(yùn)行效率比譜方法要高，但譜方法的計(jì)算速度仍然保持在數(shù)千個(gè)處理器以上[55]。

1.6 無(wú)量綱參數(shù)的標(biāo)度律

地磁場(chǎng)發(fā)電機(jī)的數(shù)值模擬還有很多問(wèn)題亟待解決，其中一個(gè)重要的障礙就是模型參數(shù)與真實(shí)參數(shù)的差異問(wèn)題[56]。受限于計(jì)算能力，發(fā)電機(jī)模型難以求解快速旋轉(zhuǎn)以及低黏性參數(shù)(相對(duì)應(yīng)較小的Ekman數(shù))，目前所取參數(shù)(E～10-5)與估計(jì)參數(shù)(E～10-15)甚至達(dá)到近10個(gè)量級(jí)的差異[57]。這一鴻溝依靠計(jì)算技術(shù)的發(fā)展短期內(nèi)仍難以解決。Davies等[58]估計(jì)即便模型所用黏性與估計(jì)黏性僅相差6個(gè)量級(jí)(E～10-9)，也需要54 000個(gè)可用的進(jìn)程運(yùn)行至少13 000天以得到一個(gè)磁擴(kuò)散時(shí)間的數(shù)值，這對(duì)于數(shù)值計(jì)算來(lái)說(shuō)是個(gè)巨大的挑戰(zhàn)。

鑒于數(shù)值計(jì)算的參數(shù)取值與實(shí)際值的巨大差距，標(biāo)度律的研究對(duì)于發(fā)電機(jī)研究十分重要。一方面，完善的標(biāo)度律是理解發(fā)電機(jī)理論的重要組成部分。另一方面，其可以預(yù)測(cè)地球不同時(shí)期的古磁場(chǎng)強(qiáng)度，或者將地球發(fā)電機(jī)的特性與其他行星上的發(fā)電機(jī)相關(guān)聯(lián)[59]。從Christensen、Aubert和Tilgner開(kāi)始，研究人員利用發(fā)電機(jī)數(shù)值模擬來(lái)推導(dǎo)(或測(cè)試)標(biāo)度律[56,60-61]。其中有關(guān)于某一個(gè)參數(shù)的標(biāo)度律，如Ra[28,62-64]，E[65]和Pm[66-67]；也有一些是系統(tǒng)性研究各個(gè)無(wú)量綱參數(shù)的標(biāo)度律[9,14,56]。研究得到許多比例關(guān)系及定量結(jié)論，例如發(fā)電機(jī)磁場(chǎng)在什么參數(shù)區(qū)間表現(xiàn)為偶極場(chǎng)占優(yōu)，何時(shí)由多極場(chǎng)主導(dǎo)[4]；如何由給定參數(shù)估算流場(chǎng)和磁場(chǎng)的強(qiáng)弱[68]，以及磁場(chǎng)強(qiáng)弱是否與各種擴(kuò)散系數(shù)相關(guān)[56]等等，具體可參考Christensen[59]的研究。

2 MoSST模型介紹

MoSST模型由Kuang和Bloxham[8]提出，其定義的無(wú)量綱參數(shù)為：表征熱驅(qū)動(dòng)大小的Rayleigh數(shù)Ra，表征黏滯效應(yīng)與旋轉(zhuǎn)效應(yīng)之比的Ekman數(shù)E，表征慣性項(xiàng)和科里奧利力項(xiàng)之比的磁Rossby數(shù)Ro，和表征熱擴(kuò)散與磁擴(kuò)散之比的磁Prandtl數(shù)qκ。MoSST模型在半徑方向采用4階緊致有限差分法求解，使用Chebyschev多項(xiàng)式展開(kāi)作為徑向網(wǎng)格配置節(jié)點(diǎn)。對(duì)于內(nèi)核、外核和D”層，選取的節(jié)點(diǎn)數(shù)分別為35，39和19。在球面上，采用球諧展開(kāi)和擬譜法求解，用FFT進(jìn)行快速譜變換。球諧展開(kāi)的截?cái)嚯A數(shù)分別為：lmax=33,mmax=21。根據(jù)截?cái)嚯A數(shù)，設(shè)置球向網(wǎng)格節(jié)點(diǎn)數(shù)分別為：θ向50個(gè)，φ向64個(gè)。模擬中所求解的網(wǎng)格數(shù)一共為277 830個(gè)。

為避免高階球諧截?cái)鄬?duì)系統(tǒng)的影響，以及減少?gòu)某跏紶顟B(tài)到最終狀態(tài)之間的過(guò)渡時(shí)間，MoSST模型中引入超擴(kuò)散[24]，定義如下:

(13)

式中l(wèi)為球諧階數(shù)。同樣的，引入熱超擴(kuò)散和磁超擴(kuò)散。該處理方法由Glatzmaier和Roberts[5]得到首個(gè)三維自洽發(fā)電機(jī)數(shù)值模型(GR95)時(shí)引入。但以上形式相對(duì)于GR95，已經(jīng)是一種弱形式。本文模擬中設(shè)定l0=4,ε=0.05。

2.1 診斷輸出

由于發(fā)電機(jī)方程組是一個(gè)各物理場(chǎng)高度耦合的非線性系統(tǒng)，當(dāng)初始狀態(tài)輸入以后，在各無(wú)量綱參數(shù)的約束下，系統(tǒng)需要一段時(shí)間進(jìn)行演化迭代以達(dá)到與所選參數(shù)相對(duì)應(yīng)的最終狀態(tài)。判斷系統(tǒng)是否已經(jīng)演化到終態(tài)，依賴診斷輸出而實(shí)現(xiàn)。圖1給出一個(gè)診斷輸出圖的例子，橫坐標(biāo)為演化時(shí)間，以磁自由擴(kuò)散時(shí)間τf=r02/π2η為度量尺度(取r0=3.5×106m,η=2 m2/s,τf≈2×104a)；圖1的縱坐標(biāo)分別表示速度場(chǎng)(VF)、磁場(chǎng)(MF)和溫度場(chǎng)擾動(dòng)(TP)、以及3個(gè)物理場(chǎng)所對(duì)應(yīng)的典型空間尺度:速度場(chǎng)典型尺度(VFLS)、磁場(chǎng)典型尺度(MFLS)、溫度擾動(dòng)典型尺度(TPLS)。

圖1 診斷輸出示例Fig.1 Diagnostic outputs example

診斷輸出中的物理量為整個(gè)外核內(nèi)的均方根結(jié)果，以外核流速為例：

(14)

其典型尺度也是均方根結(jié)果，由以下公式計(jì)算得到：

(15)

結(jié)合圖1，判斷穩(wěn)定輸出的根據(jù)是：所有物理量及其典型空間尺度都已穩(wěn)定達(dá)到一個(gè)τf以上。最終取的輸出結(jié)果為一個(gè)穩(wěn)定的τf以上時(shí)間范圍內(nèi)(不包含開(kāi)始的過(guò)渡階段)的平均結(jié)果，由圖中直線所示。本文所采用的MoSST模型是串行計(jì)算，計(jì)算時(shí)間依參數(shù)選擇、初始狀態(tài)、超擴(kuò)散大小以及硬件性能而定。一般而言，得到圖1所示的穩(wěn)定輸出結(jié)果需要4～5天時(shí)間，但在極端情形下，得到一組穩(wěn)定輸出需要兩周甚至更長(zhǎng)時(shí)間[69]。

2.2 流場(chǎng)和磁場(chǎng)形態(tài)

圖2為地磁場(chǎng)發(fā)電機(jī)流場(chǎng)形態(tài)圖，圖中給出的是子午面剖視圖，其包含左右兩個(gè)部分，其中左半邊代表軸對(duì)稱(m=0)環(huán)型流場(chǎng)，右半邊的線條代表軸對(duì)稱極型流場(chǎng)。每幅圖都是一段時(shí)間內(nèi)的輸出平均后的結(jié)果。環(huán)型流場(chǎng)(uT)和極型流場(chǎng)(uP)由下式定義：

(16)

剖視圖的左半邊是速度的φ向分量uφ，右半邊是流線圖(Ψu)，定義如下：

(17)

式中s=r·sinθ，為場(chǎng)點(diǎn)到地球自轉(zhuǎn)軸的距離。

圖2 軸對(duì)稱流場(chǎng)形態(tài)圖Fig.2 Axisymmetric velocity field

由圖2可知，不論是極型場(chǎng)還是環(huán)型場(chǎng)，地磁場(chǎng)發(fā)電機(jī)的流動(dòng)主要集中在正切圓柱(軸向與自轉(zhuǎn)軸重合，與內(nèi)核在ICB赤道面內(nèi)切，與CMB在近極區(qū)外接的柱體)內(nèi)，并且流動(dòng)尺度大小不一。從左側(cè)環(huán)型場(chǎng)來(lái)看，流動(dòng)除了活躍在正切圓柱內(nèi)，在圓柱體的上下兩個(gè)底面流動(dòng)也很顯著。由于MoSST模型中超擴(kuò)散的使用，壓制了小尺度的對(duì)流，使得數(shù)值結(jié)果為尺度較大的運(yùn)動(dòng)，真實(shí)地核內(nèi)的流動(dòng)是非常復(fù)雜的。

同樣圖3給出軸對(duì)稱磁場(chǎng)的分布圖。從左側(cè)環(huán)型場(chǎng)可以看出，磁場(chǎng)主要集中在正切圓柱內(nèi)，并且以赤道為對(duì)稱軸，上下半球的磁場(chǎng)極性是相反的。從極型磁場(chǎng)可以看出，雖然外核內(nèi)存在多個(gè)磁極子，但是CMB以外的磁場(chǎng)(只有極型場(chǎng)存在)是偶極占優(yōu)[69]。

圖3 軸對(duì)稱磁場(chǎng)形態(tài)圖Fig.3 Axisymmetric magnetic field

3 總結(jié)展望

雖然在可預(yù)見(jiàn)的未來(lái)，地磁場(chǎng)發(fā)電機(jī)數(shù)值模型的參數(shù)達(dá)到類地值是不可能的，但是發(fā)電機(jī)數(shù)值模型依然是研究地磁場(chǎng)的主要手段。經(jīng)過(guò)20多年的研究，地磁場(chǎng)發(fā)電機(jī)模型取得了重大的進(jìn)展，成功模擬出強(qiáng)度相當(dāng)且偶極占優(yōu)的磁場(chǎng)，重現(xiàn)了地磁場(chǎng)倒轉(zhuǎn)、西漂等地磁場(chǎng)行為。

然而我們也要看到局限的一面，基本控制參數(shù)的值不能與地球?qū)嶋H值相匹配。最早的地磁場(chǎng)發(fā)電機(jī)模型中，參數(shù)值被隨意選擇或根據(jù)數(shù)值計(jì)算的可行性進(jìn)行選擇。近年來(lái)，標(biāo)度律的建立取得了重大進(jìn)展，這表明盡管存在較大的參數(shù)差異，但數(shù)值模型可能會(huì)與實(shí)際的地磁場(chǎng)發(fā)電機(jī)處于相同的動(dòng)力機(jī)制。到目前為止，這些標(biāo)度律是經(jīng)驗(yàn)主義，其有效性的范圍還不清楚。為了更好地鞏固它們，運(yùn)用實(shí)際參數(shù)進(jìn)行模擬，則需要發(fā)展數(shù)值計(jì)算，并且運(yùn)用大量的計(jì)算資源。下一代的地磁場(chǎng)發(fā)電機(jī)模擬系統(tǒng)將在千萬(wàn)億次計(jì)算機(jī)上執(zhí)行更真實(shí)的參數(shù)。為了實(shí)現(xiàn)千萬(wàn)億次計(jì)算機(jī)的能力，2D或者3DMPI并行化或MPI-OpenMP 混合并行化是需要的[14]。

地磁場(chǎng)發(fā)電機(jī)的數(shù)值模擬仍有很多問(wèn)題懸而未決：Boussinesq近似合理還是非Boussinesq近似合理，并且會(huì)產(chǎn)生怎樣的影響？核幔邊界和內(nèi)核邊界最合適的邊界條件和耦合條件是什么，以及這對(duì)發(fā)電機(jī)有何影響?組份浮力和熱浮力(雙擴(kuò)散對(duì)流)的相互作用如何影響發(fā)電機(jī)?在地核的頂部是否有穩(wěn)定的分層以及這將如何影響發(fā)電機(jī)？黏性摩擦在發(fā)電機(jī)模型中有多重要，而Ekman數(shù)的影響是否可以忽略不計(jì)?模型計(jì)算方面的問(wèn)題：在徑向方向上分辨率不斷增大的解的收斂速度是未知的；包含湍流運(yùn)動(dòng)的模型的程序精度是否需要新的基準(zhǔn)測(cè)試？修正了的時(shí)間步是否精確有效？

地球并不是太陽(yáng)系中唯一有磁場(chǎng)的行星。研究其他行星的磁場(chǎng)發(fā)電機(jī)將會(huì)極大地促進(jìn)地磁場(chǎng)發(fā)電機(jī)模型的發(fā)展。目前，涵蓋行星磁場(chǎng)的多樣性的比較發(fā)電機(jī)理論還處于起步階段[11]。同時(shí)，地磁數(shù)據(jù)同化有助于改進(jìn)未來(lái)幾十年磁場(chǎng)演化的預(yù)測(cè)。實(shí)際的地磁觀察可以在數(shù)據(jù)同化的環(huán)境中來(lái)構(gòu)建模型，更好地表示地球深部的動(dòng)力狀態(tài)。隨著計(jì)算能力的不斷提升，不久的將來(lái)，地磁場(chǎng)發(fā)電機(jī)的數(shù)值模擬還會(huì)取得更大的進(jìn)展，從而為理解地磁場(chǎng)提供有力的幫助。

感謝焦立果副研究員提供的MoSST模型程序及悉心指導(dǎo)。感謝天河一號(hào)提供計(jì)算資源。