JFM：CFD模拟地幔对流与板壳运动，新视角解读地球物理板块构造

1. 绪论    

地球内部让几代科学家着迷（Plummer，Carlson & Hammersley 2015）。其中，Leonardo da Vinci（1452-1519）是注意到我们地球不停的地质运动的先驱之一，因为他观察到山脉中存在海洋化石。我们现在知道地球的大陆不会停留在原地，而是经历构造运动，地球地幔中的热对流被认为是这些运动的驱动力（Kious & Tilling 1996）。

当流体的温度不均匀导致密度和浮力不均匀时，就会发生热对流，因此温暖的流体上升，而寒冷的流体下沉。这里流体的定义可以非常广泛，因为现代地质学家证实，即使是地幔在大时间尺度上也像流体一样流动（Turcotte & Schubert 2002）。普朗特数（Pr）定义为地幔运动粘度和热扩散率之比，估计约为10²³（Meyers等 1987）。

地球的核心比表面温暖得多，不稳定的浮力足以驱动地幔对流。作为浮力和粘性效应之间相对强度的衡量标准，地幔中的瑞利数（Ra）约为10⁶（Selley，Cocks & Plimer，2005）。在Rayleigh–bénard对流的充分研究案例中，如此高的Ra已知会导致湍流运动（Ahlers，Grossmann & Lohse 2009）。由于地幔像流体一样对流，其表面流动运送了大陆板块，导致了它们的构造运动。

由于地球的空间尺度大，地幔对流的时间尺度长，板块构造的地球物理研究侧重于大陆的当前状态以及预测其重要后果，如地震（Plummer等 2015年）。另一方面，数值模拟（Howard, Malkus & Whitehead 1970; Whitehead 1972, 2023; Whitehead & Behn 2015; Mao, Zhong & Zhang 2019; Mao 2021）和实验室规模的实验（Elder 1967; Zhang & Libchaber 2000; Zhong & Zhang 2005, 2007a,b; Whitehead, Shea & Behn 2011）已被证明是理解板块构造动力学的有效手段。

Elder（1967）的早期实验展示了如何在实验室中复现构造运动。在实验室中，石蜡流体层被用来模拟地幔，而漂浮在顶部的塑料薄片被用作模型大陆板块。当石蜡从下面被加热时，对流发生，由于下面对流的剪切作用，板块移动。这样一个简单的实验装置也显示了大尺度的运动；正如Zhang & Libchaber（2000）所观察到的那样，板块在流体表面的两个边界壁之间周期性地移动。通过更详细的研究（Zhong & Zhang 2005，2007a,b），浮板的尺寸对板块运动有很大影响，小板块呈现周期性运动，而大板块则被困在流体表面的中间。将移动热源放置在热对流流体顶部也会产生板块运动，Howard等（1970年）和Whitehead（1972）对此进行了实验研究。尽管这些作品中呈现的几何形状、物理参数和时间尺度与地幔对流非常不同，但它们揭示了令人惊讶的动力学，最重要的是，为大陆漂移背后的流体-结构相互作用机制提供了宝贵的见解。

板块构造的数值探索在过去几十年中发展迅速。Gurnis（1988年）提供了第一个与时间相关的大陆漂移数值模拟，其中允许多个大陆合并和分离。在这项工作中，许多其他数值和理论的努力（Zhong & Gurnis 1993; Lowman & Jarvis 1993, 1995, 1999; Lowman & Gable 1999; Zhong et al. 2000; Lowman, King & Gable 2001）采用地幔的地球物理参数，使我们对地球内部的了解迅速推进。同样清楚的是，大陆板块和地幔对流之间的双向耦合导致了不同的动力学（Gurnis 1988; Zhong & Gurnis 1993; Phillips & Bunge 2005; Whitehead & Behn 2015）。最值得注意的是，观察到大板块具有更一致的运动，而小板块则倾向于偶尔运动（Gurnis 1988; Whitehead & Behn 2015）。这些观察结果最近由Mao等（2019年）和Mao（2021）通过解析数值模拟进行了检查。

上述工作证实，大陆板块不仅对下方的地幔流平流是被动的，而且还通过热毯效应影响流动结构:众所周知，与大洋地壳相比，大陆地壳的热通量要低得多，因为其地壳深度较大（Mao 2021），因此大陆板块本质上是一个毯子，可以防止热量逃逸并使下方的地幔升温。温暖而轻盈的地幔往往会上升，形成向上的对流气流。随着这种流动向地球表面移动，它会分散并在大陆板块下方产生流体作用力，将板块带走。这是目前对大陆漂移的理解，热毯效应已经在数值和实验上得到验证（Gurnis 1988; Zhong & Zhang 2005, 2007a,b）。

本文旨在通过低维模型为板块构造建模提供一个新的角度。在对二维周期域中的板块流相互作用进行直接数值模拟（DNS）后，我们提出了一个板块运动的随机模型，并展示了移动板块如何与其下方的对流流体流进行机械和热耦合。该模型仅保留了热对流的最基本物理性质，恢复了在完整DNS中观察到的动力学，并捕捉了Gurnis（1988）、Whitehead & Behn（2015）、Mao等（2019）和Mao（2021）看到的板块动力学转变。

在下文中，我们在第2章总结了方程和数值方法，并在第3章中给出了数值结果。在第4章中系统地推导了随机模型，并在第5章中讨论了它在具有各种纵横比的对流域中的应用。最后，我们在第6章中总结并讨论了数值结果。

2. 数值模型    

2.1流动方程

我们的数值模拟配置如图1所示，其中以x=x_p位置为中心的固体板漂浮在对流流体的顶部，该流体在y方向上受到限制，在x方向上具有周期性。在整个研究中，所有长度都由流体深度H重新调整，时间由扩散时间H²/κ（其中κ是热扩散率）重新调整，温度由底部和顶部自由表面之间的温差T重新调整。流体域的x方向是周期性的，周期Γ=D/H（其中D是域宽度），因此整个计算域为x∈(0, Γ)和y∈(0, 1)，如图1所示。使用Boussinesq近似，流速u=(u, v)、压力p和温度θ∈[0, 1]的所得偏微分方程为：

这里，瑞利数为Ra=αgΔTH³/νκ，普朗特数为Pr=ν/κ，其中ν、α和g分别为运动粘度、流体热膨胀系数和重力加速度。可以针对地球物理地幔对流对流动解算器进行简单修改，但由于我们希望考虑更普遍的流体-结构相互作用情况，并将我们的理论应用于未来的实验室实验，因此我们在本研究中保留了流体和固体板的惯性。

图1：移动板块和对流流体示意图

流体域以y∈(0, 1)为界，并且在x上是周期性的，宽度为d的浮动板具有中心位置x_p

2.2边界条件

在没有板的情况下，边界条件很简单：流速u=(u, v)在流体/固体边界处无滑移，在空气/流体界面处无剪切。温度θ在底部为1，在空气/流体边界为0。

这产生了底面y=0的边界条件:

在顶部表面，板有效地屏蔽了热量的逃逸，因此我们在那里取θ_n=0，同时将流动设置为相对于移动板无滑动。由此产生的边界条件为：

或者，这些条件可以强制执行为

其中1P是一个指示函数，在板下取值1，否则取值0。

2.3板块动力学

流体剪切力直接驱动板块运动，所以

这里，u_p=x_p^(导数)）是板块速度，m=ρd是具有线性密度ρ和宽度d的板块的无量纲质量，积分面积P={x|x∈(x_p-d/2, x_p+d/2)}是板块下的区域。

2.4.参数和数值方法

附录A中详细说明了用（2.4）、（2.7）和（2.8）求解（2.1）-（2.3）的数值方法。在附录A中，我们使用傅立叶-切比雪夫谱方法来获得精确的数值解。在所有模拟中，我们选择Ra=10⁶，Pr=7.9，Γ=1–16，m=4d（其中d为板宽），与实验中的对流参数相匹配（(Zhang & Libchaber 2000; Zhong & Zhang 2005）。

图2：漂浮在对流流体顶部的板块运动。（a）各种大小板块的轨迹。覆盖比是Cr=d/Γ，时间t₀标志着板块运动的开始。（b）板块的总位移，其中d_p通过d_p=|x_p^(导数)|定义。（c）当Cr>0.33时，平均板速度v_p=<|x_p^(导数)|>变高，在Cr=0.56时达到最大值。在这些模拟中，Γ=4，ρ=4，Ra=10⁶，Pr=7.9。

3. 数值结果    

图2（a）显示了不同大小的板块的几个轨迹，我们可以立即看到板块的大小如何影响其运动。我们定义了覆盖比Cr=d/Γ，以测量宽度为d的板覆盖了多少自由表面，其中对于本节中的所有结果，流体长宽比为Γ=4。对于小板来说，它们的净位移很小，从图2（b）所示的总位移可以更好地看出这一点。如图2（a）中绿色轨迹所示，当Cr大于0.33时，增加板尺寸会出现线性运动。这些轨迹会发生逆转，因为湍流作用力会产生有效噪声。随着Cr进一步增加，线性运动变得更加持久，因为当图2（a）中Cr →1时，板块运动的逆转变得罕见。我们注意到，在地球物理斯托克斯流模拟（Mao 2021）中已经观察到类似的动态行为，因此对于不同的流态，移动板块和下方流之间的耦合机制必须相似。

从总位移d_p可以看出，在Cr≈0.5时达到最大板速度，这可以通过绘制图2（c）中的时间平均板速度v_p来确认。对于小板，平均速度v_p仍然较低，但当Cr>0.33时，平均速度v_p显著增加，并在Cr=0.5时达到最大值。

为了研究动态之间的转换，流体中的典型流量和温度分布如图3所示。在图3（a）中，一个Cr=0.125的小板被放置在对流流体上，并被x_m处的下行流体中心吸引（图1），在此处表面流形成一个汇。这种下沉对板块来说是一种稳定的平衡，因为任何偏离这种下沉的情况都会导致作用在板块上的恢复流体力。图3（b）进一步显示了该流汇的结构，其中y平均温度和y平均垂直流速均达到最小值。

按照这种表面流动模式，板块位移x_p是随机的，如图3（c）所示。由于湍流的随机作用力，板块位置受到噪声的影响，如图3（d）所示，噪声会影响板块速度，其直方图呈高斯分布。板块遇到来自流动的强大“风”的情况很少见，但并非不可能，这种风可以推动板块远离流动汇，穿过流动源，然后再次回到汇，从而导致图3（c）所示的x_p跳跃。

图3：（a–d）Cr=0.125和（e–h）Cr=0.5的浮板的动力学

（a）Cr=0.125的板下流量和温度分布的典型快照

（b）对应于（a）的y平均温度和垂直流速，阴影区域表示板块位置

（c）板块位移是一个随机过程

（d）板块速度是均值为0的随机变量，其概率密度函数是高斯分布，如直方图所示

（e）在Cr=0.5时，典型的流量和温度分布显示板由表面流输送

（f）对应于（e）的y平均温度和垂直流速

（g）板块位移在时间上是线性的，表明是单向平移

（h）板块速度具有非零平均值

这些模拟的补充动图可在https://doi.org/10.1017/jfm.2023.1071.获得。

覆盖比为0.125的小板浮在对流流体上的动力学。上部面板显示了流场和温度场以及板的位置。下方面板显示了y轴平均温度和垂直流速。

覆盖比为0.5的大板在对流流体的顶部平移。显示了流场和温度场（上图）及其y平均值（下图）以及板的位置。

图3（e）显示了Cr=0.5的板的运动情况。在这种情况下，板块运动是单向的，如图3（g）和3（h）所示，向上的速度具有非零平均值。如图3（e）和3（f）所示，移动板倾向于位于流汇和流源之间。随着表面流将板块推向其汇点，流动温度的分布也会发生变化，导致移动的板块追逐移动的表面流汇点。

这是热毯效应的直接结果：当板块足够大时，其下方的温度会上升，因为热量无法从那里逃逸。这种局部升温改变了流动温度，并有效地将冷的下行流体推走，导致流动汇位置的移动。总的来说，板朝着冷流散热器移动，同时将散热器推开。因此，对于这个看似复杂的流体-结构相互作用问题，存在一个简单的动力学，我们将从这些观察中推导出一个模型。

4. 随机性模型    

如图3所示，y平均温度θ的变化强烈影响流型。为了捕捉θ的变化和周期性质，我们用其最低的非平凡傅立叶模式来近似它，

其中x_m是图1中表面流汇的位置，r=2πΓ^-1是波数，常数α衡量温度变化的强度。

由该温度分布引起的表面流速U(x, t)=U(x, 1, t)可近似为

其中β>0为表面流动强度。事实上，这种表面流动剖面在x=x_m处有一个汇点：与x_m的小偏差导致x<x_m时U>0，而x>x_m时U<0，因此在局部，流动指向x=x_m。

我们注意到，该表面流动剖面与固体/流体边界处的板块速度不匹配，u和u_p之间的不匹配允许我们估计∂u/∂y(x, 1, t)和由此产生的剪切应力，从而导致板块加速

这里，δ是动量边界层厚度（Schlichting & Gersten 2016），由Ra和Pr确定，因此我们在研究中假设其为常数。我们还包括一个标准差为σ的白噪声，代表湍流作用力。总体而言，（4.2）和（4.3）代表了同时具有大尺度环流和小尺度湍流的流动剖面，与高Ra热对流的观测结果一致（Ahlers等 2009；Huang & Zhang 2022）。

为了将移动的板块模拟为热毯，我们看一下y平均热量方程，

这里，我们忽略了流体平流，q(x, t)=(∂θ/∂y)(x, 1, t)(∂θ/∂y)(x, 0, t）是通过位置x的垂直热通量。假设离开流体-空气界面的热量遵循牛顿冷却定律，并且没有热量穿透板，我们可以将热量方程改写为

当x∈P时指示函数1_P(x）为1，否则为0，常数γ模拟冷却速率。我们现在从（4.1）中插入θ的值并在x上积分，这导致x_m的常微分方程，

（4.3）和（4.6）中的积分可以精确计算。定义相位角φ=r(x_p−x_m)，我们得到(u_p, φ)的闭合动力系统:

其中λ=Pr/(ρδ)。一旦知道(u_p, φ)的动态特性，就可以通过x^(导数)_p=u_p和x_m=x_p-r^-1φ计算出(x_p, x_m)的动态特性。

图4：随机模型的模拟轨迹

（a）小盘的动态是随机游走的

（b）中等大小的板块具有非零的平移速度，其方向可能发生逆转

（c）大板块的平移运动更持久

（d）所有模拟路径的轨迹如图2（a）

该模型中有四个参数：β是表面流动强度，λ=Pr/（ρδ）是阻尼系数，γ是表面冷却速率，σ是随机流体作用力。在物理上，表面冷却速率γ受表面流动强度β的影响，因此我们在该模型中取γ=rβ，这将导致正确的动力学。其余参数可通过数值模拟进行估算。

从φ（0）=0和随机值u_p（0）开始，图4（a）–4（c）显示了不同Cr下x_p（t）的一些典型轨迹。在图4（a）中，小板（Cr=0.2）的轨迹是噪声驱动的。对于Cr=0.5的中板，图4（b）显示其轨迹由反向线性平移组成。对于Cr=0.8的大板，图4（c）表明平移是单向的。图4（c）中的板速度与图2中的完全DNS结果相当，并且该速度随着Cr的进一步增加而降低。不同尺寸板的典型位移x_p如图4（d）所示，与图2（a）相似，在Cr≈0.3时，噪声驱动运动和线性运动之间存在过渡。因此，这个简单的模型包含了整个数值模拟的所有关键特征。

没有噪声，动力系统（4.7）和（4.8）的临界行为可以进一步分析。对于小Cr，很容易看出（4.7）和（4.8）有u_p=0，φ=2πN（其中N是整数）作为平衡，它们是稳定的，反映了板块运动的被动状态。随着Cr的增加，在Cr^*=1/3时出现新的平衡。当Cr>Cr^*时，有可能得到非零板速度u^*_p=（β/π）u^^*_p，其中

平衡φ^*可以由下式确定

这些状态代表板块的平移。我们注意到（4.9）和（4.10）只是Cr的函数，因此与该模型中假设的所有其他参数无关。为了恢复量纲板速度u^*p，只需知道流速系数β/π。

U^{^*}_p和φ^*的可能值如图5所示。对于Cr<Cr^*=1/3，u^*_p=0和φ^*=2πN是唯一可能的平衡，反映了始终被表面流汇吸引的小板的被动性质。当Cr>Cr^*时，新相位表现为φ^*=（2N+1）πarccos（[(2Cr)^-1-1](cosπCr)^-1），这是（4.10）的解，在较大Cr时变得稳定。它们表明较大的板块倾向于位于表面流汇(φ=2Nπ）和源[φ=(2N+1)π]之间，这证实了我们在图3中的观察。当表面流从其源头指向其汇点时（见图5b中的箭头），这些新相位表明了两种可能的板块速度，由（4.9）式给出，如图5（a）所示。此外，图5（a）与图2（c）相似，当Cr→Cr^*和Cr→1时，板块速度为零，并在Cr=0.5左右达到最大值。

图5：板速度和相位角的平衡值:（a）无量纲板速度up对于小Cr为0，但随着Cr增加而变成平移；（b）相位角φ。蓝线代表板块被流汇吸引的被动状态，红/橙色曲线显示平移状态的平衡阶段。表面流向标有箭头。

通过这个简单的模型，我们可以清楚地看到固体板如何与下面的流体相互作用的物理过程，覆盖比Cr可以作为热毯效应强度的衡量标准。对于小Cr，热毯效应较弱，因此板块运动对流体是被动的。当Cr超过Cr^*时，热毯效应会强烈到足以改变流量和温度分布，产生上升流，导致板块运动。图2（c）和图5（a）表明，板块速度在Cr≈0.5时达到峰值，并在Cr→1时消失，这也反映了热毯效应和流动对流之间的竞争。随着Cr的增加，更多的自由表面面积被板覆盖，下方的流体力在更大的范围内得到平均。由于该域可能涵盖上升流和下降流，总流体力受Cr影响。这可以在（4.3）式中看到，其中表面流速U提供了板块加速度。在（4.2）式中，U被模拟为正弦曲线，该曲线在（4.3）式中被积分，因此将（4.3）式中的积分面积增加到开放表面的一半（Cr=1/2）将覆盖U的最高贡献。进一步增加覆盖面积将因此减少U的贡献，因为积分域超过正弦函数的半个周期。在极端情况下，Cr=1表示U在整个周期内的积分，导致流体力为零，如图5（a）所示。

5. 纵横比效应    

到目前为止讨论的所有结果都集中在纵横比Γ=4的对流流体域，这与许多实验研究的几何形状相匹配。随着更复杂的多卷流结构的出现，改变区域纵横比肯定会影响对流流体和移动板块的动态特性（Ahlers等人，2009年）。在本节中，我们研究了Γ=2、4、8和16时的板块动力学，同时保持其他动力学参数与2–4中描述的相同。

图6：对流流体在不同纵横比下的快照:（a）Γ=2，（b）Γ=8，（c）Γ=16。在这些模拟中，所有其他参数都固定在Cr=0.2、ρ=4、Ra=106和Pr=7.9。

图7：长宽比（a）Γ=2，（b）Γ=8和（c）Γ=16的DNS模拟轨迹。在这些仿真中，ρ=4，Ra=106，Pr=7.9，如图2和3所示。（d）分离被动和平移状态的临界Cr*与Γ成2/3幂律比例

图6显示了Γ=2–16时DNS结果的一些典型温度分布。很明显，随着Γ的增加，更多的对流单元出现，出现了更复杂的流动结构。图7显示了不同Γ下的板块轨迹；一个共同的特征再次出现，当大板块平移时，小板块移动很小。

为了将我们的简单模型扩展到高展弦比的情况，我们考虑具有更复杂的空间相关性的整体温度和表面流动剖面。

其中整数k是傅立叶光谱中最主要的波数。包含波数k是受对流可能具有多卷结构这一事实的启发，因此上面的温度和流量分布图表明流体中存在2k个对流卷，具有k个表面汇和k个表面源。我们跟踪一个表面下沉x_m(t)，这也是最低流体温度的位置。

通过类似于第4章中的推导，我们可以再次获得相位φ=r(x_p-x_m）和板块速度u_p的随机动力系统:

在没有噪声的情况下，很容易验证（5.3）和（5.4）具有被动平衡u^*_p=0和φ^*=2mπ/k，其中m是任意整数。这种被动状态的雅可比矩阵是

在Cr→0的极限中求J的迹和行列式

因此J的两个本征值都是负的，表明小Cr的稳定被动态，并证实了我们的数值结果。

当J的一个或两个本征值变为正值时，被动状态的稳定性就会丧失，因此det(J)=0提供了一个确定临界Cr^*的方程。经过简化后，我们有

已知参数β和Γ，（5.8）的根可以用数值确定。在实验室实验和地球物理板块构造中，可以适当估计表面流速尺度β和通气系数Γ，并用于确定Cr^*，因此（5.8）提供了通过物理参数确定板块活动性的可能性。

我们注意到，对于较大的Γ，Cr^*通常较小，如图7所示，在此极限中（5.8）的根可以表示为

这里，我们使用了关系式r~Γ^-1和k~Γ，后者表明对流辊的数量与长宽比Γ成比例。等式（5.9）确实可以得到验证，如图7（d）所示，从DNS数据测得的临界Cr确实遵循2/3乘幂定律。

通过（5.3）和（5.4）式的动力学可以研究更多的内容，未来的实验研究肯定可以用来进一步解决板块和下面的对流流体之间的相互作用。

6. 讨论    

在这项工作中，我们从数值和理论上探讨了运动平板和其下对流流体之间的机械和热耦合。受到现在和过去作品的启发（Gurnis 1988; Whitehead & Behn 2015; Mao et al. 2019; Mao 2021），我们提出了一个随机模型，表明板块大小（覆盖比）是热毯效应强度的决定因素。对于小板块，屏蔽效应较弱，板块运动对流动结构是被动的；对于大板块来说，下方的流动变得足够温暖，并形成一个上升流中心，将板块推开并导致其平移。提议的简单模型由关于流动及其与板块的机械和热耦合的最小假设组成；然而，它能够预测板块运动的动态转变。虽然直接数值模拟（DNS）是在类似于实验室实验的参数范围内进行的（Zhong & Zhang 2005，2007a，b），但这种随机模型与雷诺数无关，因此可应用于实验室和地球物理场景。

实验室规模的实验通常在有界对流系统中进行，因此平板仅限于在两壁之间移动。在有界对流中，流动结构及其与板块的耦合非常不同，因为先前的工作表明，流动具有两个反向旋转的大尺度环流，其强度受板块位置的调节（Zhong & Zhang 2005，2007a，b）。该板块还受到固体边界的限制，因此一旦接触到墙壁，它就会停止移动。在这种情况下，对板块动力学进行建模并非易事，随机变分不等式（Huang et al. 2018）等先进工具可以作为分析这种板块-边界相互作用的数学手段。对流流体的周期域是实验验证我们的随机模型所必需的，这样的实验也将更接近地模拟地幔对流。由环形对流域组成的未来实验可以提供更多细节来验证随机模型。

虽然这里我们只研究了单个板块的动力学，但只要能够正确模拟它们之间的相互作用，我们的数值方法就能够处理多个板块。我们目前正在研究这种相互作用，这导致了更加多样化和不可预测的动态。例如，如果两个小板块的Cr<1/3，但它们的总大小达到Cr>1/3，那么我们已经看到每个单独的板块随机移动，但它们的组合“超级大陆”可以平移。在板块构造的地球物理案例中，板块相互作用是许多火山活动和山脉形成的原因，因此了解附近板块的聚散运动可能会为这些地球物理事件背后的流体力学提供新的见解。

为了简单起见，这项工作中的模拟和模型都是二维的，将我们的结果扩展到三维是当前的优先事项。使用Chebyshev–Fourier–Fourier方法，我们实现了数值求解器，用于演化位于在两个水平方向上具有周期性的三维域顶部的板。此外，通过Chebyshev–Chebyshev–Fourier方法也可以模拟球形壳上的板块构造，这是一种更接近地球物理板块构造的配置。通过分析那里的DNS结果，我们希望进一步开发我们的随机模型，并将其用于解决地球内部发生的流体-结构相互作用。

最后，我们注意到地球物理板块构造比迄今为止进行的任何实验或数值模拟都要复杂得多，因为地球内部的环境如此复杂，现代科学仍在探索之中。但是尽管目前的简单模型不能完全捕捉大陆漂移的动态，我们希望它们仍然可以提供一些关于地球物理学的流体力学见解。

翻译转载自：《Journal of Fluid Mechanics》:"Covering convection with a thermal blanket: numerical simulation and stochastic modelling"