在 COMSOL中模擬冰流

有關調查指出，地球的冰凍圈是“氣候變化的自然整合者，提供了一些最明顯的氣候變化特徵。”（參考文獻1）冰凍圈是氣候系統的一部分，包含冷凍水，占淡水的 80％。使用 COMSOL Multiphysics® 軟體，我們可以類比經典冰流動，分析冰凍圈動態，評估氣候變化的影響，比如海平面上升。

高山冰川的壽命

1773 年，瑞士博物學家 André Bordier 首次用“流體”一詞來描述高山冰川的運動。然而，科研人員花了一個多世紀才就冰川動態的統一描述達成一致。

冰川最令人困惑的一個方面是觀察到冰表現出粘性和塑性，具體取決於冰川。英國物理學家 John Glen 利用應力和應變之間的非線性關係觀察並描述了這種中間特性。這種經典特性稱為剪切稀化，適用於許多不同的流體（如番茄醬和血液）。

任何高山冰川的壽命都可以示意性地描述如下：

1.    雪堆積在氣溫較低的高海拔地區，並壓縮成冰

2.    冰在自重下開始變形並沿斜坡流下

3.    冰在氣溫較高的低海拔地區融化

因此，即使在穩定狀態下（降雪剛好抵消了融化），冰也有一個動態過程：蠕變。這個流體模型是一個標準的納維-斯托克斯方程，它有一個簡化形式：斯托克斯（低雷諾數）近似，其中忽略了平流項。雷諾數的典型值是  ，因此這一假設無疑是成立的。

葛籣流動定律

黏性流動的模擬通常假設應力和應變之間呈線性關係，這個假設描述了牛頓流體。實際上，許多流體在標準條件下是牛頓流體（例如，水和空氣）。然而，許多流體在剪切時表現出粘度的變化。一種較為普遍的方法是使用本構定律來描述粘度  隨剪切速率的某個冪變化的情況。從數學上講，  ,其中  是剪切速率，經典定義為應變率張量  的模。

要完全定義流動定律，需要評估兩個參數：

1.    一致性， 

2.    應力指數， 

對於冰，通常取  。然而，冰的粘度不僅取決於剪切速率，還取決於溫度和壓力。於是我們定義一致性來表示這些依賴關係。定義冰模型一致性的經典方法是使用阿累尼烏斯定律（參考文獻 1）：  ，其中 是理想氣體常數，  是相對於壓力熔點的溫度。

實際上，壓力依賴性反映在冰的熔點隨壓力的變化上（隨著壓力的增加而降低）。利用克勞修斯－克拉佩龍關係，我們得到  ，其中  是克拉佩龍常數。  和 的值目前存在爭議。（參考文獻 2）

黏性滑動

這種簡明的流動定律是多年來通過嚴格的實驗室工作憑經驗得出的，它無法預測在現實冰川上觀測到的高速度。人們花了很多年才明白原因到底在哪裡。20 世紀 50 年代後期，J. Weertman 提出了基底滑動的概念，並在業內形成共識。

在方程層面，基於分子相互作用的考慮，冰川的基底滑動定律與 H. Navier 一個世紀前提出的黏性滑動 概念沒有什麼不同。然而，對於冰流動而言，這個定律背後的物理過程仍然是一個爭論性問題，並不是這篇文章的主題。我們回想一下，它寫為  ，其中  是基底速度；  是粘度（這裡是非線性的）；  是基底牽引力，或基岩的剪切應力；  被稱為滑移長度（參考文獻 2 ），這個量在冰川流動建模中至關重要，這是因為它占下游品質通量的很大一部分，積少成多就產生流動運動，並且表現出類似於剛性運動的特性。

冰流動建模：現實案例

Mer de Glace，翻譯為“冰海”，是一座高山冰川，位於夏蒙尼山谷上方法國阿爾卑斯山的勃朗峰地塊。它被認為是法國最大的冰川，因其山谷冰川具有相當大的移動速度（每年約 100 米），而且在過去 80 年裡它的體積顯著減小，因此被廣泛觀察和監測。研究表明，該冰川平均每年損失 5 米厚，30 米長。

下圖（右）是山谷冰川的幾何圖形，是使用 COMSOL Multiphysics 和 CAD 幾何內核（可從 CAD 導入模組獲得）建立的。幾何圖形大致模擬了冰海冰川（左）的測量和視覺效果。

左：1909 年冰海冰川的航拍照片。圖片來自公共領域，通過 Wikimedia Commons 分享。（作者添加的注釋。）右：模型幾何結構，按厚度著色。

我們來模擬冰塊在自重和基底滑動作用下沿坡流下的非等溫流動。

就流體而言，流入和流出邊界條件是正常約束，對應於冰的外加壓力，冰不包含在域中，它僅僅對應於指定的流體靜壓（或冰靜壓）。上游邊界壓在域上，因此有助於加快流向速度，而下游邊界阻礙流動。冰川表面是自由表面。

就傳熱而言，表面被認為處於環境溫度。與基岩接觸的邊界通常受到地熱熱通量的影響，地熱熱通量可以被模擬為邊界條件。然而，由於這樣的值在空間上是變化的並且通常是未知的，因此在本例中施加了溫度。這樣，我們確保冰保持在 0°C 以下，從而避免相變和潛熱通量的貢獻。值得注意的是，可以使用具有相變介面的材料來分析這個方面。熱量可以在流入邊界和流出邊界進入和離開域。

使用與幾何結構的縱橫比一致的拉伸網格。

外部天氣條件是地球物理模擬的重要輸入資料。通過“傳熱”介面直接訪問 ASHRAE 2017 資料庫，我們可以在全球超過 6000 個氣象站輸入一年中給定時間的平均外部溫度和風速。這裡，我們使用的資料來自瑞士阿爾卑斯山的聖伯納德 站，該站距離冰海冰川 16 公里，位於 2 月 1 日中午的同一高度。環境溫度施加在冰川表面，風速用於模擬表面的對流熱通量。

模擬結果和討論

首先，我們在沒有基底滑動的情況下運行模擬，看看黏性流動對觀測到的冰川速度有多大影響。預計結果是冰川頂部每年約 120 米，冰川末端每年約 90 米。

正如我們在圖的左側所看到的，僅僅根據此處描述的黏性定律，我們只能得到預期速度的 50％。

我們可以引入滑移長度  的黏性滑移，然後再次運行模擬。下面，我們繪製了兩種情況下沿著冰川中心流線表面的速度。

現在，全球範圍內冰川消失速度要高得多，也更符合預期的幅度。值得注意的是，黏性滑動不僅引入了速度的純粹變化。實際上，作為非線性粘度的函數，它增加了非線性貢獻，因此它不是純粹的剛性。對於滑移長度這個值，滑動貢獻了表面速度的大約 60％。

接下來，我們繼續討論溫度對冰流動的影響，這是氣候變化研究中的重要耦合問題。為了量化全球變暖對冰川的影響，我們設想如下一個實驗。資料顯示，1940 年至 1970 年間，全球氣溫一直處於穩定狀態，因此我們可以假設冰川在此期間達到穩定狀態。測量結果表明，全球平均氣溫在 50 年間上升了約 1 度。因此，我們可以模擬這 50 年間平均氣溫穩定上升 1 度情況下冰的瞬態流動。

為了瞭解這種溫度變化的影響，我們可以繪製表面和下游邊界總品質流出量（單位為兆噸/年）的演變圖。

觀察氣溫開始上升和冰川反應之間的延遲非常有趣。線性升溫始於 1970 年，線性升溫的品質通量始於 1985 年左右。在 1970 年至 1985 年間，品質損失的加速度緩慢增加，直至達到每年約 13 千噸的穩定值。這種延遲主要是由於地表溫度變化傳播到整個冰川的時間（從而增加了整個冰川的平均溫度）。結果，輸出品質通量在此期間增加了大約 10％，導致在此期間淨增加了 15 兆噸的冰損失（與溫度保持穩定的情況相比）。

我們將結果與之前討論的冰海冰川資料進行比較。如果我們所在地區的冰川厚度每年減少 5 米（平均 5500 米長，600 米寬），我們每年將損失 15 兆噸的冰。即使假設所有這些冰通量都會在較低的海拔處融化（事實並非如此），2017 年計算的每年 750 萬噸遠小於過去幾十年冰海冰川的實際品質損失。這是因為模擬沒有考慮負表面品質平衡（積雪減去表面冰的融化）。

就建模而言，表面品質平衡是資料登錄，本身就是複雜物理學的產物。例如，較熱的夏季對冰川有很大的影響，這是因為通常夏季降落的少量雪可以起到防禦太陽輻射的作用，從而保護冰川在夏季不會大量融化。如果夏季沒有降雪，那麼冰川融化量會大得多。這種額外的融化導致液體水通過裂縫大量滲透，最終形成冰下水文網路，主要通過冰-基岩介面的潤滑和水壓“抬升”冰川，在基底滑動中發揮重要作用。

由於幾何形狀會影響動力學，因此我們針對冰川的給定幾何形狀執行模擬，忽略了幾何形狀因表面品質平衡和動力學而發生的演變，這也是很重要的一個方面。

關於冰流動建模的總結性思考

本文介紹了如何用 COMSOL Multiphysics 建立和求解一個簡單的冰川流動模型。COMSOL® 軟體為此類建模中涉及的大多數問題提供了專用功能。本例和一般冰川學的主要局限是資料；通常包括地形資料、基底滑動長度、表面品質平衡、堆積和融化。

參考文獻

1.    Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, T.F. Stocker et al. (eds.), Cambridge University Press, 2013.

2.    R. Greve and H. Blatter, “Dynamics of Ice Sheets and Glaciers,” Springer, 2009.