利用穆顿成像技术检测和监测地质构造中的二氧化碳羽流：一项可行性研究

《International Journal of Greenhouse Gas Control》：Investigation of Muography for CO2 plume detection and monitoring in geological formations: A feasibility study

　　哈米德·巴希里（Hamid Basiri）| 松岛淳（Jun Matsushima）| 小玉雅志（Masashi Kodama）| 横田敏之（Toshiyuki Yokota）| 穆罕默德·Y·阿里（Mohammed Y. Ali）| 法特赫·布查阿拉（Fateh Bouchaala）| 艾哈迈德·埃莱斯拉姆布卢利（Ahmed Eleslambouly）| 田中宏之辅（Hiroyuki K.M. Tanaka）| 埃杜阿多·科蒂纳·吉尔（Eduardo Cortina Gil）| 安德烈亚·贾曼科（Andrea Giammanco）
东京大学前沿科学研究生院环境系统系，日本东京

**摘要**
碳捕获与封存（CCS）通过将二氧化碳注入深层储层来减少温室气体排放，但追踪二氧化碳羽流范围和孔隙空间饱和度仍具有挑战性。基于斯莱普纳（Sleipner）CCS项目（北海）的理想化可行性研究，我们评估了宇宙射线μ成像技术在受控模拟条件下检测二氧化碳羽流层和区分不同饱和度水平的物理敏感性。我们提出了一种两阶段蒙特卡洛工作流程：MUSIC（μ子模拟代码）根据实验推导的海平面μ子谱生成并传输μ子，使其穿过上覆岩层；PHITS（粒子与重离子传输代码系统）随后在乌蒂拉沙（Utsira Sand）的九层羽流模型中传播存活的μ子。与单阶段全传输模拟相比，该工作流程将运行时间减少了近三个数量级，使得千米级μ子生成和传输模拟成为可能。我们模拟了25%、50%和100%的层平均二氧化碳饱和度，并发现μ子传输随饱和度升高而单调增加，而路径平均密度图能够恢复多层几何结构。在理想化假设下（仅垂直入射μ子、800米半径的虚拟探测器具有完美效率、均匀饱和度），100%和50%的饱和度情况下，在穿过羽流模型的10^7个μ子内即可实现≥3σ的统计分离；区分50%和25%的饱和度需要更多的μ子样本。真实探测器的采集时间与几何接受度和效率成反比。这些结果表明，在理想化条件下，μ成像技术对孔隙流体替换引起的整体密度变化具有物理敏感性，并为未来的实验设计和与时间序列地震监测的比较建立了计算框架。

**1. 引言**
碳捕获与封存（CCS）被广泛认为是减少温室气体排放的有效措施。在CCS过程中，二氧化碳在主要点源被捕获并作为超临界流体注入深层地下储层（Bui等人，2018）。确保这些储层的长期完整性需要频繁且准确地评估其中的二氧化碳含量。传统的地震调查可以监测和界定二氧化碳羽流，但由于振幅-饱和度关系高度非线性，量化饱和度仍具挑战性。地震方法对低自由气体比例（百分之几的量级）非常敏感，但在中等至高二氧化碳饱和度下则逐渐失去诊断能力，导致反演结果不唯一（Landr?等人，2001；Grude等人，2013；Landr?和Zumberge，2017；Knight等人，1998）。此外，海洋地震调查成本高昂、具有间歇性，并可能影响海洋生态系统（Mitra等人，2023；Lumley和Behrens，1998；McCauley等人，2000；Gordon等人，2003；Nowacek等人，2013；Pirotta等人，2014）。基于钻孔的技术，如脉冲中子测井（Müller等人，2007）、交叉孔地震（Ajo-Franklin等人，2013）、带有分布式声学感应的垂直地震剖面（VSP）（Correa等人，2019）、交叉孔电阻率层析成像（ERT）（Schmidt-Hattenberger等人，2016）及钻孔重力测量（Alyousuf等人，2024），虽然具有较高的垂直分辨率，但水平覆盖范围有限，仅限于井筒或井间走廊，因此提供的空间信息有限。

在这种情况下，我们提出了μ成像技术，也称为宇宙射线μ子射影或层析成像（Tanaka等人，2023），作为一种补充的地下监测技术。μ成像通过分析宇宙射线μ子穿过感兴趣体积时的衰减（吸收）或多次库仑散射来获取图像（Bonechi等人，2020）。与钻孔方法不同，单个μ子探测器提供了宽的立体视角和储层规模的水平覆盖范围。原位空间分辨率首先由探测器的角度分箱和目标距离（Δx≈LΔθ）决定，在典型的储层条件下，每个分箱可达到数十米的分辨率（Bonechi等人，2020；Tanaka等人，2014；Oláh等人，2019）。在吸收μ成像中，检测到的μ子通量变化直接与μ子路径上的积分密度相关。由于GeV级μ子能够穿过数百米的岩石，μ成像提供了其他地球物理方法难以实现的成像能力。在过去的十年中，这项技术已在多种情境下得到研究，包括火山学（Tanaka，2019）、考古学（例如金字塔的内部成像Alvarez等人，1970；Morishima等人，2017）、土木工程（Guardincerri等人，2017；Oláh等人，2023）、边境控制（Barnes等人，2023）、文化遗产研究（Giammanco等人，2025）以及磁场监测（Basiri等人，2022；Kin等人，2023）。最近的实验室规模实验表明，μ成像能够检测多孔介质中流体饱和度变化引起的微妙整体密度变化（Eleslambouly等人，2025；Eleslambouly等人，2024）。

**2. 探测器技术和分析方法的最新进展**
探测器技术和分析方法的进步提高了深层地下μ成像的可行性。紧凑型便携式探测器（Chaiwongkhot等人，2022）和适用于钻孔的检测系统（D’Errico等人，2023；Thompson等人，2024；IDEON Technologies，2025）解决了高压高温环境中的部署挑战。机器学习算法可以从低μ子计数数据中提取有意义的信号（Galgóczi等人，2024；Lefevre等人，2024；O’Donnell等人，2025；Paccagnella等人，2025）。这些进展促使人们将μ成像与地震方法结合使用，作为互补工具（Matsushima等人，2024）。

**1.1. 本研究的先验工作和动机**
一些开创性研究探讨了μ成像在CCS监测中的应用，证明了概念的可行性，并 identification 了定量应用的关键挑战。Kudryavtsev等人（2012）首次提出使用MUSIC（μ子模拟代码）/MUSUN（地下μ子模拟）蒙特卡洛代码可以检测深层咸水含水层中的二氧化碳羽流。Klinger等人（2015）使用GEANT4模拟了776米深度处的储层密度变化，并在不到50天内检测到羽流。Thompson等人（2022）在布尔比矿（Boulby Mine）的恶劣地下环境中开发并测试了基于闪烁体的探测器原型。Steer等人（2024）在斯韦尔维克CO2现场实验室（挪威）首次进行了现场部署，使用直径60毫米的钻孔探测器在38-65米深度处成功检测到二氧化碳羽流，结果与地震测量结果非常吻合（Steer等人，2024；Jordan等人，2022）。尽管取得了这些进展，但之前的研究主要属于定性研究，重点在于能否检测到羽流，而非如何系统地区分不同的饱和度水平（Matsushima等人，2025）。它们没有量化μ子传输如何随离散的二氧化碳饱和度水平（例如25%、50%、100%）变化，也没有研究中等饱和度是否会产生统计上可区分的特征。此外，这些研究采用了简化的单层羽流几何结构，如均匀圆柱体（Kudryavtsev等人，2012）或径向扩展的羽流（Klinger等人，2015），而不是实际CCS项目（如斯莱普纳）观察到的多层结构（Furre等人，2017；Cavanagh和Haszeldine，2014）。

计算和实际限制进一步限制了定量评估。早期的模拟使用了仅限于均匀单材料传输的MUSIC/MUSUN（Kudryavtsev等人，2012）。这些简化模型无法表示复杂的异质地质结构或真实的探测器几何形状。Klinger等人（2015）通过使用GEANT4模拟了多层上覆岩层，但他们的二氧化碳羽流几何结构仍然是单层径向扩展的。这些模拟研究使用了假设的大面积探测器阵列（1000平方米）进行概念验证（Zhong等人，2015；Thompson等人，2022），并承认实际钻孔探测器受到井筒尺寸的限制。

**1.2. 斯莱普纳CCS项目和研究目标**
斯莱普纳CCS项目位于北海距岸约250公里处（图1）。自1996年以来，每年大约有100万吨二氧化碳被注入高孔隙度的乌蒂拉沙层，到2020年已成功储存超过2000万吨（Furre等人，2017；Izadian，2024）。乌蒂拉沙层是一个主要由未固结沙组成的厚储层单元，具有高孔隙度和渗透性。注入的二氧化碳被夹在页岩盖层下，起到上部区域密封的作用，防止向上泄漏。我们选择斯莱普纳作为测试地点，是因为它处于CCS μ成像近乎最佳的深度窗口内。首先，大规模地质二氧化碳储存要求注入的二氧化碳处于超临界状态（临界点：Tc=31.1°C，Pc=73.8巴），以实现足够的储存效率；在典型的地热梯度下，这种超临界条件大约在800米深度达到（Bui等人，2018）。这为可行的CCS储层深度设定了热力学下限。其次，μ子穿过上覆岩层的衰减随深度急剧增加，为μ成像的可行性设定了实际上限：超过约1000米后，实际探测器的采集时间变得不可行。斯莱普纳的储层深度（海平面下800-1000米）位于这两个限制的交叉点（接近超临界二氧化碳储存物理可行的最浅深度），并且在μ成像的实际范围内。相比之下，In Salah（约1800米，White等人，2014）和Quest（约2000米，Rock等人，2017）的深度远远超出μ成像的实际范围，μ子衰减使计数率相对于斯莱普纳深度降低了几个数量级。在斯莱普纳的较浅深度，减少的μ子衰减使得在实际采集时间内能够积累足够的统计数据（Crozat等人，2022）：使用6平方米的探测器，在500米深度采集需要10天，在1000米深度需要282天，在1500米深度需要2470天，这促使我们将CCS μ成像的重点放在1000米以下的深度（Matsushima等人，2025）。对于本研究而言，最重要的是，斯莱普纳具有通过重复的4D地震调查清晰描绘的九层羽流结构（Cavanagh和Haszeldine，2014；Chadwick等人，2005），使得能够系统评估多层敏感性，这是之前的单层羽流研究无法评估的。值得注意的是，斯莱普纳的4D地震监测数据集可通过CO2 DataShare存储库公开访问（CO2 DataShare，2024），使本研究中采用的羽流几何参数任何读者都可以独立验证。这种开放访问特性也使斯莱普纳成为未来研究的理想测试平台：两阶段的MUSIC→PHITS工作流程可以直接应用于从公开存档的地震解释中获得的更详细的流动模拟派生的饱和度场。

**2. 本研究的主要成果**
1. **定量饱和度区分**：我们系统地评估了μ成像是否能够区分25%、50%和100%的二氧化碳饱和度水平，并使用统计假设检验量化在理想化模拟条件下实现3σ分离所需的采集时间。
2. **高效的两阶段计算工作流程**：我们将用于千米级上覆岩层传输的MUSIC与用于高分辨率储层建模的PHITS（粒子与重离子传输代码系统）相结合，实现了近三个数量级的计算效率提升（约700倍）。该工作流程使得参数化研究公里级几何结构成为可能，否则在计算上不可行。
3. **特定地点的多层二氧化碳羽流几何结构**：我们使用成熟的地震解释（Cavanagh和Haszeldine，2014；Chadwick等人，2005）对斯莱普纳的九层二氧化碳羽流进行了建模，采用了现实的深度依赖密度分布和盖层相互作用，使地质真实性超越了理想化的单层羽流模型。

总之，之前的CCS-μ成像研究回答了μ成像是否能够检测到羽流的定性问题，而本研究则回答了我们是否能够区分不同饱和度水平的定量问题。从检测到区分的转变，使μ成像成为连续、定量CCS储层监测的潜在补充工具。两阶段架构将计算时间相对单阶段全传输方法减少了近三个数量级（约700倍）：第一阶段计算一次并在所有场景中重复使用；只有第二阶段在每个饱和情况下都会重复计算。2. 方法我们采用两阶段蒙特卡洛工作流来模拟μ子通过Sleipner地质环境的传输（图2概述）。MUSIC（Kudryavtsev，2009）将μ子从海平面传输到储层深度，然后PHITS（粒子与重离子传输代码系统）（Sato等人，2024）将它们传输通过CO2羽流层到达探测器平面。尽管PHITS本身可以模拟μ子的产生和传播，但我们的初步测试显示由于模拟尺寸较大（直径为1.6公里的圆柱体和高度为1公里），导致计算时间过长。这种工作流将计算时间减少了近三个数量级（约700倍），同时保持了物理上的准确性。为了量化这一点，我们使用相同的硬件（Apple M4处理器，单核）对1000个能量范围在500-2500 GeV之间的μ子通过圆柱形岩体积（半径800米，厚度700米）进行了基准测试。PHITS运行使用了随PHITS包提供的推荐μ子传输输入文件（phits/recommendation/muon/muon.inp）作为配置参考。MUSIC大约在2秒内完成了等效的传输，而PHITS则需要1400秒。这两种代码都随μ子数量线性扩展，因此这个比例与样本大小无关，并且在生产模拟中使用的10^7 μ子规模下仍然成立。2.1. 地质和羽流模型我们采用了简化的圆柱几何形状来表示Sleipner CCS站点的地质环境，其中注入发生在海平面下大约1012米的深度（Eiken，2019）。模拟域是一个半径为800米、总高度约为1050米的垂直圆柱体。从顶部到底部，圆柱体包括80米的海水、700米的岩石覆盖层，以及Utsira地层中的9个离散CO2羽流层（Eiken，2019）。PHITS中的几何设置如图3(a)所示；完整的单元定义、表面卡和坐标规范在补充材料S1中提供。所定义的几何形状保留了μ子相互作用的物理特性，同时去除了不必要的结构细节。这些简化减少了计算复杂性，使得千米级的μ子传输模拟能够在实际时间内完成，从而使得模拟适合进行可行性研究。主要的物理过程，如μ子通过电离和辐射过程的能量损失，以及通过羽流层的路径长度分布，在这种几何形状中得到了保留。岩石覆盖层使用了McConn等人（2011）提供的标准岩石元素组成（见补充材料，表S1）。九个CO2羽流层的尺寸、深度和属性是根据2002年7月的报告地质数据获得的，当时大约有5百万吨CO2被注入Sleipner（Cavanagh和Haszeldine，2014）。每个羽流层被近似为一个圆柱区域，其体积、厚度、温度、压力和CO2密度在表1中总结，并在图3(b)中示意性地展示。下载：下载高分辨率图像（530KB）下载：下载全尺寸图像图3. 基于Sleipner的模拟域的PHITS几何实现。(a) 完整的圆柱域（半径800米，高度约1050米）从PHITS输入文件中渲染出来，从海平面延伸到储层下方；完整的单元定义和表面卡在补充材料S1中提供。(b) 九层CO2羽流结构，近似于实际的Sleipner羽流几何形状。表1. 基于Cavanagh和Haszeldine（2014）的CO2羽流层总结。层深度厚度体积压力温度CO2密度空单元(mbsl*) (m) (10^6 m^3) (MPa) (°C) (kg/m^3)空单元顶端平均空单元空单元空单元空单元空单元空单元L1 95 5960 10.5 2.7 89.7 139.5 616L2 920 926 12.0 4.2 99.3 738.6 606L3 90 590 99.0 2.5 19.2 038.1 602L4 890 895 9.5 1.18 9.0 637.7 597L5 875 880 13.0 10.4 08.9 137.3 589L6 860 866 10.5 1.29 8.7 736.9 582L7 845 849 10.0 2.6 08.6 036.4 571L8 825 832 14.5 3.34 8.4 35.9 551L9 795 80 110.0 1.4 08.1 235.0 355mbsl* = 海平面以下米。2.1.1. 层平均饱和度参数化我们规定每个羽流层的CO2饱和度在空间上是均匀的，并考虑了三种情况（25%、50%、100%）。通常，层平均饱和度定义为空间内层面积的平均值，(1)S?CO2=∫_Alayer SCO2(x,y)dA，其中Alayer是羽流层的水平横截面积。在我们的模拟中，SCO2(x,y)=SCO2在每层内是空间恒定的，因此方程（1）简化为分配的值本身。这相当于通过方程（2）为每个圆柱形PHITS单元指定一个体积密度，其中SCO2∈{0.25,0.50,1.00}在给定场景的所有九层中均匀应用。这种理想化使得可以在明确定义的饱和状态下系统地比较μ子的衰减。实际上，由于浮力驱动的迁移、毛细捕获和地质异质性，CO2饱和度在每层内是空间变化的。然而，对于严格垂直穿过整个羽流堆的μ子，测量的传输只对储层厚度上累积的柱平均密度敏感。在保持这个柱平均值不变的情况下重新分配各层之间的饱和度，在我们的垂直入射近似下将产生基本相同的μ子传输。因此，我们基于场景的方法系统地改变了总路径积分密度（μography直接测量的量），使得这种参数化适用于建立μography对由CO2注入引起的体密度变化的物理敏感性。我们注意到，即使在像Sleipner这样被广泛监测的站点，定量三维CO2饱和度分布也不是直接观察到的，而是通过反推4D地震振幅数据并受到储层流动模拟的限制（Cavanagh和Haszeldine，2014，Furre等人，2017）推断出来的；在中等到高饱和度下，地震振幅-饱和度关系的非唯一性（Landr?，2001）意味着没有一个独立的真实异质饱和度场可以用来进行基准测试。因此，我们基于场景的方法提供了一个透明的、可重复的路径积分密度对比参数化，这是μography直接测量的。这种理想化对实际羽流条件的影响在第4.4节中讨论。2.1.2. 体积密度计算为了估计不同CO2饱和度场景下每个羽流层的密度，我们使用了孔隙率加权的密度方程：(2)ρ_tot=(1??)ρ_grain+?SCO2ρ_CO2+(1?SCO2)ρ_water，其中：•?=0.37是砂岩层的孔隙率，•ρ_grain=2.65 g/cm^3是砂岩颗粒的密度，•ρ_CO2按层变化（见表1），•ρ_water=1.0 g/cm^3是水的密度，•SCO2表示在方程（1）中定义的层平均饱和度。图4展示了25%、50%和100% CO2饱和度场景下计算出的层密度，突出了各羽流层之间的密度对比。下载：下载高分辨率图像（148KB）下载：下载全尺寸图像图4. 25%、50%和100% CO2饱和度场景下九个羽流层的估计密度剖面。2.1.3. 宇宙射线μ子谱和能量阈值为了模拟真实的宇宙射线μ子源，我们采用了CosmoALEPH合作组织（Schmelling等人，2013）在海平面测量的垂直μ子谱和电荷比。ALEPH探测器位于大约120米的标准岩层下方，记录了动量范围在100 GeV/c到2.5 TeV/c内的垂直μ子通量；这些测量结果已经校正了回海平面的能量损失。我们在模拟中使用了他们校正后的差分通量Φ_total(p)和电荷比R(p)=?+(p)/??(p)。为了将总通量分为正μ子和负μ子的贡献，我们为?+(p)和??(p)在实验动量范围内创建插值函数。图5显示了测量得到的通量以及一个单幂律拟合，(3)Φ_fit(p)=10^3.85?3.20log10(p)，用于扩展超出表格范围的采样。为了避免传播无法穿透覆盖层的μ子，我们设定了一个较低的动量截止值。Utsira Sand储层位于大约80米的海水和700米的岩石下方，相当于总共730米的标准岩层覆盖。根据Particle Data Group的μ子范围表（Groom等人，2001），μ子需要大约700 GeV的动能才能穿过这个厚度（图6）。我们采用了500 GeV的保守截止值，以考虑岩石密度和深度的不确定性，确保所有有现实机会到达探测器的μ子都被包括在内。超过这个阈值的μ子占总入射通量的3.78%。下载：下载高分辨率图像（272KB）下载：下载全尺寸图像图5. 基于CosmoALEPH数据的垂直μ子通量及其统计和系统误差条以及幂律拟合（方程（3））。下载：下载高分辨率图像（421KB）下载：下载全尺寸图像图6. 标准岩层中的μ子能量与动量关系。超过706 GeV的μ子能穿透730米的岩石；模拟的截止值为500 GeV。2.1.4. μ子源生成程序对于超过500 GeV阈值的μ子，我们通过积分通量来构建累积分布函数（CDF）并对其进行逆向处理以实现随机动量采样。每个μ子根据（4）Prob(μ+)=R(p)1+R(p)，Prob(μ?)=11+R(p)来分配电荷，其中R(p)是动量p处的电荷比。生成程序如下进行：1. 从CosmoALEPH数据获取Φ_total(p)和R(p)。2. 将总通量分解为?+和??。3. 对?±(p)进行插值，或者使用幂律拟合以减少插值开销。4. 构建?++??的CDF并通过CDF逆向采样动量。5. 根据R(p)在每个采样的动量处分配电荷。图7显示了结果的能量和空间分布，并确认生成的分布与CosmoALEPH输入谱一致，验证了实现的正确性。下载：下载高分辨率图像（861KB）下载：下载全尺寸图像图7. 在海平面生成的宇宙射线μ子的初始属性。(a) 空间源分布。(b) 用于覆盖层传输的能谱。2.2. MUSIC通过覆盖层的传输在模拟工作流的第一阶段，10^7个垂直μ子（θ=0°）在海平面生成，并使用MUSIC（MUon SImulation Code）蒙特卡洛代码（Kudryavtsev，2009）通过覆盖层进行传输。MUSIC通过模拟电离能量损失（连续减速）、制动辐射、对产生和光核相互作用以及多次库仑散射引起的角度偏转来传播μ子。完整的传输设置和材料定义在补充材料S2中提供。限制为垂直入射有三个理由：垂直μ子主导了海平面通量（微分通量遵循cos2θ，在θ=0°时达到峰值（Gaisser等人，2016）；它们穿过的覆盖层路径长度最小，因此在储层深度的存活概率最高；在垂直近似下，每个水平羽流层的路径长度等于其真实厚度，使得柱密度积分在分析上变得可行，密度重建明确。这种理想化是特意为这项可行性研究采用的，以隔离μography对CO2引起的密度变化的物理敏感性，而不增加角度反卷积的复杂性。实际上，接近垂直的入射可以在实验中实现：多层跟踪探测器可以重建单个μ子轨迹，并允许离线角度选择|θ|≤θ_max，使得这种理想化在隧道部署或表面阵列配置中可以直接实现（Bonechi等人，2020，D’Errico等人，2023，Thompson等人，2024）。这种限制的影响在第4.4节中讨论。由于MUSIC一次只传输一种均匀材料，因此覆盖层模拟分两个步骤进行。首先，10^7个μ子通过80米的海水（ρ=1.00 g/cm3；柱密度8.0×10^3 g cm?2）传输，这产生的衰减可以忽略不计（<0.1%），但通过电离和辐射过程导致平均能量下降了大约43 GeV（从807 GeV降至764 GeV）。然后，存活的μ子通过700米的标准岩石（ρ=2.65 g/cm3；柱密度1.855×10^5 g cm?2）传输，这产生了主要的衰减。在出口平面（780米深度），48.8%的μ子存活下来。这些存活的μ子作为PHITS通过储层层的传输的源输入。图8显示了存活μ子的MUSIC后分布。能量谱（图8a）显示了显著的衰减：平均能量从807 GeV降至177 GeV，中位数从676 GeV降至103 GeV，这是由于电离（dE/dx≈2 MeV cm?2）和辐射过程（dE/dx∝bE，b≈4×10?? cm2 g?1），后者在临界能量Ec≈500 GeV以上占主导（Groom等人，2001）。天顶角分布（图8b）保持紧密聚焦（θrms=0.82°，与Highland公式（Highland，1975）对700米标准岩石在177 GeV时的预测≈0.70°一致；超过99%的μ子在|θ|<5°范围内），验证了随后储层规模传输的接近垂直近似。MUSIC采用基于时间的随机数种子，由于第一阶段计算一次并在所有饱和场景中重复使用其输出，这不会影响场景间比较的可重复性。下载：下载高分辨率图像（288KB）下载：下载全尺寸图像图8. 780米深度的MUSIC后μ子分布（第一阶段存活比例为初始的48.8%）。(a) 能量谱：平均值177.3 GeV，中位数103.5 GeV，显示了从输入平均值807.5 GeV的显著衰减。(b) 天顶角：θrms=0.82°，超过99%的μ子在|θ|<5°范围内，验证了接近垂直的近似。图表放大以强调接近垂直的峰值。2.3.水库规模的传输和探测器建模（PHITS）在模拟工作流程的第二阶段，PHITS将存活的μ子传输通过九个离散的CO2羽流层，并在位于最深羽流层（第9层，平均深度801米）下方约50米的虚拟探测器平面上记录撞击事件，总深度约为海平面下1021米（PHITS表面坐标z=?102,145厘米）。模拟分别针对25%、50%和100%的层平均CO2饱和度进行，并与无羽流的基线情况进行了比较，在无羽流情况下，所有九个羽流层都被赋予了砂岩的密度。完整的PHITS输入设置，包括粒子传输参数、能量截止值、物理模型、源配置和计数规范，都在补充材料S1中提供。为了在保持运行时间可控的同时最大化统计精度，探测器被建模为一个完美的高效圆形虚拟平面，半径为800米，直接位于羽流堆下方（图9）。μ子撞击事件使用PHITS中的T-Cross计数器记录，该计数器记录穿过表面的所有粒子，并且只计算初级μ子（排除次级粒子）。探测平面被划分为20×20的空间网格，每个网格覆盖80×80平方米的面积。这种配置是理想化的，意味着探测平面没有材料厚度、死区或效率损失，并且能够拦截比任何实际钻孔或隧道基础探测器阵列更多的μ子。这种便利性允许在不产生过高计算成本的情况下区分不同的饱和度情况，并能够清晰地评估μ成像对CO2引起的密度变化的敏感性。这种理想化配置与实际探测器部署（包括将采集时间缩放到有限探测器面积、立体角接受度和效率）之间的关系在第2.5节中进行了推导。我们模拟了107个μ子历史记录，分为100个批次，每个批次包含105个μ子。这种批次结构有两个目的：它提供了增量检查点，以便在中断时可以从最后一个完成的批次继续模拟，而无需丢弃全部107个μ子的运行数据；并且它产生了每个批次的探测器撞击总数{Nk}k=1100，这些总数直接用作累积Z分数分析（第2.4.2节）的样本序列，避免了任何事后采样。所有PHITS运行都使用了一个固定的初始种子（rseed = 3），每个后续批次都是从上一个批次的最终随机数状态确定性地种子化的。因为所有饱和度场景都读取相同的MUSIC输出文件并共享相同的种子序列，所以相同的107个μ子历史记录（具有相同的入射坐标、方向和能量）会通过每个羽流配置。因此，不同批次之间的计数差异仅反映了CO2密度变化的物理效应。MUSIC和PHITS阶段的全随机数生成器规范，包括RANLUX（James, 1994）和RANMAR（Marsaglia和Zaman, 1991）算法，都在补充材料S2中提供。下载：下载高分辨率图像（134KB）下载：下载全尺寸图像图9. 在PHITS模拟中使用的理想化探测平面。垂直线表示代表性的μ子轨迹；括号内的堆表示九个CO2羽流层。20×20像素网格（400像素，每个80×80平方米）用于传输映射和密度重建。2.4. 数据处理和统计分析分析分为两个独立的部分。第一部分将原始像素撞击计数转换为平滑的传输图、差分传输图像和带有相关不确定性的路径平均密度重建（第2.4.1节）。第二部分直接使用原始批次计数来通过累积Z分数评估不同饱和度场景之间的统计可分离性（第2.4.2节）。2.4.1. 传输、平滑和密度重建传输定义。μ子传输比率按逐个像素的基础定义如下：(5)T(x,y)=Nscenario(x,y)/Nsource(x,y)，其中Nsource(x,y)是在MUSIC源平面（海平面，500–2500 GeV）生成的初始宇宙射线μ子的计数，Nscenario(x,y)是在给定储层情景下到达探测器平面（像素(x,y)的μ子计数。因此，传输代表了从海表面通过上覆层（MUSIC阶段：80米海水＋700米岩石）和储层层（PHITS阶段）到探测器的总路径衰减。所有后续处理步骤都使用这种逐像素传输作为主要观测值。高斯平滑。原始传输图包含由有限统计量引起的像素级泊松噪声。我们对每个20×20的传输图应用二维高斯平滑，使用以下核函数：(6)G(i,j)=1/(2πσ^2)exp(-i^2+j^2/2σ^2)，其中σ=1.5个像素间距（≈120米），规范化使得∑i,jG(i,j)=1，并在3σ（≈5个像素）处截断。平滑后的传输图通过离散卷积获得：(7)Tsmooth(x,y)=∑i=??3σ??3σ?∑j=??3σ??3σ?Traw(x+i,y+j)?G(i,j)。在像素边界处，应用最近邻外推（scipy.ndimage.gaussian_filter中的模式=‘nearest’），这复制了最近的有效像素值，而不是在探测器边缘引入人为的零值。这一步将外部像素噪声从大约10%减少到大约3%，同时保持了中心羽流的特征（第3.2节）。差分传输。为了将羽流引起的信号与空间均匀的背景衰减区分开，我们从每个平滑的羽流场景中减去平滑的无羽流基线：(8)ΔT(x,y)=Tplume(σ)(x,y)?Tbase(σ)(x,y)，其中正值（ΔT>0）表示相对于无羽流情况的μ子传输增加，反映了由于超临界CO2替代了更密集的地层盐水而导致的路径积分密度降低。密度重建。根据Minato的衰减公式（Minato, 1986），路径平均密度与传输的关系如下：(9)ρ=?ln(T/κL)。在我们的实现中，T被归一化到无羽流基线而不是初始源，这消除了共同的上覆层衰减，并引入了一个密度偏移ρ0：(10)ρij=ρ0?ln[Tplume(i,j)/T0(i,j)]κL，其中ρij是像素(i,j)中的路径平均体密度；ρ0=2.65克/厘米^3是宿主砂岩的基线体密度；TPLume(i,j)和T0(i,j)分别是羽流和无羽流情况下的平滑传输图；L=185米是μ子垂直穿过的九层羽流堆的总厚度；κ=3.9×10^?6厘米^2克^1是质量衰减系数。因为基线归一化Tplume/T0消除了共同的上覆层衰减（80米海水＋700米岩石），L仅代表羽流堆的路径长度，而不是从海平面到探测器的总距离（约1021米）。κ的值是通过对比无羽流基线经验验证的：将公式(10)应用于平滑的无羽流图证实κ=3.9×10^?6厘米^2克^1可以再现已知的砂岩密度ρ0=2.65克/厘米^3，在中心像素(ic,jc)=(10,10)处——选择该点是因为其信噪比最高（最大垂直路径长度，最小边缘效应）在0.1%以内。校准是全局的（单个κ均匀应用于所有像素和场景）；其残余不确定性被吸收到下面的系统术语中。因为校准使用了无羽流基线，所以羽流和参考情况共有的系统偏差在相对传输比率中大致抵消了。相对密度亏损为：(11)Δρ=ρij?ρ0/ρ0，其中正值表示由于存在低密度的超临界CO2而导致路径平均密度降低。重建密度的不确定性来自两个独立的来源。来自泊松计数的统计不确定性为：(12)σρ,stat(i,j)=1/κL√Nplume(i,j)+1/Nref(i,j)，来自残余校准和路径长度误差的系统不确定性为：(13)σρ,sys(i,j)=|ρij|Δκ/κ^2+ΔL/L^2，其中Δκ/κ=0.05是5%的残余校准不确定性（Minato (1986)对异质地质介质的典型±15%分布的保守比例），ΔL/L=0.02反映了由于表1中的层深度不确定性而产生的有效路径长度的2%不确定性。总不确定性为：(14)σρ,tot(i,j)=σρ,stat^2+σρ,sys^2，通过公式(11)传播以获得Δρ/ρ0的置信区间。2.4.2. 统计显著性测试为了独立于任何传输或密度反演来量化不同饱和度场景之间的统计可分离性，我们直接从原始批次计数计算累积双样本Z分数。对于每对场景A和B，N批次后的累积Z分数为：(15)Z(N)=∑k=1^N(NA,k?NB,k)∑k=1^N(σA,k^2+σB,k^2)，其中σA,k和σB,k是批次k的泊松计数不确定性。在高斯近似下对泊松统计量的处理是合理的，因为中心羽流区域每个批次的每像素计数超过50（平均每个像素每个批次约240次计数），满足标准正态近似条件。使用Z≥2σ和Z≥3σ的阈值来识别区分不同饱和度水平所需的最少μ子历史记录数量。3σ阈值对应于一尾p值约为0.0013；2σ阈值对应于p≈0.023。这些等价关系在物理文献中是标准的。2.5. 采集时间扩展到实际探测器配置本研究中报告的采集时间基于一个理想化的模拟构建：严格的垂直μ子入射（θ=0°）和半径为800米的完美高效虚拟计数表面（Aideal≈2.0×10^6平方米）。这个大面积不代表一个可部署的仪器；选择它是为了在模拟中高效地积累大量数据。模拟μ子计数与实际世界采集时间之间的对应关系取决于角接受度标准。使用半径为800米的圆形探测表面和CosmoALEPH在500 GeV以上的垂直μ子谱（Schmelling等人，2013年），我们计算了不同半角截止θmax的入射μ子率（表2）。我们采用θmax=2°作为参考，因为它保留了垂直μ子的近似（在图8b中得到验证），同时提供了31.0 s^-1的μ子率，相当于大约3.8天可以积累1×10^7个检测到的μ子（即每105个μ子的批次约54分钟）。为了将模拟计数与实际世界探测器的接受度相关联，我们计算了θmax=2°圆锥所张开的立体角ΔΩideal：(16)ΔΩideal=2π(1?cos^2θmax)≈3.83×10^?3sr。检测到的μ子率与有效探测器面积Adet、立体角接受度ΔΩdet和效率??成线性比例。因此，实际探测器实现与理想化配置相同统计精度所需的采集时间为：(17)?real=tideal×Aideal/Adet×ΔΩideal/ΔΩdet×1?。公式(17)提供了一个关于曝光时间的一阶近似。虽然与面积（A）和效率（??）的缩放是严格的线性关系，但与立体角接受度（ΔΩ）的逆向缩放假设探测器锥体内的μ子通量是均匀的。实际上，在储层深度，偏离轴线的μ子穿过上覆层的斜向深度（h/cosθ）显著更长，因此被不成比例地衰减。因此，从狭窄的垂直圆锥（θ≤2°）线性外推到宽角度接受度会高估真实的宽角度计数率，这意味着从这个方程得出的treal是一个乐观的下限。表2. 半径为800米的圆形表面和500 GeV能量阈值的μ子率及所需采集时间，这些数据来自CosmoALEPH（Schmelling等人，2013年）测量的垂直μ子通量。参考情况θmax=2°（粗体）保留了垂直入射的近似。对于θmax≥5°的比率是从立体角缩放∝(1?cosθmax)外推的，代表了上限，因为斜向深度的衰减在780米深度时减少了净增益（见第4.4节）。尽管如此，公式(17)中的主导因素是面积比Aideal/Adet：单个钻孔探测器的Adet约为1–10平方米，比理想化计数表面小五到六个数量级。因此，实际部署要么需要延长的积分周期，要么需要探测器阵列，或者两者兼而有之。Crozat等人（2022年）对6平方米探测器的独立深度缩放计算提供了一个实际参考：在500米深度大约需要10天，在1000米深度需要约282天，在1500米深度需要约2470天，以达到3σ的检测置信度，这说明了探测器面积、角接受度和储层深度共同决定了采集的可行性。3. 结果本节展示了从107个μ子模拟中提取的关键观测值，并量化了区分不同CO2饱和度的统计能力以及检测羽流区域的可行性。3.1. 传输图逐像素的传输比率（公式(5)）是主要的观测值。图10显示了三种CO2饱和度情景下探测器深度的传输图。在羽流核心区域（r?400米）内，传输增强效果清晰可见，随着饱和度从25%增加到100%，传输比率呈单调增加。第3.1.1节分析了核心区域的定量平均值及其空间依赖性。信号在r≈400米之外的径向衰减反映了羽流边界外未受扰动的储层密度的过渡，这些图为第3.2节中的差分传输和密度重建提供了定性预览。3.1.1.1.### 定量信号分析
为了量化羽流探测信号的空间依赖性，我们定义了三个同心分析区域（如图11a所示）：
- **内层区域（r≤200米）**：最大信号区域，对应于沿近似垂直轨迹的最高积分CO2柱密度。
- **核心区域（r≤400米）**：以羽流为主的区域，与名义上的羽流足迹半径相匹配。该区域提供了主要的检测指标，将真实的羽流特征与背景稀释区分开来。
- **整个探测器覆盖区域（r≤800米）**：包括羽流核心（占面积的25%）和非羽流背景（占面积的75%）在内的整个圆形孔径。在该区域上进行平均仅用于参考，以说明当忽略空间选择性时信号会如何稀释；它不代表推荐的分析策略。

表3总结了每个区域的空间平均μ子计数和饱和度水平。核心区域平均（r≤400米）显示出明显的单调缩放关系：25%饱和度时增强率为（5.50±0.11）%，50%时为（5.97±0.11）%，100%时为（6.79±0.12%）。内层区域（r≤200米）在100%饱和度时显示出高达（10.74±0.15%）的峰值对比度。作为参考，整个探测器平均——均匀地对待所有像素，因此没有利用已知的羽流足迹几何形状——将100%饱和度的信号稀释为（2.15±0.06）%，这是由于包含了75%没有CO2扰动的背景像素所导致的。实际应用中，通过选择感兴趣的空间区域（相当于核心区域平均）可以利用已知的羽流足迹来恢复完整的信号对比度。图11展示了探测信号的空间分布和径向依赖性。图11(a)显示了100%饱和度的传输图，并叠加了三个分析区域，直观地证实了增强传输集中在羽流足迹内。图11(b)量化了信号随径向距离的衰减：所有三种饱和度情景都在r≈100米处显示出峰值增强，随着μ子轨迹离开CO2饱和体积，信号逐渐降低到基线水平。随着饱和度水平的增加（25% < 50% < 100%），这种单调缩放关系在所有径向位置上都成立，证明了CO2的定量敏感性。一条水平参考线表示100%饱和度的核心区域平均值（+6.79%）。整个探测器平均（+2.15%，虚线）仅用于参考，以说明当忽略空间信息时信号会稀释三倍；具有空间分辨率的实时探测器会应用感兴趣区域的选择来利用已知的羽流足迹几何形状。

### 地图3. 三个分析区域的空间平均μ子计数
核心区域（r≤400米）的值将羽流特征与背景稀释区分开来，并作为主要的检测指标。整个探测器值仅用于参考，以说明在未应用空间选择性时的稀释情况。不确定性代表了通过空间平均传播的泊松计数统计。

| 情景 | 整个探测器 | 核心区域 | 内层区域 | 空单元 |
|------------------|---------|---------|----------|---------|
| 无羽流 | 10,231 | 11,070 | | |
| 25% CO2 | 10,413 | 11,386 | +5.50% | |
| 50% CO2 | 10,424 | 11,119 | +5.97% | |
| 100% CO2 | 10,451 | 11,205 | +6.79% | |

**下载：下载高分辨率图像（593KB）**
**下载：下载全尺寸图像**

### 图11. 羽流探测信号的空间分布和径向依赖性
(a) 增强的μ子传输集中在羽流足迹内（r≤400米）。
(b) 随着CO2饱和度的增加，所有径向位置都表现出单调缩放；核心区域（+6.79%）与整个探测器（+2.15%）平均值之间的三倍差异表明了忽略空间选择性时的信号稀释。所有结果：107个μ子历史数据，在理想化的垂直入射条件下。

### 3.2. 差分传输图
为了将羽流引起的信号与空间均匀的背景衰减和像素级统计波动区分开来，每个原始传输图都与高斯核（σ=1.5个bin ≈120米）进行卷积，并减去平滑后的无羽流基线（ΔT，公式（8））。正值（ΔT>0）表示μ子传输增加，反映了由于CO2存在而导致的路径积分密度降低，因为CO2替代了更密集的流体。

图12显示了三种饱和度情景下的ΔT图。正值清楚地勾勒出了羽流足迹。随着饱和度的增加，羽流层重叠较少的区域（图3(b)）变得越来越可区分，显示出更强的局部密度对比。峰值从25%饱和度的ΔT≈0.05增加到100%饱和度的ΔT≈0.08，证实了CO2的定量敏感性，并在空间上恢复了Utsira层的九层堆叠几何结构。

**下载：下载高分辨率图像（322KB）**
**下载：下载全尺寸图像**

### 3.3. 密度重建
通过公式（10）从差分传输图重建路径平均密度，其中ρ0=2.65 g cm?3，κ=3.9×10?6 cm2 g?1，L=185米。图13显示了三种情景下的百分比密度亏损图（Δρ/ρ0）。每张图都显示出一个随饱和度单调增加的中心峰值亏损，与CO2逐渐取代更密集流体的过程一致，并在空间上恢复了九层堆叠的羽流几何结构。表4量化了羽流核心区域（r≤400米）的重建路径平均密度，并考虑了传播的统计和系统不确定性（公式（12）、（13）、（14））。100%与50%饱和度之间的密度对比（Δρ=0.019 g cm?3）超过了2σtot=0.014 g cm?3，使用107个μ子获得了统计上显著的区分性。50%与25%之间的对比（Δρ=0.011 g cm?3）位于2σtot范围内，与第3.4节中量化的部分分布重叠一致。基于密度的分离（2.0σ和1.1σ）比Z分数结果（≥3σ，第2.4.2节）要弱，因为两种测试使用了不同的统计池：Z分数在整个探测器孔径上积分了原始的总计数（约3.3×10^6个μ子，无系统误差项），而表4的不确定性是从受δκ/κ和δL/L（公式（12）、（13）、（14）影响的核心区域像素平均值得出的。两个结果在内部是一致的，并且相互补充：Z分数量化了全局区分性，而密度对比量化了羽流核心中的空间信号强度。

**下载：下载高分辨率图像（302KB）**
**下载：下载全尺寸图像**

### 表4. 羽流核心区域的路径平均密度（r≤400米）
统计不确定性（σstat≈0.004 g cm?3）主要受固定数量（107个μ子样本）的影响。系统不确定性来源于衰减系数校准（δκ/κ5%）和羽流厚度不确定性（δL/L=2%），通过公式（12）–（14）传播。95%置信区间为±2σtot≈±0.014 g cm?3。

| 沉饱度 | ρ | σstat | σsys | σtot | 95%置信区间 |
|-------|--------|--------|------|---------|
| 100% | 2.487 | 0.004 | 0.006 | ±0.014 |
| 50% | 2.506 | 0.004 | 0.006 | ±0.014 |
| 25% | 2.517 | 0.004 | 0.006 | ±0.014 |

**密度对比：Δρ**
| | | | | |
| 100% vs 50% | -0.019 | --- | --- |
| 50% vs 25% | -0.011 | --- | --- |

### 3.4. 全局计数统计和区分性
表5总结了整个孔径（r≤800米）上的积分μ子计数。相对泊松不确定性1/N≈5.5×10^-4与情景间的差异1.74–2.12%相比可以忽略不计，验证了后续所有测试中使用的高斯近似。图14显示了这些积分计数及其1σ误差条。100%情景与所有较低饱和度情景明显分开，而25%和50%情景有部分重叠。

表5. 所有情景下的积分探测器计数（r≤800米）。泊松不确定性σ=N。相对变化是相对于无羽流基线计算的。

| 情景 | 总计数 | σ（泊松） | 相对变化 |
|------------|---------|---------|---------|
| 无羽流 | 3,294,308 | 181 |
| 25% CO2 | 3,351,691 | 183 |
| 50% CO2 | 3,355,269 | 183 |
| 100% CO2 | 3,364,021 | 183 |

重叠积分OAB=∫fmin(n)dn得出O100,50=0.017，O50,25=0.329，O100,25=0.001，分别对应于1.7%、33%和0.1%的误分类概率。图15可视化了成对的高斯分布。单个107个μ子的测量（在理想化配置下约4天）在98.3%的情况下正确区分了100%和50%的饱和度。为了达到95%的置信度（O50,25<0.05），我们需要将Z分数分离从当前的Zcur≈1.95σ（由O50,25=0.329通过2Φ(?Z/2)=0.329得出）增加到Ztarget≈3.92σ（满足2Φ(?Z/2)=0.05的值）。由于Z∝N，所需的μ子计数按（Ztarget/Zcur）^2=(3.92/1.95)^2约为16天，在相同的理想化假设下。

### 4. 讨论
本研究将CCS微告age从定性的羽流检测推进到了定量的饱和度区分。在理想化的模拟条件下（仅垂直入射的μ子、效率完美的800米半径虚拟探测器以及均匀的层平均饱和度），μ子传输随CO2饱和度单调增加，路径平均密度图恢复了Utsira层的九层堆叠几何结构。100%与50%的情景在大约4天的实时时间内以≥3σ的显著性区分开；区分50%与25%则需要更长的曝光时间，这反映了较小的密度对比（表4）。这些发现表明，在第2.1.4节详细阐述的理想化假设条件下，吸收微告age能够以物理上可解释且可量化的方式响应孔隙流体的替代。

### 4.2. 统计性能和实际扩展性
由于累积Z分数随Z∝N（公式（15）增长，曝光时间、探测器面积和立体角接受度可以通过公式（17）中的缩放关系进行互换。在理想化配置下，100%与50%的比较在大约7.8×10^6个μ子后达到3σ显著性，而100%与25%的对比在相似时间范围内得到解决。区分50%与25%在107个μ子后仍然低于2σ；要将O50,25的误分类阈值（95%置信度）从0.329降低到0.05，需要大约4倍的计数，即在相同的理想化假设下约16天。对于实际的钻孔探测器，限制性能的主要实际因素按重要性顺序为：(i) 有效的探测器面积和立体角接受度，它们直接影响计数率；(ii) 随深度变化的垂直μ子通量衰减，这显著增加了获取时间（Crozat等，2022）；(iii) 大气压力和温度变化在小时到天的时间尺度上调节μ子通量（Gaisser等，2016；Barrett等，1952）；以及(iv) 地质异质性和与探测器相关的噪声源（第4.4节讨论）。

### 4.3. CCS监控中的微告age
除了物理敏感性之外，选择CCS监控策略还必须明确考虑部署成本、后勤复杂性以及可实现的空间分辨率。高分辨率的地球物理方法，如4D地震和4D VSP，属于高成本领域，需要大量的操作努力——重复采集、专门的团队和严格的基线控制——然而它们提供了对羽流几何结构和应力诱导的速度变化的稳健成像（Arts等，2004；Chadwick等，2009；Lumley和Behrens，1998）。斯莱普纳（Sleipner）的实地经验表明，时间延迟地震（time-lapse seismic）技术可以描绘羽流（plume）的演变过程并限制二氧化碳（CO2）的质量分布，但这需要较高的运营成本和财务投入（Arts等人，2004年；Chadwick等人，2009年）。相比之下，地球化学监测、井口监测和示踪技术相对成本较低且操作流程较为简单，然而它们主要对泄漏路径或间接的储层响应敏感，而非对体积饱和度的变化敏感（Benson和Cole，2008年；Ringrose和Meckel，2019年）。在这几种方法中，μ射线成像（muography）在成本效益（图17a）和成本操作复杂性（图17b）方面都处于中间位置。其对于由二氧化碳替代引起的总体密度差异的敏感性提供了对质量重新分布的直接限制，提供了一种基于密度的观测指标，这种指标不受到基于速度的方法中常见的饱和度地质学复杂性的影响（Kudryavtsev等人，2012年；Gluyas等人，2019年）。从操作角度来看，μ射线成像需要部署探测器并持续较长时间进行数据采集，但避免了重复的主动源观测活动，因此相对于4D地震来说具有中等的操作要求（Gluyas等人，2019年）。尽管其空间分辨率通常低于反射地震学（由探测器几何形状和角度接收范围而非波长决定，Bonechi等人，2020年），但它提供的体积密度信息超过了基于点采样的方法的诊断能力（Kudryavtsev等人，2012年）。因此，μ射线成像最好被视为综合监测组合中的一个补充工具，特别是在中等至高饱和度条件下，此时地震振幅-饱和度关系变得不清晰且无法进行诊断（Landr?，2001年；Knight等人，1998年），而地球化学方法无法解析体积变化（Benson和Cole，2008年；Ringrose和Meckel，2019年）。

为了对这些方法进行半定量比较：在斯莱普纳进行的重复海底重力测量每年可检测到大约10^8公斤的二氧化碳质量变化，横向分辨率为约200-500米（Alnes等人，2008年；Alnes等人，2011年；Nooner等人，2007年）；4D地震可以在数十米的横向尺度上解析羽流的几何形状，但在二氧化碳饱和度超过20-30%时分辨率降低且非唯一性增加（Landr?，2001年；Arts等人，2004年）；示踪和地球化学方法只能提供点检测，无法解析羽流的几何形状；而本研究中的μ射线成像在185米的羽流范围内实现了平均路径密度敏感度约为Δρ≈0.06克/立方厘米（ΔT≈4×10^-3，对应ΔSCO2=0.5），其空间分辨率主要由探测器几何形状决定。因此，在中等至高饱和度条件下，当地震响应饱和且不再具有诊断能力时，μ射线成像具有独特的优势（Landr?，2001年；Knight等人，1998年），而地球化学方法无法解析体积变化（Benson和Cole，2008年；Ringrose和Meckel，2019年）。为了将这些比较放在半量化的基础上：在斯莱普纳，重复的海底重力测量每年可检测到大约10^8公斤的二氧化碳质量变化，横向分辨率为约200-500米；4D地震可以在数十米的横向尺度上解析羽流的几何形状，但在二氧化碳饱和度超过20-30%时分辨率降低且非唯一性增加；示踪和地球化学方法只能提供点检测，无法解析羽流的几何形状；而本研究中的μ射线成像在185米的羽流范围内实现了平均路径密度敏感度约为Δρ≈0.06克/立方厘米（ΔT≈4×10^-3，对应ΔSCO2=0.5），其空间分辨率主要由探测器几何形状决定。因此，在中等至高饱和度条件下，当地震响应饱和且不再具有诊断能力时，μ射线成像具有独特的优势，它能够提供分层的密度信息，而地球化学方法无法提供这种信息。

4.4. 限制性因素与实地部署路径
这项可行性研究存在一些固有的限制，在实际部署之前需要在未来的工作中加以解决。仅限于接近垂直方向的μ子（θ≤2°）的假设忽略了偏轴通量的可用性；如表2所示，如果将角度接受范围从θmax=2°扩展到3°，表面μ子率大约会翻倍（从31 s^-1增加到69 s^-1），尽管这是一个乐观的上限，因为偏轴μ子会穿过倾斜的覆盖层路径h/cosθ，并且在780米深度时衰减更为严重。将实际的角度分布纳入正向模型将增加计数率并减少采集时间，但这需要路径长度依赖的传输校正，这超出了本可行性研究的范围。相反，如果从操作上希望仅使用垂直入射的μ子，可以通过多层跟踪探测器在硬件中实现这一要求，这些探测器能够重建和过滤单个μ子的轨迹（Bonechi等人，2020年；Thompson等人，2024年）。理想的800米半径虚拟探测器（Aideal≈2×10^6平方米）的面积比任何可部署的钻孔阵列都要大几个数量级；如公式（17）所示，一个10平方米的钻孔探测器如果效率达到80%，且实际的角度接受范围约为1弧秒（sr），在这种深度下需要大约10年的时间才能完成采集。需要注意的是，这种线性的立体角缩放提供了一个乐观的下限，因为它忽略了偏轴μ子在更大倾斜深度下的严重衰减（h/cosθ）。这说明了优化阵列几何形状和探测器技术的重要性。将相同的饱和度分配给所有九层并不符合典型的碳捕获与封存（CCS）羽流演变（Cavanagh和Haszeldine，2014年；Furre等人，2017年）；然而，由于μ子的衰减与路径积分的密度有关，只要不同情景之间的总密度差异保持不变，这里建立的单调传输-饱和度关系仍然成立。这里的假设忽略了覆盖层和储层中的地质不均匀性，这可能会在传输图中引入空间变化的偏移，这些偏移可能模仿羽流引起的信号；通过注入前的基线采集并结合4D地震数据进行联合反演是将地质背景与时间延迟的流体替换信号分离的主要策略。大气通量的变化（数小时到数天内约0.1%-1%）可以通过连续的压力记录结合气压校正（Barrett等人，1952年）、地面参考探测器和时间平均来减轻。探测器噪声、电子死时间和对准不确定性超出了这项物理敏感性研究的范围，但在任何实地试验之前必须在硬件层面进行表征。

尽管如此，深部μ射线成像的实际可行性正在迅速提高。包括Ideon Technologies（IDEON Technologies，2025年）和Geoptic（Thompson等人，2024年；D’Errico等人，2023年）在内的公司已经开发出了与钻孔兼容的探测器，后者首先在STFC的Boulby地下实验室进行了演示，并随后在挪威的一个浅层CCS现场试验中部署（Steer等人，2024年；Jordan等人，2022年），这些探测器的外径小于90毫米，与标准钻孔套管兼容，其压力-温度等级适合储层深度。这些探测器的有效面积比这里假设的理想化值Aideal≈2×10^6平方米小几个数量级；如公式（17）所示，对于斯莱普纳深度100%与50%的区分度，这样的探测器在理想情况下需要数年甚至数十年的采集时间。斯莱普纳现有的基础设施，包括观测井15/9-A-23（Chadwick等人，2005年；Falcon-Suarez等人，2018年），在满足访问和完整性要求的条件下，可以为吸收型μ射线成像提供部署几何结构，而无需钻专门井的成本。这里建立的MUSIC→PHITS工作流程非常适合在硬件投入之前对这种几何结构进行预期信号的正向建模，并与4D地震数据结合进行联合反演，以利用互补的敏感性。推荐的下一步包括：（i）将实际的角度μ子分布和有限的钻孔-探测器几何形状纳入模拟；（ii）使用斯莱普纳基准模型对异质、基于流动模拟的饱和度场进行建模（Cavanagh和Haszeldine，2014年；Furre等人，2017年）；（iii）通过注入前的基线采集来表征大气和地质噪声源；（iv）评估在≤1000米深度下平衡空间分辨率与采集时间的多探测器阵列配置。

5. 结论
这项理想化的可行性研究评估了宇宙射线μ射线成像对二氧化碳饱和度变化的物理敏感性，模型基于斯莱普纳CCS站点，采用了两阶段蒙特卡洛工作流程（MUSIC用于覆盖层，PHITS用于羽流层）。重要的是要区分这里展示的基于场景的区分能力与完整的饱和度估计：可观测的量是沿垂直μ子轨迹的路径积分的总体密度差异，而对空间变化的二氧化碳饱和度分布进行连续的定量反演超出了当前理想化几何和技术假设的范围，需要额外的角度多样性、探测器密度和经过验证的正向模型。μ射线成像在这个背景下的基本物理优势在于它直接对由于孔隙-流体替换引起的总体密度差异敏感，这一量随着二氧化碳饱和度的变化而单调变化，即使在地震振幅-饱和度关系变得不清晰且无法进行诊断的中等至高饱和度条件下也保持线性诊断能力（Landr?，2001年；Landr?和Zumberge，2017年）；这使得μ射线成像成为时间延迟4D地震监测的物理补充方法，而不是竞争方法。所有定量发现都基于三个理想化的模拟条件：仅限垂直入射的μ子、一个效率完美的800米半径虚拟探测器（Aideal≈2×10^6平方米），以及对所有九层施加均匀的层平均饱和度。这些结果不应被解释为可部署系统性能的预测。

μ子传输随着二氧化碳饱和度的增加而单调增加。在上述条件下，100%和50%的饱和度情景可以实现≥3σ的统计分离。要在相同的理想化假设下区分50%和25%，大约需要4倍的μ子数量。这种敏感性体现在路径积分的总体密度差异层面；精确的逐层饱和度量化超出了当前配置的分辨率。基于斯莱普纳地震解释得出的九个离散羽流层（Cavanagh和Haszeldine，2014年；Chadwick等人，2005年）实施了具有特定层的深度、厚度、压力、温度和二氧化碳密度，使得地质真实性超过了以往研究的均匀单层模型。路径平均的密度图恢复了Utsira地层的叠加层几何形状，证明了吸收型μ射线成像能够以物理上可解释的方式响应孔隙-流体替换。这些图中的明显空间划分不应被解释为实地部署中的实际成像分辨率：这是由施加的圆柱形羽流几何形状和平滑核（σ=1.5个bin，约120米）所驱动的，而不是由实际探测器的固有角度分辨率所驱动的。

MUSIC→PHITS架构相对于单阶段的全传输计算时间减少了近三个数量级。第一阶段（MUSIC）计算一次，然后在第二阶段（PHITS）中重复使用780米深度下的存活μ子群体。这种效率提升使得在其他计算成本禁止的情况下也能对公里级几何结构进行参数研究，并为未来包含异质羽流几何形状或实际探测器响应函数的 ?al??malar 提供了一个可重用的框架。对于代表性的钻孔探测器（A≈10^-2–10^-1平方米），这意味着每个探测器的采集时间根据饱和度比较可能需要几个月到几年的时间，强调了多钻孔阵列和优化的井布置是实际现场应用的先决条件，而不是根本的物理障碍。总的来说，这些结果确立了μ射线成像对斯莱普纳相关深度处孔隙-流体替换引起的总体密度差异具有物理敏感性，MUSIC→PHITS框架为未来的实验设计提供了计算效率高的、可重用的基础。该方法的被动、非侵入式采集方式及其对总体密度的直接敏感性（而不是声波速度）使其成为在4D地震振幅响应变得不清晰的饱和度范围内的补充监测工具。在第4.3节讨论的主要实际限制和实地部署前的具体下一步将在4.4节中详细讨论。最关键的未解决问题是实际探测器噪声和钻孔探测器的小检测面积、大气通量的变化以及不均匀的饱和度分布。推荐的下一步包括：（i）将实际的角度μ子分布和有限的钻孔-探测器几何形状纳入模拟；（ii）使用斯莱普纳基准模型对异质、基于流动模拟的饱和度场进行建模；（iii）通过注入前的基线采集来表征大气和地质噪声源；（iv）评估在≤1000米深度下平衡空间分辨率与采集时间的多探测器阵列配置。

在已建立的CCS测量、报告和验证（MRV）框架内，包括欧盟指令2009/31/EC（欧洲议会和理事会，2009年）和ISO27914:2017（国际标准化组织，2017年），μ射线成像在技术上与所有三个监测阶段都兼容。在注入前的基线采集期间，μ射线成像调查提供了未充气储层的独立密度图，所有注入后的传输变化都可以据此进行归一化，这是无论后续监测频率如何通常都需要的一个步骤。在主动注入期间，其被动和连续的特性使得在时间分辨率上优于间歇性的4D地震活动，特别是在上述中等至高饱和度条件下具有最高的区分能力。在关闭后的验证期间，永久部署的钻孔阵列消除了重复的主动源动员成本，并提供了羽流稳定性的连续证据。第3.5节中确定的定量可检测性阈值为实现特定场地MRV计划中的μ图监测要求提供了物理基础；将这些阈值转化为监管合规标准被确定为未来工作的课题。

作者贡献声明：
Hamid Basiri：撰写 – 审稿与编辑、撰写 – 初始草稿、可视化、验证、软件开发、方法论、调查、形式分析、概念化。
Jun Matsushima：撰写 – 审稿与编辑、项目监督、方法论、资金筹集、概念化。
Masashi Kodama：撰写 – 审稿与编辑。
Toshiyuki Yokota：撰写 – 审稿与编辑。
Mohammed Y. Ali：撰写 – 审稿与编辑、概念化。
Fateh Bouchaala：撰写 – 审稿与编辑。
Ahmed Eleslambouly：撰写 – 审稿与编辑、可视化。
Hiroyuki K.M. Tanaka：撰写 – 审稿与编辑、项目监督、资金筹集。
Eduardo Cortina Gil：撰写 – 审稿与编辑、方法论。
Andrea Giammanco：撰写 – 审稿与编辑、验证、方法论、概念化。